On a non-linear finite element method for piezoelectric structures made of hysteretic ferroelectric ceramics

Computational Materials Science 19 (2000) 81±86

www.elsevier.com/locate/commatsci

On a non-linear ®nite element method for piezoelectric structures made of hysteretic ferroelectric ceramics Marc Kamlah *, Ulrich B ohle, Dietrich Munz Forschungszentrum Karlsruhe, Institut f ur Materialforschung II, Technik und Umwelt, Postfach 3640, D-76021 Karlsruhe, Germany

Abstract A simple macroscopic constitutive law for ferroelectric and ferroelastic hysteresis e€ects of piezoceramics is presented. This constitutive model has been implemented in the public domain ®nite element code PSU of Stuttgart University (Germany). The fully coupled electro-mechanical boundary value problem for our hysteresis model is solved by incrementation of the loading history. In order to verify the capabilities of our tool, a multilayer-like actuator geometry is analyzed. Concerning the residual stresses present after poling, a tensile stress ®eld perpendicular to the direction of the electrodes can be found in the passive region of the actuator where so-called poling cracks are known to occur. It is the key feature of our ®nite element tool that it allows the consideration, in structural mechanical analyses, of such phenomena as they are induced by the remanent hysteresis properties of piezoceramics. Ó 2000 Elsevier Science B.V. All rights reserved. Keywords: Piezoelectricity; Ferroelectricity; Finite element method; PZT-ceramics; Piezoceramic

1. Introduction In order to assess the reliability of a component, one has to know as precisely as possible the mechanical stress state which, in turn, depends strongly on the constitutive properties of the material in question. Besides their reversible piezoelectric behavior, PZT-ceramics exhibit, due to domain switching processes, signi®cant irreversible ferroelectric and ferroelastic hysteresis properties if loaded by electric or mechanical loads of sucient magnitude. For recent experimental investigations of the macroscopic phenomena induced by microscopic domain switching see [1±3].

*

Corresponding author. Fax: +49-7247-822-347. E-mail address: [email protected] (M. Kamlah).

As long as the loadings are small enough, the response behavior of piezoceramics can be approximated by the well-established classical linear piezoelectric theory. However, the restriction to linear piezoelectric behavior is no longer justi®ed in general for present day applications involving complicated geometries and severe loadings. Thus, for a representative stress analysis, the equilibrium condition and the Gaussian law have to be solved in conjunction with an appropriate constitutive law relating, on a macroscopic level, stress and polarization to the histories of strain and electric ®eld. For microscopically motivated approaches to describing the constitutive behavior of piezoceramics see [4±11]. A purely phenomenological approach has ®rst been applied by Chen et al. [12] followed by the work of Bassiouny and Maugin

0927-0256/00/$ - see front matter Ó 2000 Elsevier Science B.V. All rights reserved. PII: S 0 9 2 7 - 0 2 5 6 ( 0 0 ) 0 0 1 4 2 - 7

82

M. Kamlah et al. / Computational Materials Science 19 (2000) 81±86

[13,14]. A simple phenomenological model for the non-linear, electro-mechanically coupled constitutive behavior of ferroelectric and ferroelastic piezoceramics was proposed in [15±17]. In the this paper, we summarize the simple constitutive law just mentioned in Section 2. Then, we present the ®nite element implementation of the constitutive law and discuss the capabilities of the ®nite element tool by considering a multilayer-like geometry in Section 3. In the following analysis, which is restricted to a geometrically linear setting, all component representations of tensors are referred to a Cartesian coordinate system (summation convention). Firstorder tensors (vectors) are denoted by upright ~ and secondletters with superscript arrows (~ a, A) order tensors (tensors) by bold italics …a; A; a†. A dot between tensors indicates contraction relative to one index, for example, the inner product between vectors, i.e., ~ a
~ b ˆ ai bi , the composition of two tensors, i.e., A
Aÿ1 ˆ I (Aÿ1 : inverse of A; I: identity tensor), the inner product between tensors, i.e., A
B ˆ tr …A
BT † ˆ Aij Bij (tr A: trace of A; B T : transpose of B) and the linear mapping of a vector by a tensor, i.e., A
~ a ˆ Aij aj . AD ˆ A ÿ …1=3†…tr A†I p  stands for the deviator of A, and kAk ˆ A
A is the norm of the tensor A. ~ a=k~ ak is thepunit ea ˆ ~ vector in the direction of the a
~ a: norm of ~ a). The symbol

vector ~ a (k~ ak ˆ ~ indicates the dyadic product, e.g., ~ a ~ b yields a second-order tensor with components ai bj . … _ † ˆ d… †=dt denotes time derivative of a ®eld ( ). Further mathematical de®nitions will be given where they are needed.

The basic assumption of our model is to introduce the irreversible or remanent polarization ~ Pi i and the irreversible or remanent strain S as internal variables by the additive decompositions ~ Pi ; Pˆ~ Pr ‡ ~

…1†

S ˆ Sr ‡ Si

…2†

of the polarization and the strain, respectively [13]. The macroscopic internal variables ~ Pi and S i must be understood as averages of the corresponding microscopic quantities, i.e., the spontaneous polarization and the spontaneous strain of a domain. Please note that the notion of irreversibility is introduced here as a synonym of remanence and without explicit reference to a thermodynamic framework. The parts ~ Pr and S r are assumed to be reversible in the sense that they vanish if the loads ~ E and T are zero. As a simple choice, we adopt for them equations of linear piezoelectric structure: ~ E; Pr ˆ d : T ‡ 
~

…3†

E: S r ˆ Cÿ1 : T ‡ dT
~

…4†

2. A simple constitutive model for ferroelectric and ferroelastic piezoceramics

, d and C are the (second-order) tensor of dielectric permittivities, the (third-order) tensor of piezoelectric moduli and the (fourth-order) tensor of elastic moduli, respectively. In order to keep the model as simple as possible, we neglect the in¯uence of the loading history on the anisotropy of  and C and assume isotropic representations, where , l, and m are the dielectric constant, shear modulus and Poisson ratio, respectively (I; I: identity tensors of second and fourth order, respectively):   m I I : …5†  ˆ I; C ˆ 2l I ‡ 1 ÿ 2m

For a macroscopic theory of a deformable dielectric, a relation for the dependence of the electric ®eld ~ E ˆ ÿgrad u and the stress tensor T on the histories of the polarization ~ P and the strain T tensor S ˆ …1=2†…grad~ u ‡ …grad~ u† † is needed (u; ~ u: electric potential and displacement vector, respectively). In this section, we summarize our proposal for such a constitutive model, as it has been motivated in full detail in [15±17].

Note that this restriction is not essential for our approach (for similar assumptions see [4] and related papers). The situation is completely di€erent for the piezoelectricity tensor. The history dependence of the anisotropy cannot be neglected here, since the piezoelectric properties are not just modi®ed as is the case with dielectricity and elasticity. Rather, the phenomenon of macroscopic piezoelectricity is

M. Kamlah et al. / Computational Materials Science 19 (2000) 81±86

absent in the unpoled state, i.e., if there is no macroscopic remanent polarization ~ Pi . Furthermore, in the fully poled state, we have transversely isotropic piezoelectricity where the axis of anisotropy coincides with the direction ~ ePi of poling. Therefore, we assume the representation  k~ Pi k Pi † ˆ dk ei ej ek ‡ d? …dij ÿ ei ej †ek dkij …~ Psat   1 ‡ dˆ ……dki ÿ ek ei †ej ‡ …dkj ÿ ek ej †ei † 2 …6† for the piezoelectricity tensor, where the constants dk ; d? ; dˆ correspond to the well-known parameters d33 ; d31 ; d15 , respectively, and Psat is the maximum remanent polarization (dkij ˆ …d†kij ; ei ˆ ePi ). …~ ePi †i is the ith component of ~ The macroscopic remanent polarization ~ Pi stands for the net polarization of a neighborhood of many polarized domains with possibly di€erent orientations and is a measure of the extent of alignment of the domains in a certain direction. The change of alignment of the domains causes a change of the remanent strain of this neighborhood and it may be represented as uniaxial strain state in the direction of ~ Pi . Thus, we relate the i remanent polarization ~ P to the part   i 3 k~ Pk 1 p ~ ePi ~ …7† e Pi ÿ I S ˆ Ssat 2 Psat 3 of the remanent strain (Ssat : maximum remanent polarization). Since experimental ®ndings indicate that changes of the macroscopic strain state due to domain switching are volume preserving, we constructed S p such that it yields a deviatoric uniaxial strain state in the direction of ~ Pi [1]. Domain switching can also be induced by mechanical stresses. In the case of an initially unpoled material, the state of macroscopic polarization remains unchanged, while purely ferroelastic irreversible strain changes occur. Therefore, we also have to consider a ferroelastic contribution S f to the remanent strain, the latter consisting in the general case of a superposition of both parts: i

p

f

S ˆS ‡S :

…8†

83

Since, according to (7), S p is simply a function of ~ Pi , the constitutive law has to be completed by evolution laws for the remanent polarization ~ Pi f and for the ferroelastic part S of the remanent strain, the latter allowing for deviatoric ferroelastic strain changes only. The evolution laws are de®ned on the basis of four di€erent loading criteria: Pi k ÿ Ec ; E; ~ Pi † ˆ k~ E ÿ cp~ f p …~

…9†

^ sat …T; ~ E; ~ Pi † ˆ k~ Pi k ÿ P E; ~ Pi †; …10† hp …T; ~ r 3 D k…T ÿ cf S f † k ÿ r^c …~ E; ~ Pi ; S f † ˆ E; ~ Pi †; f f …T; ~ 2 …11† ! r r 2 f 2 p kS k ÿ Ssat ÿ kS k : Pi ; S f † ˆ …12† hf …~ 3 3 The ®rst and the third functions indicate the onset of domain switching due to electric and mechanical loading, respectively. The second and fourth criteria characterize the fully switched ferroelectric and ferroelastic domain states, respectively. Here, the ®eld-dependent coercive stress is given by * + ~ E i E; ~ P † ˆ rc ‡ n
~ ei ; …13† r^c …~ Ec P while the history-dependent saturation value of the remanent polarization is de®ned by ^ sat …T; ~ E; ~ Pi † ˆ P    1 3 D e i
T
~ ePi ÿ r^c ÿ ~ ‡ Pd : …Psat ÿ Pd † 1 ÿ m 2 P …14† As rough approximation for a soft PZT, the material constants can be chosen according to Table 1. Introducing the notations   1; x P 0 1; x > 0 bxe ˆ ; dxc ˆ …15† 0; x < 0 0; x 6 0 and

d p f ˆ f ; dt ~P_i ˆ0  p

d p h ˆ h ; dt ~P_ iˆb f p ebf p ekp of p =o~E  p

f

…16†

84

M. Kamlah et al. / Computational Materials Science 19 (2000) 81±86

Table 1 Material constants chosen as rough approximation for a typical soft PZT 1:5  10ÿ8 3:5  10ÿ10 1:0  106 5:0  107 5:0  107

 dk Ec rc n



ff ˆ

C/(V m) m/V V/m N=m2 N=m2

d f f ; dt S_ f ˆ0

as well as i

i

 i

F ˆ bf ebf e;



hf ˆ

Y d? Psat Ssat m

8:0  1010 )1:4  10ÿ11 0.3 0.002 1:5  108

d f h ; …17† dt S_ f ˆbf f ebf f ekf of f =oT f

   i i i H ˆ bh edh c ‡ df c ; i

i ˆ p; f;

…18†

and ®nally p

oh ~ bˆ i; o~ P

Nf ˆ

ohf =oS i ; kof f =oS i k

…19†

the evolution laws for the internal variables can be written as   p _i ohp p p p of ~ eb ~ ‡ H p kph i eb †
F k f P ˆ …I ÿ H ~ o~ E o~ P …20† and S_ f ˆ …I ÿ H f N f N f † :

  of f ohf F f kff ‡ H f kfh f : oT oS …21†

The factors of proportionality kpf , kph , kff and kfh are determined from the consistency conditions f_ p ˆ 0, h_ p ˆ 0, f_ f ˆ 0 and h_ f ˆ 0, respectively. For a detailed motivation and discussion of these relations see [15±17]. 3. Finite element implementation The constitutive law for ferroelectric and ferroelastic hysteresis e€ects sketched in the previous section has been implemented in the public domain non-linear multipurpose ®nite element code PSU provided by Stuttgart University (for more infor-

N=m2 m/V C=m2 ) N=m2

m dˆ cp cf Pd

0.35 5:2  10ÿ10 2:0  106 5:0  109 0.1

) m/V V m/C N=m2 C=m2

mation see the web page http://www.isd.uni-stuttgart.de/arbeitsgruppen/psu-www/index.html). The constitutive model of Section 2 can be written as a system of ordinary di€erential equations with respect to time. The factors of proportionality kpf , kph , kff and kfh were calculated analytically from the corresponding consistency conditions for the case of prescribed electric ®eld ~ E and strain S. For the numerical integration of the resulting system of di€erential equations, we chose a highly accurate, explicit Runge±Kutta scheme. The points of change from one loading condition to another are identi®ed by cutting back subincrements iteratively until the time instant of the change is identi®ed with sucient accuracy. There are two options available for passing the material-dependent contribution to the tangent sti€ness operator to the non-linear ®nite element solver. The ®rst one is obtained by numerical differentiation of the constitutive law by symmetric di€erence quotients. As the second approximate choice, by means of the tensors , C and d according to Eqs. (5) and (6), respectively, a tangent sti€ness operator of the structure of linear piezoelectricity is constructed. For a detailed description see [19]. We now want to apply our ®nite element tool to the example structure shown in Fig. 1. Exploiting symmetries, the model represents a single layer of a multilayer geometry. The dimensions (width 410 lm, height 57:5 lm, lower electrode terminates at the middle of the model) were chosen as adapted from the literature [18,20]. Fig. 1 also shows the boundary conditions applied and the ®nite element mesh. Due to symmetry conditions, the upper electrode must not deform or rotate and is only allowed to translate parallel to its starting position. Concerning the strain state, the plane strain assumption has been employed. In order to simulate a poling process, the electric potential

M. Kamlah et al. / Computational Materials Science 19 (2000) 81±86

85

Fig. 1. Boundary conditions and ®nite element mesh for the electro-mechanical analysis of a single actuator layer.

Vappl …t† is prescribed as triangular loading history in the following manner: it ®rst is increased from zero to )170 V until t1 ˆ 100 s and then reduced to zero until t2 ˆ 200 s. We consider the residual stress state after unloading induced by the poling process. In Fig. 2, we see the normal stresses perpendicular to the electrodes plotted along the lower border of our model. As can be expected, stresses are stronger at the instant t ˆ t1 of maximum voltage than they are after unloading (t ˆ t2 ). However, even after removing the external electric load, we observe a signi®cant residual stress state. As a potential danger for the actuator, a tensile stress state of over 20 MPa is found in the passive region towards the free right border. This is exactly the region where the so-called poling cracks are known to occur. In this context please note that the tensile strength of piezoceramics may be as low

as the order of magnitude of these stresses induced by poling. The tensile stresses in the poled region are caused by the vertical extension of the area between the electrodes during poling. 4. Conclusion A simple phenomenological model for the ferroelectric and ferroelastic properties of piezoceramics has been presented. The model has been implemented in the non-linear ®nite element code PSU. As demonstration example, the ®nite element tool is applied to a multilayer-like structure. It turns out that the response of the structure is dominated by the presence of ferroelectric and ferroelastic large signal hysteresis e€ects. A signi®cant tensile stress state remains after poling in the passive region between the electrode tip and the free border of the actuator. Even if, after poling, the actuator is operated in the small signal regime with small stress changes, the fact that they are superposed on the residual stress state will increase the danger of accelerated fatigue. Acknowledgements Support by the Deutsche Forschungsgemeinschaft is gratefully acknowledged. References

Fig. 2. Normal stresses perpendicular to the electrodes plotted along the lower border of the layer (the horizontal line represents the lower electrode). Dotted line: t ˆ t1 , solid line: t ˆ t2 .

[1] H. Cao, A.G. Evans, Nonlinear deformation of ferroelectric ceramics, J. Am. Ceram. Soc. 76 (1993) 890±896.

86

M. Kamlah et al. / Computational Materials Science 19 (2000) 81±86

[2] A. Schaufele, K.H. H ardtl, Ferroelastic properties of lead zirconate titanate ceramics, J. Am. Ceram. Soc. 79 (1996) 2637±2640. [3] C.S. Lynch, The e€ect of uniaxial stress on the electromechanical response of 8/65/35 PLZT, Acta Mater. 44 (1996) 4137±4148. [4] S.C. Hwang, C.S. Lynch, R.M. McMeeking, Ferroelectric/ ferroelastic interactions and a polarization switching model, Acta Metall. Mater. 43 (1995) 2073±2084. [5] Y. Huo, Q. Jiang, Modeling of domain switching in polycrystalline ferroelectric ceramics, Smart Mater. Struct. 6 (1997) 441±447. [6] W. Chen, C.S. Lynch, A micro-electro-mechanical model for polarization switching of ferroelectric materials, Acta Mater. 46 (1998) 5303±5311. [7] Th. Michelitsch, W. Kreher, A simple model for the nonlinear material behavior of ferroelectrics, Acta Mater. 46 (1996) 5085±5094. [8] Th. Steinkop€, Finite-element modelling of ferroic domain switching in piezoelectric ceramics, J. Eur. Ceram. Soc. 19 (1999) 1247±1249. [9] S.C. Hwang, R.M. McMeeking, A ®nite element model of ferroelectric polycrystals, Ferroelectrics 211 (1998) 177± 194. [10] S.C. Hwang, R.M. McMeeking, A ®nite element model of ferroelastic polycrystals, Int. J. Solids Structures 36 (1999) 1541±1556. [11] J.E. Huber, N.A. Fleck, C.M. Landis, R.M. McMeeking, Constitutive model of ferroelectrics, J. Mech. Phys. Solids 47 (1999) 1663±1697. [12] P.J. Chen, Three dimensional dynamic electromechanical constitutive relations for ferroelectric materials, Int. J. Solids Structures 16 (1980) 1059±1067.

[13] E. Bassiouny, A.F. Ghaleb, G.A. Maugin, Thermomechanical formulation for coupled electromechanical hysteresis e€ects ± I. Basic equations, II. Poling of ceramics, Int. J. Eng. Sci. 26 (1988) 1279±1306. [14] E. Bassiouny, G.A. Maugin, Thermomechanical formulation for coupled electromechanical hysteresis e€ects ± III. Parameter identi®cation, IV. Combined electromechanical loading, Int. J. Eng. Sci. 27 (1989) 975±1000. [15] M. Kamlah, C. Tsakmakis, Phenomenological modeling of the non-linear electro-mechanical coupling in ferroelectrics, Int. J. Solids Structures 36 (1999) 669±695. [16] M. Kamlah, U. B ohle, D. Munz, Ch. Tsakmakis, Macroscopic description of the non-linear electro-mechanical coupling in ferroelectrics, in: V.V. Varadan, J. Chandra (Eds.), Proceeedings ± Smart Structures and Materials 1997: Mathematics and Control in Smart Structures, vol. 3039, SPIE, Bellingham, 1997, pp. 144±155. [17] M. Kamlah, U. B ohle, Finite element analysis of piezoceramic components taking into account ferroelectric hysteresis behavior, Int. J. Solids Structures 38 (2000) 605± 633. [18] X. Gong, Z. Suo, Reliability of ceramic multilayer actuators: a nonlinear ®nite element simulation, J. Mech. Phys. Solids 44 (1996) 751±769. [19] U. B ohle, Ph anomenologische Modellierung und FiniteElemente-Simulationen von nichtlinearen elektromechanischen Vorg angen in ferroelektrischen Materialien, dissertation thesis, University of Karlsruhe, 1999. [20] C.L. Hom, N. Shankar, A ®nite element method for electrostrictive ceramic devices, Int. J. Solids Structures 33 (1996) 1757±1779.