PDF

Hysteresis in Piezoelectric Ceramics
Introduction
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:
Ps
6.4105 V/m
α
4.2·106 m/F
cr
kp
1·106 V/m
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 give in the table below (Ref. 2):
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.
Application Library path: MEMS_Module/Piezoelectric_Devices/piezoelectric_hysteresis
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  3D.
2
In the Select Physics tree, select Structural Mechanics>Electromagnetics-Structure Interaction>Ferroelectroelasticity.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Stationary.
6
Global Definitions
Parameters 1
Define parameters for the geometry, material properties, applied voltage, and mechanical load.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Definitions
Variables 1
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
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)
1
In the Geometry toolbar, click  Block.
Because of the symmetry, it is sufficient to model one quarter of the actuator.
2
In the Settings window for Block, locate the Size and Shape section.
3
In the Width text field, type W/2.
4
In the Depth text field, type D/2.
5
In the Height text field, type H/2.
6
Click  Build All Objects.
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.
1
In the Model Builder window, under Component 1 (comp1)>Solid Mechanics (solid) click Linear Elastic Material 1.
2
In the Settings window for Linear Elastic Material, locate the Linear Elastic Material section.
3
From the Solid model list, choose Anisotropic.
4
From the Material data ordering list, choose Voigt (11, 22, 33, 23, 13, 12).
Electrostatics (es)
Charge Conservation, Ferroelectric 1
1
In the Model Builder window, under Component 1 (comp1)>Electrostatics (es) click Charge Conservation, Ferroelectric 1.
2
In the Settings window for Charge Conservation, Ferroelectric, locate the Ferroelectric Material Properties section.
3
Select the Hysteresis JilesAtherton model check box.
Add Material
1
In the Home toolbar, click  Add Material to open the Add Material window.
2
Go to the Add Material window.
3
In the tree, select Built-in>Lead Zirconate Titanate (PZT-5H).
4
Right-click and choose Add to Component 1 (comp1).
5
In the Home toolbar, click  Add Material to close the Add Material window.
Materials
Lead Zirconate Titanate (PZT-5H) (mat1)
Define the remaining ferroelectric properties for the material using the parameters.
1
In the Settings window for Material, locate the Material Contents section.
2
Solid Mechanics (solid)
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Boundary Load 1
1
In the Physics toolbar, click  Boundaries and choose Boundary Load.
2
3
In the Settings window for Boundary Load, locate the Force section.
4
Specify the FA vector as
Electrostatics (es)
In the Model Builder window, under Component 1 (comp1) click Electrostatics (es).
Ground 1
1
In the Physics toolbar, click  Boundaries and choose Ground.
2
Electric Potential 1
1
In the Physics toolbar, click  Boundaries and choose Electric Potential.
2
Because of the symmetry, the voltage at the horizontal symmetry plane is half of that applied at the bottom surface.
3
In the Settings window for Electric Potential, locate the Electric Potential section.
4
In the V0 text field, type V0/2.
Multiphysics
Electrostriction 1 (efe1)
Study the electrostriction in the material using a fully coupled model.
1
In the Model Builder window, under Component 1 (comp1)>Multiphysics click Electrostriction 1 (efe1).
2
In the Settings window for Electrostriction, locate the Coupling Type section.
3
From the list, choose Fully coupled.
Because of a certain symmetry in the material microstructure, you need three parameters to characterize the electrostrictive coupling.
4
Locate the Electrostriction section. From the Solid model list, choose Cubic crystal.
5
In the Q11 text field, type Q11.
6
In the Q12 text field, type Q12.
7
In the Q44 text field, type Q44.
Mesh 1
Mapped 1
1
In the Mesh toolbar, click  Boundary and choose Mapped.
2
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
From the Predefined list, choose Coarse.
Swept 1
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 1 and choose Distribution.
2
In the Settings window for Distribution, locate the Distribution section.
3
In the Number of elements text field, type 2.
4
Click  Build All.
Study 1
Step 1: Stationary
In the first study, no mechanical load is assumed, so the entire excitation is caused by the applied electric field.
1
In the Model Builder window, under Study 1 click Step 1: Stationary.
2
In the Settings window for Stationary, click to expand the Study Extensions section.
3
Select the Auxiliary sweep check box.
4
Compute three full cycles for the applied electric potential for each given maximum value.
5
6
Click the  Show More Options button in the Model Builder toolbar.
7
In the Show More Options dialog box, select Study>Batch and Cluster in the tree.
8
9
Batch Sweep
1
In the Study toolbar, click  Batch and choose Batch Sweep.
2
In the Settings window for Batch Sweep, locate the Study Settings section.
3
4
5
Locate the Batch Settings section. Find the Before sweep subsection. Clear the Clear meshes check box.
6
Clear the Clear solutions check box.
7
Locate the Advanced Settings section. In the Number of simultaneous jobs text field, type 3.
Batch Data
In the Study toolbar, click  Compute.
Results
Polarization
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Polarization in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1/Parametric Solutions 1 (sol2).
4
Locate the Grid section. Select the Manual spacing check box.
5
In the y spacing text field, type 0.1.
6
Locate the Legend section. From the Position list, choose Lower right.
Point Graph 1
1
Right-click Polarization and choose Point Graph.
2
3
In the Settings window for Point Graph, locate the y-Axis Data section.
4
In the Expression text field, type es.PZ.
5
Locate the x-Axis Data section. From the Parameter list, choose Expression.
6
In the Expression text field, type es.EZ.
7
From the Unit list, choose MV/m.
8
Click to expand the Legends section. Find the Include subsection. Clear the Point check box.
9
Select the Show legends check box.
10
In the Polarization toolbar, click  Plot.
Electrostriction
1
In the Model Builder window, right-click Polarization and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Electrostriction in the Label text field.
3
Locate the Grid section. In the y spacing text field, type 0.001.
Point Graph 1
1
In the Model Builder window, expand the Electrostriction node, then click Point Graph 1.
2
In the Settings window for Point Graph, locate the y-Axis Data section.
3
In the Expression text field, type efe1.emZZ.
4
In the Electrostriction toolbar, click  Plot.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies.
Add one more stationary study to analyze the mechanical load effect.
4
In the Select Study tree, select General Studies>Stationary.
5
6
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Stationary
1
In the Settings window for Stationary, locate the Study Extensions section.
2
Select the Auxiliary sweep check box.
3
4
Batch Sweep
1
In the Study toolbar, click  Batch and choose Batch Sweep.
2
In the Settings window for Batch Sweep, locate the Study Settings section.
3
4
5
Locate the Batch Settings section. Find the Before sweep subsection. Clear the Clear meshes check box.
6
Clear the Clear solutions check box.
7
Locate the Advanced Settings section. In the Number of simultaneous jobs text field, type 3.
Batch Data
In the Study toolbar, click  Compute.
Results
Polarization 1
1
In the Model Builder window, right-click Polarization and choose Duplicate.
2
In the Settings window for 1D Plot Group, locate the Data section.
3
From the Dataset list, choose Study 2/Parametric Solutions 2 (sol7).
4
In the Polarization 1 toolbar, click  Plot.
Strain
1
In the Model Builder window, right-click Electrostriction and choose Duplicate.
2
In the Settings window for 1D Plot Group, type Strain in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2/Parametric Solutions 2 (sol7).
Point Graph 1
1
In the Model Builder window, expand the Strain node, then click Point Graph 1.
2
In the Settings window for Point Graph, locate the y-Axis Data section.
3
In the Expression text field, type solid.eZZ.
4
In the Strain toolbar, click  Plot.