COMSOL 6.4 - Hysteresis in Piezoelectric Ceramics

Many piezoelectric materials are ferroelectric. Ferroelectric materials exhibit nonlinear polarization behavior such as hysteresis and saturation at large applied electric fields. In addition, the polarization and mechanical deformations in such materials can be strongly coupled due to the electrostriction effect. This model uses the Ferroelectroelasticity interface to analyze a simple actuator made of PZT piezoelectric ceramic material, which is subjected to applied electric field and mechanical load.

Model Definition

The direct electrostrictive effect for a material of arbitrary symmetry can be represented as the following additive contribution to the strain:

which is quadratic in polarization P. Due to the symmetry, the fourth order tensor Q can be effectively represented by a 6-by-6 coupling matrix. For piezoelectric ceramics, the matrix can be characterized by three independent components: Q11, Q12, and Q44.

For ferroelectroelastic materials, the polarization vector is nonlinear function of the electric field and possible mechanical stress in the material. The Jiles–Atherton model is available in COMSOL Multiphysics for modeling ferroelectric hysteresis. It assumes that the total polarization can be represented as a sum of reversible and irreversible parts. The polarization change is computed from the following incremental equation:

where the reversibility is characterized by the parameter cr, and the anhysteretic polarization is found from a relation:

where Ps is the saturation polarization. The polarization shape is characterized by the Langevin function

where a is a material parameter called the domain wall density.

The effective electric field is given by

(1)

where E is the applied electric field, α is a material parameter called the inter-domain coupling, and the mechanics stress is computed assuming mechanically linear material as

where C is the fourth order elasticity tensor. The last term in Equation 1 represents the inverse electrostrictive effect.

Finally, the change of the irreversible polarization is computed from the following incremental relation:

where the pinning loss is characterized by the parameter kp.

The ferroelectroelastic actuator in this model example is a rectangular plate with dimensions of 1.5 in-by-0.25 in-by-0.015 in, which is composed of PZT-5H piezoelectric ceramic material. The following polarization parameter values have been estimated in Ref. 1 based on experimental data:

Table 1: Material properties of PZT-5H.
Material property	Value	Description
Ps	0.425 C/m2	Saturation polarization
a	6.4105 V/m	Domain wall density
α	4.2·106 m/F	Inter-domain coupling
cr	0.2	Polarization reversibility
kp	1·106 V/m	Pinning loss

The mechanical properties for PZT-5H are available in the Material Library of COMSOL.

The coupling coefficients for PZT ceramics can vary with the material composition and temperature. The reference values used in this example are given in the table below (Ref. 2.

Table 2: Electrostrictive coupling coefficients.
Material property	Value
Q11	3.579·10-2 m4/C2
Q12	-5.33510-3 m4/C2
Q44	1.923·10-2 m4/C2

The upper surface of the actuator is grounded, while the lower one is subjected to an electric potential that can cyclically vary in small increments between −Vmax and +Vmax.

The actuator can be subjected to a compressive stress by applying boundary loads of various magnitude.

Because of the symmetry, it is sufficient to model one quarter of the actual geometry.

Figure 1: Model geometry.

Because of the large aspect ratio of the actuator and the unidirectional nature of the electrical and mechanical loading, a coarse mesh can be used for the discretization.

Figure 2: Model mesh.

Results and Discussion

Three full cycles have been computed for each value of Vmax. The variation of polarization and electrostrictive strain is studied at the point in the middle of the actuator. The first cycle includes the initial transient; see Figure 3 and Figure 4. The hysteresis loops become fully established after two full cycles; see Figure 5 and Figure 6.

Finally, Figure 7 and Figure 8 show the effect of the applied compressive stress.

Figure 3: Polarization hysteresis loop including the initial transient for the maximum applied voltage of 600 V.

Figure 4: Electrostrictive strain hysteresis loop including the initial transient for the maximum applied voltage of 600 V.

Figure 5: Polarization hysteresis loops fully established after two initial cycles for different values of the maximum applied voltage.

Figure 6: Electrostrictive strain hysteresis loops fully established after two initial cycles for different values of the maximum applied voltage.

Figure 7: Fully established polarization hysteresis loops for different values of the mechanical load and maximum applied voltage of 1200 V.

Figure 8: Fully established strain hysteresis loops for different values of the mechanical load and maximum applied voltage of 1200 V.

Notes About the COMSOL Implementation

In this example, you study the hysteresis with respect to the incremental variation of the applied electric potential using a stationary parametric study. The same hysteresis model can be also used for time dependent studies.

References

1. R.C. Smith and Z. Ounaies. “A Domain Wall Model for Hysteresis in Piezoelectric Materials,” J. Int. Mat. Sys. Struct., vol. 11, no. 1, pp. 62–79, 2000.

2. B. Völker, P. Marton, C. Elsässer, and M. Kamlah, “Multiscale modeling for ferroelectric materials: a transition from the atomic level to phase-field modeling,” Continuum Mech. Thermodyn., vol. 23, pp. 435–451, 2011.

MEMS_Module/Piezoelectric_Devices/piezoelectric_hysteresis

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

Define parameters for the geometry, material properties, applied voltage, and mechanical load.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
t	0[s]	0 s	Time parameter
W	1.5[in]	0.0381 m	Actuator width
D	0.25[in]	0.00635 m	Actuator depth
H	0.015[in]	3.81E-4 m	Actuator height
alpha	4.2e6[m/F]	4.2E6 m/F	Interdomain coupling
a	6.4e5[V/m]	6.4E5 V/m	Domain wall density
c	0.2	0.2	Polarization reversibility
k	1e6[V/m]	1E6 V/m	Pinning loss
Ps	0.425[C/m^2]	0.425 C/m²	Saturation polarization
Q11	3.579e-2[m^4/C^2]	0.03579 m4/C²	Electrostriction coupling parameter
Q12	-5.335e-3[m^4/C^2]	-0.005335 m4/C²	Electrostriction coupling parameter
Q44	1.923e-2[m^4/C^2]	0.01923 m4/C²	Electrostriction coupling parameter
Vmax	1200[V]	1200 V	Maximum applied voltage
F0	0[MPa]	0 Pa	Applied mechanical load

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
V0	Vmaxsin(2pi*t[1/s])	V	Applied voltage

This variation of the potential with respect to the parameter at one of the actuator boundaries will cause the electric field within the material to gradually change between -Vmax and Vmax.

Geometry 1

Block 1 (blk1)

In the toolbar, click

Because of the symmetry, it is sufficient to model one quarter of the actuator.

In the window for , locate the section.

In the text field, type W/2.

In the text field, type D/2.

In the text field, type H/2.

Click

Solid Mechanics (solid)

Linear Elastic Material 1

Prescribe the material stiffness using data available in the material library for PZT-5H, which is represented by the full elasticity matrix.

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Electrostatics (es)

Charge Conservation, Ferroelectric 1

In the window, under > click .

In the window for , locate the section.

Select the checkbox.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Right-click and choose .

In the toolbar, click

to close the window.

Materials

Lead Zirconate Titanate (PZT-5H) (mat1)

Define the remaining ferroelectric properties for the material using the parameters.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Saturation polarization	Psat	Ps	C/m²	Ferroelectric
Interdomain coupling	alphaJAe_iso ; alphaJAeii = alphaJAe_iso, alphaJAeij = 0	alpha	m/F	Ferroelectric
Domain wall density	aJAe_iso ; aJAeii = aJAe_iso, aJAeij = 0	a	V/m	Ferroelectric
Pinning loss	kJAe_iso ; kJAeii = kJAe_iso, kJAeij = 0	k	V/m	Ferroelectric
Polarization reversibility	cJAe_iso ; cJAeii = cJAe_iso, cJAeij = 0	c	1	Ferroelectric