DOI : 10.17577/IJERTV8IS100144
- Open Access
- Authors : Rabenjarivelo R. P. H , Ranaivoarison Z. N , Andrianaharison Y
- Paper ID : IJERTV8IS100144
- Volume & Issue : Volume 08, Issue 10 (October 2019)
- Published (First Online): 23-10-2019
- ISSN (Online) : 2278-0181
- Publisher Name : IJERT
- License:
This work is licensed under a Creative Commons Attribution 4.0 International License
Finite Element Modelling of Ferroelectric using Quadratic Mesh
Rabenjarivelo R. P. H
PhD. Student University of Antananarivo in STII
Antananarivo, Madagascar
Ranaivoarison Z. N
Ph.D
University of Antananarivo, Antananarivo, Madagascar
Andrianaharison Y
Full Professor
University of Antananarivo in ESPA Antananarivo, Madagascar
AbstractPiezoelectric and ferroelectric active materials are widely used in many fields of technology and science. Using such materials requires numerical tools taking into account this coupling. This work focus on the formulation of a solid finite element that integrates a constitutive law for ferroelectric and ferroelastic piezoceramics. It has been implemented into the Matlab code. A linear and quadratic mesh structure is presented and compared in this paper to model these materials using the finite element method.
KeywordsPiezoelectric materials, ferroelectric, finite element method, hysteresis.
-
INTRODUCTION
Active materials such as Ceramic Lead Zirconate Titanate (PZT) play an increasing role in actuator and sensor applications in intelligent structures. It is important to know the phenomenology of these materials for example: coupling, non- linearity and hysteresis to model, optimize, design and development of smart devices. It is why several constitutive models, are considered in different aspects of the actual behavior of ferroelectric ceramics have been developed. In addition to that phenomenological models is too important to understand the macroscopic behavior on the basis of microscopic properties.
The purpose of this article is to use the structure using quadratic mesh. The choice of the meshes is an essential question, the finer the mesh, the less it contributes to the differences between modeling and reality, but also the cost of calculation increases. This method achieves to the improvement of the distribution of the stresses, electric displacement field and a decrease of the average stress in the interesting part, however it leads to the appearance of a too important displacement at the junction between the zone of design and the non-design area near the recess.
Switching of an electric domain may occur when a ferroelectric is under the action of an applied electric field. This process is not continuous, but only occurs when the strength of the applied electric field exceeds certain critical value, causing the polarization intensity P lag behind the applied electric field
-
It is generally agreed that such nonlinear relation of hysteresis results from the switching of electric domain. When the electric field changes periodically, electric hysteresis loop
is formed between P and E, which is an important sign of ferroelectricity. A typical electric hysteresis loop of a ferroelectric is shown in Fig.1.
Fig.1. Hysteresis loop.
: Electric coercive,
: Polarization remnant,
: Spontaneous polarization.
Fig. 1. Typical hysteresis loop of ferroelectric ceramics
Many methods can be used to have hysteresis loop of ferroelectric material. Here, we introduce volumic fraction in the finite element code.
-
-
PIEZOELECTRIC EQUATIONS
-
State equation
Piezoelectricity is the result of a coupling between the electrical energy and the mechanical energy of a material neglecting the thermal energy and the pyroelectric effect. The fundamental quantities of piezoelectricity are for the electrical part: electrical field E and electrical displacement D, and for the mechanical part: stress and the strain .
These four quantities can be described as the partial derivatives of the Gibbs thermodynamic potential:
=
(1)
: Gibbs thermodynamic potential,
U: Internal energy.
For an adiabatic and reversible transformation, we have:
=
(2)
From where
=
(3)
If we develop the Gibbs potential in the first order, we can write:
= ( ) + ( )
(4)
Which by identification gives:
= ( ) = ( )
(5)
Thus, if we express the total derivatives of stress and induction (electrical displacement), we obtain:
where the exponents l and s respectively correspond to the reversible and remnant or spontaneous parts. The linear parts can be described by the piezoelectric equation:
= = ( )
(15)
= + = = ( ) +
(16)
-
Piezoelectric coefficients:
-
= =
(17)
= =
(18)
= =
(19)
= =
(20)
= =
(17)
= =
(18)
= =
(19)
= =
(20)
The preceding equations can determine the relations between the piezoelectric coefficients.
2 2
= ( ) ( )
(6)
2 2
= ( ) ( )
(7)
2
= ( )
(8)
2
= ( )
(9)
2
= ( )
(10)
2 2
= ( ) ( )
(6)
2 2
= ( ) ( )
(7)
2
= ( )
(8)
2
= ( )
(9)
2
= ( )
(10)
With
, = 1 6 and , = 1 Ã 3
where :
: Piezoelectric coefficient,
: Compliance,
: Permittivity Electric. From where
= +
(11)
= +
(12)
The same development can be applied to express other pairs of independent variables. We obtain a linear system of 8 tensor equation given in table 1.
TABLE I. DIMENSIONS OF TOP AND BOTTOM PIEZOELECTRIC
Variable
Electrical part
Mechanicalpart
,
= +
= +
,
=
= +
,
= +
=
,
=
=
C. Finite element method
-
Basic equation:
Basic equations governing the electro-elastic continuity behavior of piezoelectricity are presented, which are the starting point for finite element formulation. The virtual work and the relationships of the energetic formulation are then established. Some specifications related to the problems of the modeling of intelligent structures, such as: the electromechanical coupling and induced potential are discussed.
The electroelastic response of a piezoelectric body of volume and regular boundary surface S, is governed by equations of mechanical, dynamic and electrostatic:
, + =
(21)
, = 0
(22)
Where , , are components of the mechanical force, electric body charge and mass density, respectively , and
, are the Cauchy symmetric stress tensor and electric displacement vector components.
They are connected to those of the symmetric linear Lagrange tensor and the electric field vector . Equations of linear piezoelectric (direct effect) are given by:
=
(23)
= +
(24)
= ( )
(25)
= ( ) + +
(26)
= ( )
(25)
= ( ) + +
(26)
Ferroelectric equations are:
Due to the directional dependence of material properties (anisotropy), parameters of the constituent equations of piezoelectricity must be represented by tensor.
B. Non-linear piezoelectric equation (ferroelectric)
-
Constitutive law for ferroelectric and ferroelastic ceramics:
In this part, we summarize the key elements of constitutive law proposed by M.Kamlah, S.C. Hwang, C. S. Lynch [1,2 et
, , denote the elastic, piezoelectric and dielectric constants of the material.
, are the remnant of the deformation and the polarization.
The strain tensor and the vector components of the electric field are related to displacement components and the electric potential field , by the following relationships,
1
= 2 (, + ,)
(27)
= ,
(28)
1
= 2 (, + ,)
(27)
= ,
(28)
3], the stress and the electrical displacement are given as
follows:
= +
(13)
= +
(14)
1
= 2
(38)
1
= 2
1
2 + +
(39)
1
= 2
(38)
1
= 2
1
2 + +
(39)
The piezoelectric body , could be subject to either essential or natural mechanical and electric boundary conditions, or a combination of them, on its boundary ,
and are respectively the kinetic energy and the extension of the potential energy including the electrical contribution, defined by the following expressions:
where 0
= 0 = in
(29)
= 0 = , in
(30)
= 0 = in
(29)
= 0 = , in
(30)
is the prescribed stress on and 0
is the
prescribed displacement taken to be zero on . With + = .
The electrical boundary conditions are:
= 0 = ,
(31)
= = 0 ,
(32)
= 0 = ,
(31)
= = 0 ,
(32)
where , , , , are respectively components of displacement, force, electrical potential, surface charge and
If the behavioral equ
atio
ns (23)
and
-
are substituted in
normal unit outward vector components.
-
Variational piezoelectric equations:
-
For arbitrary space-variables and the admissible virtual displacement and potential , the equation. (21) and
(22) are equivalent to:
(, + ) + (, ) = 0
(33)
The partial integration of this equation and the use of the divergence theorem leads to:
, + +
, + = 0
(34)
Using the stress tensor symmetry property, the natural boundary conditions (29), (31) and the electric potential field relationship (28) result:
+ +
+ = 0
(35)
Substituting the constitutive equations (23) and (24) in equation (35) carries to the electric potential (or field) based in variational principle, which is the starting point of finite
the first and fourth integrals of the equation (39), the Hamilton
principle (37) is reduced to the functional Lagrangian =
= = 0
(40)
= = 0
(40)
, for an arbitrary and possible,
Introducing the following electromechanical energy , and external mechanical and electrical body work and surface forces and charges,
1
= 2 ( )
(41)
= +
(42)
The following relation between the extended potential energy, the electromechanical energy and work of external mechanical and electric loads is obtained,
= +
(43)
This leads to the more common form of the variational equation (40) when relations (23) and (24) are used in Eq. (41); i.e., for admissible and ,
+ = 0
(44)
/tr>
A widely used to assumption is the through-thickness linear variation of the electric potential. Consequently, the induced potential is systematically neglected. To illustrate the influence of this hypothesis on the piezoelectric coupling, the electric potential is decomposed into a linear part 0, known from the prescribed potentials, and an unknown part , representing the induced potential [7],
= 0 +
(45)
= 0 +
(45)
element formulations using independent variables and . In
this case, this could be seen as the sum of the conventional
principle of virtual displacements (first line of (35)), and that of the virtual electric potential (second line of (35)), as suggested in [4-6].
Let us suppose that and depend on time, the following expression holds:
Using this decomposition, together with the constitutive equations (23) and (24), and after some manipulations, Eq.
[ ] () + () + =
+ + + +
+ [(0) + (0)]
(46)
[ ] () + () + =
+ + + +
+ [(0) + (0)]
(46)
(35) becomes,
and when
0
1 1
1
= ( )
2
(36)
1 1
1
= ( )
2
(36)
it is
0
used in equatio
n (35), we obtain the extended
Hamilton principle, for an arbitrary space
in time and disappearing at 0 and 1:
and variables
1
( ) = 0
0
(37)
1
( ) = 0
0
(37)
() is negligible, since for actuation problems only, the variations of the electric field components are zero and electric charges are often not considered, the above equation reduces to
+ =
+ + (0)
(47)
Starting from the virtual work principle, the basic equation of the finite element method is given by
( ) = ( )
(48)
1
{2} = []{}
3
(49)
1
{2} = []{}
3
(49)
In matrix notation, the components of the displacement and the components of the electric potential can be written as:
with
{} = []{}
(50)
Where {} and {} represent the nodal values of displacement and electric potential respectively. Matrices [] and [] include the shape functions.
{11 22 33 23 13 12} = {} = []{}
(51)
{1 2 3} = {} = []{}
(52)
Hence the interpolation matrix [] and [] comes from the differentiation [] and []. In a discretized form, equation
-
becomes:
([ ]{} [ ]{})
= ([]{} [] )
(53)
where {} contains the stress components, {} electrical displacement components, and {} compression components. Using the equations constituting ferroelectricity (25), (26) by deleting the piezoelectric response terms the linear form of finite element resolution:
[] {} = {} + {} (54)
[] {} = {} + {} (55)
where [] is the stiffness matrix for nodal displacement {}
9
(, ) = (, )
=1
9
(, ) = (, )
=1
9
(, ) = (, )
{ =1
(63)
9
(, ) = (, )
=1
9
(, ) = (, )
=1
9
(, ) = (, )
{ =1
(63)
is the nodal load due to the tractions T on and {} is the
[] is the stiffness matrix for nodal potentials, {} is the vector from the surface charge density, {} is another vector due to the presence of spontaneous polarization.[] = [][][] (58)
{} = [] {}
(59)
{} = [ ][][]
(60)
III- HIGH PRECISION ELEMENT OF TYPE LAGRANGE.
The choice of meshes in order to obtain the precise results is a primordial question in modeling. In this article, we use the complete quadratic element (quadrilateral, 9 nodes), this element is often used in fluid mechanics, and now we use to solve the electromechanical problem: active materials (piezoelectric / ferroelectric).
Approach used to generate quadratic mesh is a posteriori approach. The first step is to generate a linear mesh. The second step is to create the middle node by interpolation in the space of parametric coordinate. Interpolation do not guarantee that the node interpolated are in the middle of the ridge of quadrilateral. So, we add another step to check if the node is in the middle of the ridge after interpolation and adjust it if it is not in the middle.
To simplify analytical definition of complex elements, we introduce reference elements and let be (, ) reference space of parametric coordinate.
Relation between the two coordinate espace is given by transformation below :
= (, )
{ = (, )
= (, )
(61)
(, ) = ((, ), (, ), (, ))
(62)
(, ) = ((, ), (, ), (, ))
(62)
With And
additional nodal load due to the spontaneous str
ain.
The stiffness matrix [], column vectors {} and {} are given as :
[] = [][][] (56)
{} = []{}
(57)
{} = [ ][]{}
(58)
Where are quadratic interpolation fonction of reference quadrilateral.
Fig.2. Complete quadratic element (quadrilateral, nine-node) [1]
() = ((), (), ())
(64)
() = ((), (), ())
(64)
Generally, rdge equation of quadrilateral can be definide by : With
TABLE II. Characteristic of PZT-5H:
1 1
() = ( 2) + ( + 2) +
(1 2)
Piezo ceramique PZT-5H
Elastic Coefficient [Pa]
11 = 22 = 1.27205 1011
33 = 1017436 1011
44 = 55 = 2.29885 1010
66 = 2.34742 1010
12 = 21 = 8.02122 1010
13 = 31 = 23 = 32
= 8.46702 1010.
Permittivity Dielectric
11 = 22 = 1704.4
33 = 1433.6
Piezoelectric constants [2]
13 = 6.62281
33 = 23.2403
15 = 17.0345
Piezo ceramique PZT-5H
Elastic Coefficient [Pa]
11 = 22 = 1.27205 1011
33 = 1017436 1011
44 = 55 = 2.29885 1010
66 = 2.34742 1010
12 = 21 = 8.02122 1010
13 = 31 = 23 = 32
= 8.46702 1010.
Permittivity Dielectric
11 = 22 = 1704.4
33 = 1433.6
Piezoelectric constants [2]
13 = 6.62281
33 = 23.2403
15 = 17.0345
2 2
1 1
() = ( 2) + ( + 2) +
(1 2) (65)
2
1
2
1
() =
( 2) +
2
( + 2) + (1 2)
2
Where (xn, yn, zn), (xm, ym, zm) et (xt, yt, zt) are respectively coordinate of first node, middle node and the second node of ridge sont respectivement Length of arc, denote by S, of
vectorial curve p(u) are :
= ,
(66)
() () ()
= ( , , ),
(67)
() 2 () 2 () 2
= ( ) + ( ) + ( )
(68)
Where
B. Results of the simulations
The voltage applied to the terminals of a material PZT-5H is an amplitude A = 100 [V]
Electrical displacements are given as follows.
() 1 1
= 2 (1 2) + 2 (1 + 2) (2)
() 1
=
1
(1 2) +
(1 + 2)
(2) (70)
()
In
2
1
2
1
= 2 (1 2) + 2 (1 + 2) (2
the follow, let be :
() =
)
(71)
+1 +
() = ( + ) .
1 2 2 2
(72)
= 2 (1 2) + 2 (1 + 2) (2
the follow, let be :
() =
)
(71)
+1 +
() = ( + ) .
1 2 2 2
(72)
Using Gauss quadrature in [-1 1], we have :
+
() = ( + ).
2 2 2
=1
(73)
Mesh quadratic Q9 is shown in fig.3.
Fig.3. Mesh quadratic Q9-2D
-
-
MODELING RESULTS
A. Piezoelectric material
The material we used on modeling and optimization in structure is the PZT-5H, often used as an actuator. For simplicity, a 2D finite element model is proposed in this paper by applying structure optimization.
Fig.4. Electric displacement PZT-5H ELEM Q4
Fig.5. Electric displacement PZT-5H ELEM Q9.
Q4: linear element (quadrilateral 4-nodes), Q9: quadratic element (quadrilateral 9-nodes).
We see that using quadratic mesh Q9 give better repartition of electrical displacement D in fonction of electrical field, in the part of solicitation. So with the meshes Q9, the result electric displacement is more precise compared to the result obtained by the mesh Q4.
-
Ferroelectric material
The piezoelectric ceramic is prepared from ferroelectric ceramic grains. Polycrystalline ceramic materials consist of a large number of grains randomly oriented in three-dimensional space. A local coordinate system (, , ) parallel to the network axis can represent the orientation of each grain. The material property and the spontaneous state are built using by the local coordinate system. A global coordinate system is necessary for the response to be measured. We use the Euler angles to connect the three reference frames and the rotating tensor components of the unit between the three coordinate systems:
(74)
The components of the permittivity are transformed into:
3
(, , ) =
,=1
(75)
3
(, , ) =
(76)
3
(, , ) =
(76)
The components of the piezoelectric tensor are transformed into:
,,=1
The tensor components of the elastic modulus are transformed
The voltage applied to the terminals of a material ferroelectric ceramic 8/65/35 PLZT is triangular shape of amplitude A = 7000 [V] give an electric field
= ±0.7 [. 1]
Fig.6. Electric field applied in ferroelectric.
The nonlinear constitutive behavior of the ferroelectric body subjected to the axial electromechanical loading is studied in this section. The comparison of simulated and measured results is shown in Fig.7 and Fig.8. The solid lines correspond to the experimental results of Lynch [8].
into:
3
(, , ) =
,,,=1
(77)
-
Boundary conditions
In Lynchs experimental work [8], a sample of 10 mm cube was used to obtain the experimental data. Hence strictly speaking, 3D finite element simulation should be adopted. But the 3D finite element simulation costs too much computation time and the efficiency of the 3D finite element simulation is quite low due to the complex distribution of the polarization direction. For simplicity, a 2D finite element model is proposed in this paper. Suppose that the polarizations of all the domains in the specimen are parallel to the plane. The polarizations of all the domains will not change along the y direction.
In the calculation, we use ferroelectric ceramic material titanate, lanthanum doped lead zirconate, which Lynch was used in the experimental study, includes a substitution of 8% lanthanum atom in a 65% lead zirconate, 35% lead titanate, solid solution, rated 8/65/35 PLZT, with a ypical grain size of 5m. This composition is very close to a morphotropic rhombo-hetragonal limit. The parameters are in the table III.
|
ferroelectric ceramic 8/65/35 PLZT |
|
|
Young's modulus [Pa] |
E = 68 109, |
|
Coefficient of Poisson |
= 0.3 |
|
Ferroelectric coefficient [] |
33 = 0.9 109, 13 = 4.5 1010, 15 = 1.545 109. |
|
Coercive electric field [MV/m] |
= 0.36 |
|
Spontaneous polarization[C/2] |
= 0.302 |
|
Spontaneous deformity |
= 0.203% |
|
ferroelectric ceramic 8/65/35 PLZT |
|
|
Young's modulus [Pa] |
E = 68 109, |
|
Coefficient of Poisson |
= 0.3 |
|
Ferroelectric coefficient [] |
33 = 0.9 109, 13 = 4.5 1010, 15 = 1.545 109. |
|
Coercive electric field [MV/m] |
= 0.36 |
|
Spontaneous polarization[C/2] |
= 0.302 |
|
Spontaneous deformity |
= 0.203% |
TABLE III. Parameters ferroelectric ceramic 8/65/35 PLZT
Fig.7. Hysteresis loops D-E step time dt=1.
The first result is the loop hysteresis obtain by application of electric field = ±0.7 [. 1]. The step time used is 1 second. Initial loading from the reference state with a controlled voltage ramp starts at the center of the hysteresis loop. As the electric field is increased, the electric displacement increases non-linearly. When the electric field reaches 0.4 MV/m, the sample switches to the polarized state.
Fig.8. Hysteresis loops D-E step time dt=0.1.
When the electric field is zero, the sample has remanent polarization of +0.25[. 2].
A negative electric field (opposite the direction of polarization) is applied. We can see the electric field reached the coercive field, -0.36 MV/m. Here is the direction of polarization of the sample begins to switch.
We see also the polarization reversed and is aligned with the negative electric field.
IV- CONCLUSIONS
Here, finite element formulation is developped. The code was implemented in Matlab. We see that using quadratic meshes give better precision that using linear meshes. Plot of hysteresis loop show that result obtain by model we developped agree with literature result.
REFERENCES
-
M. Kamlah, A.C. Liskowsky, R. M. McMeeking, and H. Balke, "Finite Element Simulation of a Polycrystalline Ferroelectric Based on a Multidomain Single Crystal Switching Model", International Journal of Solids and Structures, vol. 42, pp. 29492964, 2005.
-
S.C. Hwang, C.S. Lynch, and R.M. McMeeking,"Ferroelectric/Ferroelastic Interactions and a Polarization Switching Model", Acta metall. mater., vol. 43, pp. 2073-2084, 1998.
-
J.E. Huber, N.A. Fleck, and R.M. McMeeking, "A Crystal Plasticity Model for Ferroelectrics", in 29th Conference of the UK Dielectrics Society incorporating the 4th UK Transducer Materials and Transducers Workshop, 6-9 April 1998, Switzerland, pp. 39-52, 1999.
-
P. Gaudenzi, K.J. Bathe. An iterative finite element procedure for the analysis of piezoelectric continua. International Journal of System and Structures, 1995.
-
K. Ghandi, N. W. Hagood. Nonlinear finite element modeling of phase transitions in electro-mechanically coupled material, Smart Structure Material. 1996.
-
K. Ghandi, N. W. Hagood. A hybrid finite element model for phase transition in nonlinear electro-mechanically coupled material., Smart Structure Material. 1997.
-
M. Rahmoune, D. Osmont, A. enjeddou, R. Ohayon. Finite element modeling of a smart structure plate system. Seventh International Conference Adapt. Struct. , 1997.
-
C.S. Lynch, S.C. Hwang, R. M. McMeeking. Ferroelectric/ferroelastic interactions and a polarization switching model. Acta metall. mater., 1995