- Email: [email protected]
Nonlinear finite element modeling of polycrystalline ferroelectrics based on constrained domain switching
Computational Materials Science 44 (2008) 322–329
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Nonlinear finite element modeling of polycrystalline ferroelectrics based on constrained domain switching F.X. Li, R.K.N.D. Rajapakse * Department of Mechanical Engineering, The University of British Columbia, 2054-6250 Applied Science Lane, Vancouver, Canada V6T 1Z4
a r t i c l e
i n f o
Article history: Received 9 December 2007 Received in revised form 16 March 2008 Accepted 25 March 2008 Available online 21 May 2008 PACS: 46.15.Cc 77.80.Fm 77.84.s Keywords: Constitutive models Finite element method Ferroelectrics Domain switching Charge screening
a b s t r a c t A review of the nonlinear finite element formulations for ferroelectrics shows that instability related to domain switching may arise in formulations using a phenomenological constitutive model or micromechanical domain switching criterion. To overcome the instability problems and to predict the experimentally observed incomplete non-180° switching in ferroelectrics, a constrained domain switching model for polycrystalline ferroelectrics is proposed in this paper and implemented in a finite element code. In the proposed model, a ferroelectric polycrystalline is made up of numerous random oriented grains, each of which contains multi-domains and is represented by a finite element. Charge screening effect in real ferroelectrics is also taken into account in the model, thus 180° domain switching is not constrained. While the internal stress induced by non-180° switching cannot be compensated in a similar manner and is considered naturally by the finite element method. The fraction of non-180° switching in each grain is obtained through a dichotomy approximation. Selected numerical results show that during electrical poling, 180° switching is usually complete while 90° switching is only about 2% for both BaTiO3 and tetragonal PZT ceramics, which agree well with the experimental results. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Ferroelectric materials such as Barium Titanate and PZT ceramics are widely used in smart structures and devices as actuators, sensors, ultrasonic motors, etc. [1]. Domain switching is now accepted as the predominant source of the nonlinear response of these ceramics [2]. Ferroelectric ceramics are also very brittle and are subjected to complex electromechanical loading in practical applications that may cause domain switching and cracking. It is therefore important to develop nonlinear finite element models to analyze ferroelectric materials under complex loading. Since the original work of Allik and Hughes [3], linear finite element analysis of piezoelectric materials has received considerable attention. Hom and Shankar [4], and Gong and Suo [5] independently presented nonlinear finite element formulations to study multi-layer actuators made of relaxor ferroelectric ceramics. A fully coupled nonlinear finite element model for phase transition of electromechanically coupled materials was proposed by Ghandi and Hagood [6] by using nodal displacements and electric potential as degrees of freedoms. They later presented a finite element model based on the Hu–Washizu functional, taking the electric displace* Corresponding author. Present Address: Institute for Computing, Information and Cognitive Systems, The University of British Columbia, Canada. E-mail address: [email protected] (R.K.N.D. Rajapakse). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.03.040
ment as nodal variables [7]. Recently, Kamlah and Bohle [8] presented a nonlinear finite element model based on a bilinear phenomenological constitutive model. They first solved a pure dielectric boundary-value problem for the history of the electric potential and thereafter prescribing this electric potential, an electromechanical stress analysis with the mechanical boundary conditions was conducted. Kamlah et al. [9] also conducted finite element simulations by using a multi-domain single crystal switching model. The above mentioned nonlinear finite element (FE) models are based on explicit or implicit phenomenological constitutive models. There are only few experimental studies for ferroelectric materials under multi-axial loading [10–13], and the FE models are therefore often used to simulate one-dimensional problems for the purpose of comparison with experiments [8]. Since the development of the first domain switching criterion for ferroelectric materials by Hwang et al. [14], several nonlinear FE models accounting for domain switching have appeared in the literature [15–18]. Hwang et al. [15,16] formulated two nonlinear FE models for ferroelectric and ferroelastic materials by using an improved domain switching criterion based on potential energy [19]. To simplify the modeling, they solved the mechanical and dielectric boundary-value problems separately, neglecting the piezoelectric effect. Recently, Li and Fang [20] presented an improved domain switching criterion based on the Gibbs free energy and
323
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
developed a fully coupled, three-dimensional (3D) nonlinear FE model to study the effects of intergranular interactions during domain switching. The nonlinear FE models accounting for domain switching usually allow only one domain to switch in each iteration step [15,16,20], thus leading to quite low numerical efficiency. Another problem with the existing FE models is that the energy barrier for 90° (or non-180°) domain switching is not well known, though the energy barrier for 180° switching is widely accepted to be equal to 2P0EC (P0 is the magnitude of the spontaneous polarization and EC the magnitude of coercive field). This paper first examines the so called ‘‘switching-related instability” problem which is likely to be encountered in the FE models accounting for domain switching. A general nonlinear finite element formulation is summarized in Section 2. In Section 3, it is shown that for the nonlinear FE models based on a phenomenological constitutive law, the Newton–Raphson iteration is usually unstable if the electric potential is taken as a nodal variable. In the case of FE models using a micromechanical domain switching criterion, the eigen changes in spontaneous polarization/strain caused by domain switching are found to be equivalent to a rather large depolarization field/stress on the switched domain, which makes the domain switch back to its original state, thus causing switching instability. To overcome these problems and to predict the observed incomplete non-180° domain switching in ferroelectrics [21,22], a constrained domain switching model and its 2D finite element implementation are presented in Section 4. Selected numerical results are presented in Section 5 and the conclusions are in Section 6.
To examine the causes of domain switching related instability, a general nonlinear finite element formulation is summarized in this section [20]. For a ferroelectric material with a remnant polarization and remnant strains, the constitutive equations can be written as
Dj ¼ ejmn ðemn ermn Þ þ
e kjm Em
þ
ð1Þ Prj
ð2Þ
Di;i ¼ 0
ð3Þ
Let ui and / denote the components of the displacement vector and electric potential, respectively. The strains and electric field can be expressed as eij ¼
1 ðui;j þ uj;i Þ; 2
Ei ¼
o/ oxi
uz
/ T
ð9Þ
The corresponding generalized stress, strain and prescribed traction vectors are defined as R ¼ ½rx ; ry ; rz ; rxy ; rxz ; ryz ; Dx ; Dy ; Dz T
ð10Þ
C ¼ ½ex ; ey ; ez ; 2exy ; 2exz ; 2eyz ; Ex ; Ey ; Ez
T
ð11Þ
T ¼ ½tx ; ty ; tz ; pe T
ð12Þ
The constitutive law can be written in the following incremental form. dR ¼ C dC þ DCðC Cr Þ þ dRr
ð13Þ
where C is the generalized constitutive matrix which includes elastic, piezoelectric and dielectric coefficients. Following incremental form of displacement interpolation can be defined for a typical finite element. Du ¼ N Dae
ð14Þ e
where Du; Da and N denote the incremental generalized displacement vector of an element, incremental generalized nodal displacement vector and the interpolation matrix, respectively. The incremental generalized strain vector of an element is given by ð15Þ
where L is the linear differential matrix operator and B = LN. The incremental global finite element equations can be expressed as K Dae ¼ DP
ð16Þ
where X T e K¼ G K G
ð17Þ
e
where rij , emn are the components of the stress tensor and the strain tensor, respectively; Em and Dj are the components of the electric field vector and the electric displacement vector, respectively; e C Eijmn , ejmn and kjm are the elements of the elastic stiffness tensor, piezoelectric tensor and the dielectric permittivity tensor, respectively; and ermn and Prj are the components of the remnant strain tensor and remnant polarization vector, respectively. Consider a piezoelectric body that occupies a volume, V, with a closed surface S. If body forces and body charges can be neglected, the equilibrium equations can be written as, rij;j ¼ 0;
uy
u ¼ ½ ux
DC ¼ LDu ¼ LNDae ¼ BDae
2. Nonlinear finite element formulation
rij ¼ C Eijmn ðemn ermn Þ emij Em
kk denotes the jump across a boundary; pe is the surface charge per unit area; ui denotes the components of the prescribed displacement; and / is the prescribed electric potential. The generalized displacement vector of a finite element is taken as
ð4Þ
The boundary conditions can be expressed as rij nj ¼ ti on St
ð5Þ
kni Di k ¼ pe on SD
ð6Þ
ui ¼ ui on Su
ð7Þ
/ ¼ / on S/
ð8Þ
where nj are the components of the outward unit normal vector of the surface S and t i denotes the components of the traction vector;
and Ke ¼
Z
BT CBdve
ð18Þ
ve
DP ¼
Z
NT DTds þ St þSD
Z
X e
GT
Z
BT CDCr dve þ ve
BT DCðC Cr Þdve ve
X e
GT
Z
X
GT
e
BT DRr dve
ð19Þ
ve
where G is a coordinate transformation matrix that relates generalized element stiffness in local coordinates to the global coordinates and the generalized element nodal vector in local coordinates to the global coordinates. It should be noted that for ferroelectric materials during domain switching, the second term on the right side of Eq. (19) corresponds to the contribution due to the change of the remnant strain and the third term corresponds to the contribution due to the change of the constitutive matrix C, i.e., change of materials properties. The last term of Eq. (19) represents the contribution due to the change of remnant polarization. 3. Domain switching related instabilities in FE modeling of ferroelectrics In this section, the iteration instabilities which might be encountered in nonlinear FE analysis based on phenomenological
324
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
constitutive models and the ‘‘back-switching” phenomenon in nonlinear FE models using a micromechanics based domain switching criterion are discussed. 3.1. Switching induced iteration instability in FE models based on phenomenological constitutive laws For FE models based on phenomenological material models, the material response to applied fields, such as the D–E hysteresis loop, stress–strain curve, can be obtained directly from the material models, explicitly or implicitly, without the help of the FE method. Domain switching (or phase transitions) is not involved in the nonlinear FE formulation, i.e., the second and the last terms on the right hand side of Eq. (19) vanish. The FE formulation of such models is similar to the nonlinear FE formulation of conventional elastoplastic materials except for the additional terms involving dielectric and piezoelectric effects. An incremental nonlinear iterative process is needed in the FE analysis due to the history-dependent behavior of ferroelectric materials. Among the nonlinear iterative methods, the Newton– Raphson method is common and convenient with good accuracy and convergence [23]. Note that in conventional elastoplastic analysis, the Newton–Raphson method is usually stable [24]. However, for ferroelectric materials due to the characteristics of domain switching, the slope of D–E curves near the coercive field is much larger than that of the linear dielectric part. The Newton–Raphson iteration can become unstable due to such abrupt changes of slope. To illustrate this iterative instability, the poling curve of a typical PZT ceramic is shown in Fig. 1. For 1D dielectric problems with the electric potential as the nodal unknowns with prescribed electric displacement D0 , it can be seen from Fig. 1a that Newton– Raphson iteration to solve for the exact solution Ex is unstable, i.e., Ei ! Eiþ1 ! divergence, where Ei and Eiþ1 are the electric field computed at the ith and i + 1th iteration. Landis [25] indicated that by taking the vector electric displacement instead of the scalar electric potential as nodal unknowns, a standard iteration process is stable for ferroelectric materials with a simple bilinear D–E response. However, if the saturation part of the poling curve (BC in Fig. 1b) is also taken into account, even using the electric displacement as nodal unknowns, the Newton–Raphson iterative process for a prescribed electric field may still be unstable. As shown in Fig. 1b, Di is the ith iteration and it is obvious that Di ! Diþ1 ! divergence. This can be true even if the ith iteration, such as Di1 in Fig. 1b, is very close to the exact solution Dx , that ! divergence, and Di ! Diþ1 ! divergence. is Di1 ! Diþ1 1 The above iteration instabilities are caused by the steep slope of the curve AB in Fig. 1a and the curve BC in Fig. 1b associated with domain switching. The D–E curve of a non-polar dielectric material without switching is exactly like the stress–strain curve of a con-
a
ventional elastoplastic material, and the Newton–Raphson iteration is always stable [24]. Therefore, this type of iteration instability is called ‘‘domain switching related iteration instability” to distinguish it from the direct instability of the solution due to domain switching that is discussed in the ensuing section. It should be noted that such ‘‘iteration instabilities” are usually encountered in incremental FE models using phenomenological constitutive laws. It arises from the iterative process when solving for the unknown electric field (electric displacement) with the prescribed electric displacement (electric field). In phenomenological models using numerical methods other than FEM [26–28], the iterations may not be involved in the solution process and above mentioned ‘‘iteration instabilities” may not exist. Moreover, even in FE models using the iterative processes, if the nonlinear items, such as the remnant polarization and strain, are separately taken into account as internal variables [29], ‘‘iteration instabilities” can be avoided by assuming that the linear properties do not change with the internal variables or by using a simplified bilinear D–E assumption [8]. 3.2. Back-switching in FE models utilizing micromechanics based domain switching criteria In the case of nonlinear FE models utilizing micromechanics based domain switching criteria [15–18,20], the equivalent nodal load induced by the change of spontaneous strain and spontaneous polarization during domain switching, i.e., the second and the last terms of the right hand side of Eq. (19), can be rather large compared to the contribution due to moderate electromechanical loading. Consider the case of 180° domain switching without a change in spontaneous strain (Fig. 2a). When the polarization of a domain reverses, a large depolarization field Ed is generated. In the absence of free charge compensation, the depolarization field Ed is on the order of DP=k (10–50 kV/mm), where DP is the magnitude of the change of spontaneous polarization during switching and k is the dielectric constant. Such a large depolarization field (usually higher than 10 times the coercive field) makes the domain switches back to its original state, thus domain switching instability arises. In the case of 90° domain switching (Fig. 2b), in addition to the large depolarization field Ed caused by the spontaneous polarization change, a large mechanical field (tensile stress r1 and compressive stress r2 ) is also generated due to the spontaneous strain change during domain switching. This mechanical field also tends to make the switched domain reverts to its original strain state. r1 and r2 are both of the order of Dc=s (1 GPa, where Dc is the magnitude of the spontaneous strain change and s is the elastic compliance), so that ‘‘back-switching” takes place. Therefore in an unmodified nonlinear FE formulation, domain switching is always unstable due to the terms involved in Eq. (19).
b
D
E
C
C
E
B
D
0
0
B A
A O
i
E
Ex
E
i+1
E
O Di
D
i+1
i D1 Dx D1i+1
D
Fig. 1. Illustration of Newton–Raphson iteration instability: (a) electric potential as nodal unknowns; (b) electric displacement as nodal unknowns.
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
Fig. 2. Depolarization electric field and mechanical field induced by (a) 180° switching or ferroelectric switching and (b) 90° switching or ferroelastic switching.
Domain switching instability associated with the finite element formulation has been noticed in previous modeling studies [16,17,20,25]. To weaken the high electric and mechanical fields induced by ferroelectric and ferroelastic switching in an inclusion model which is equivalent to FE models, Hwang et al. [30] employed a much smaller Young’s modulus (about 20% of the true value) and a rather larger dielectric constant (about 14 times the true value) in their simulations. Chen and Lynch [17] noticed in their modeling, domain switching instability caused by the high local electric field due to the jump of polarization during switching. To reduce this effect, they introduced a factor kð0 < k < 1Þ to the last term of the right hand side of Eq. (19). Kessler et al. [31,32] discussed the unstable switching at the initiation of domain switching in their numerical simulations (not FE model). They suggested that switching instability can be removed by introducing ferroelectric ‘‘hardening” in the material model, which is equivalent to the simplified bilinear D–E curve assumption of Kamlah and Bohle [8]. However, the switching instability addressed in their work is associated with a ‘‘jump” in polarization or strain which is different from the ‘‘back-switching” addressed in this paper. 4. FE model based on constrained domain switching In this section, a constrained domain switching model free from the above mentioned instabilities in the finite element analysis is proposed and implemented in a finite element code. Furthermore, the present model can also predict the experimentally observed incomplete non-180° domain switching in ferroelectric ceramics [21,22,33]. As discussed above, the cause of instability in FE models based on domain switching criteria is the rather large depolarization field/stress induced by the domain switching associated terms. In real ferroelectric ceramics, the large depolarization field by mismatched polarization is fully or partially compensated by free charges inside the material, which is called ‘‘charge screening” [34,35]. The internal stress caused by non-180° switching, however, cannot be compensated in a similar way. The non-180° switching induced constraint in polycrystalline ferroelectrics are actually crystal symmetry dependent. According to the Taylor rule of plasticity [36], a crystal must have at least five slip systems for a
325
polycrystalline to be ductile. In terms of the deformation modes, a ferroelectric crystal is similar to a ferroelastic crystal or a shape memory alloy [37] which has two independent slip systems for the tetragonal and three for the rhombohedral case. Thus in tetragonal or rhombohedral polycrystalline ferroelectric ceramics, non180° domain switching generates large internal stresses, i.e., they are constrained by neighboring grains and may not occur. While in orthorhombic ferroelectric ceramics or PZT ceramics near the MPB, where tetragonal and rhombohedral phases coexist, there are five or more independent deformation modes, which make the materials ductile. Then non-180° domain switching does not generate very large internal stresses from neighboring grains and such switching can almost be accomplished completely in these types of polycrystallines [37,38]. In the present study, we do not study these ‘‘ductile” polycrystallines but focus only on the constrained domain switching phenomena in tetragonal and rhombohedral ferroelectric ceramics. So far, only a few theoretical studies have addressed the experimentally observed incomplete non-180° domain switching [21,22,33] except the work by Li et al. [38] and Wan et al. [39]. Li et al. [38] deduced through a rigorous mathematical analysis that the lower limit of switching strains that can be achieved in tetragonal and rhombohedral ferroelectric ceramics. In the phenomenological model of Wan and Bowman [39], they used an adjustable parameter to address the depolarization field by 90° switching and fit the model to experimental results. As to the FE models, the models using phenomenological material models as well as those based on domain switching criteria with the single-domain grain assumption cannot address the incomplete non-180° switching. The FE simulation by Kamlah et al. [9], which is based on a multi-domain single crystal switching model, could technically predict the fraction of non-180° switching, but they did not specifically address this issue. In the constrained domain switching model proposed in this paper, a ferroelectric polycrystalline is made up of numerous random oriented grains, each of which contains N types of domains (N = 6 for the tetragonal ferroelectrics and N = 8 for the rhombohedral) and is represented by a finite element. In an unpoled ceramic, the fraction of each type of domains in a grain is equal to 1/N, thus both the remnant polarization and strain of each grain is zero. Inside each grain, the electromechanical field is assumed to be uniform. Charge screening effect is taken into account in this model, i.e., the mismatched polarization by domain switching is to be compensated by free charges. Thus, an 180° switching does not generate a depolarization field and not get constrained. The large depolarization stress caused by the spontaneous strain change during non-180° switching cannot be compensated in a similar way and is taken into account by the FE method naturally. As the depolarization stress increases with the fraction of non-180° switching, the achievable fraction of non-180° switching increases with the applied field/stress and is approached by a mathematical approximation. In this model, it is also assumed that a ferroelectric polycrystalline is linear dielectric, elastic and piezoelectric unless domain switching occurs, i.e., all nonlinear polarization and strains are caused by domain switching.
4.1. A simplified 2D finite element model In this paper, the constrained domain switching model is implemented in a 2D finite element code. In the 2D case of tetragonal materials, there only exist four types of domains inside each grain, as shown in Fig. 3a. For a textureless polycrystalline ferroelectrics, the orientations of the grains, which can be denoted by an orientation angle h (0 6 h < 2p) in Fig. 3b, distribute randomly in all directions, i.e., h distributes uniformly in ½0; 2p. Note that in the 2D case,
326
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
a
b
y a
y
a c
1 4
ing cases are considered without loss of any generality. At the bottom face, both the vertical displacement u2 and the electric potential / are fixed for all the nodes. The horizontal displacement of the left node on the bottom face is also fixed to prevent rigid body motion. At the top face, both the electric potential and tractions are prescribed according to the applied electromechanical loading. Both lateral boundaries are traction and charge free.
X2
4
1
θ
x 2
X1 4.4. Domain switching criterion
3
2
3
x
Fig. 3. Domain and grain orientations for tetragonal polycrystalline ferroelectrics. (a) domain orientation in local grain coordinates x–y; (b) grain orientation in sample coordinates X1–X2.
the remnant strain for an unpoled polycrystalline ferroelectric is not zero, but exx ¼ eyy ¼ ðc=a 1Þ=6 ¼ c0 =4
ð21Þ
4.2. Material constants during domain switching Theoretically, all the material constants of a ferroelectric polycrystalline change either substantially or slightly with domain switching. The intrinsic electromechanical properties can be obtained by using the single crystal properties and a domain orientation distribution function (ODF) [40]. However, as the intrinsic effect only contributes about 30% to the dielectric and piezoelectric constants at room temperature [41], it is far from accuracy to simply take into account the intrinsic effect. In this 2D FE model, both the dielectric and elastic constants are assumed to be isotropic and not to vary with domain switching [20]. As to the piezoelectric constants, it is assumed that the properties are proportional to the remnant polarization [14,17]: jPr j s d d¼ P max
ð22Þ s
where Pr is the remnant polarization vector, P max and d are maximum polarization and piezoelectric tensor for the fully poled ceramics, respectively. 4.3. Boundary conditions Fig. 4 shows the boundary conditions and domain orientations in a 2D 4 5 rectangular FE mesh of a polycrystalline ferroelectric. Note that in this paper, only the uni-axial electromechanical load-
traction and charge free on lateral boundaries
In this paper, the change of Gibbs free energy of a specific domain before and after switching is considered as the driving force for domain switching. The Gibbs free energy (per unit volume) of a ferroelectric can be defined as [42]: 1 1 r gðr; EÞ ¼ r : es þ E Ps þ E d : r þ E k E þ r : SE : r 2 2 ð23Þ Bearing in mind the charge screening effect, the electromechanical field inside a grain (element) is assumed to be unaffected by the 180° switching within the grain. The 180° domain switching criterion can then be written as g bef g aft P W 180 bef
ð24Þ
aft
where g and g are Gibbs free energy of the specific domain before and after switching, respectively; and W 180 is the energy barrier (per unit volume) for 180° switching and is taken equal to 2P0 EC [14–20]. Since 180° switching is not constrained by neighboring domains, complete 180° switching occurs inside a grain once Eq. (24) is satisfied. As 90° switching is constrained by the neighboring domains and the degree of constraint is switching-fraction dependent, the 90° switching criterion is also switching-fraction dependent and is expressed as: f90 g bef ðE; rÞ f90 g aft ðE; r; f90 Þ P f90 W 90
ð25Þ
where f90 is the fraction of the completed 90° switching; and W 90 is the energy barrier (per unit volume) for 90° switching. pffiffiffi Referring to Huber et al. [10], W 90 ¼ 2P0 EC to activate 90° switching at coercive field during electric poling. As the constraint by 90° switching increases monotonously with f90 , it is reasonable to assume that f90 increases monotonously with the applied field/ stress. Thus for 90° switching, the problem is to solve the maximum f90 which can satisfy Eq. (25) under a specific electromechanical loading. In this study, f90 , is estimated by a dichotomy approximation. It should be noted that in this paper, the domain switching criterion is based on Berlincourt and Krueger’s definition [21] of
T = T0 and Φ = Φ 0 at the top
X2
y
θ
X1
u1 = 0
x u 2 = 0 and Φ = 0 at the bottom Fig. 4. Boundary conditions and domain orientations in the 2D FE model.
Fig. 5. 2D illustration of complete 180° domain switching based on Berlincourt and Krueger’s definition [21] in tetragonal ferroelectrics during electric poling.
327
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
complete 180° switching during electric poling. As shown in Fig. 5, during electric poling, complete 180° switching means that all the domains in Zone III switch to Zone I via 180° switching. Domains in Zone II and Zone IV can only switch to Zone I via 90° switching. No 180° switching can occur between domains in Zone II and Zone IV in this case. 4.5. FE Algorithm The analysis process of the 2D FE model is described as below. 1. Assemble the global stiffness matrix using the material constants. 2. At each step of loading, calculate the global traction matrix to get the electric field and stress inside each grain. 3. For a specific grain, calculate the Gibbs free energies for all the 4 types of domains to determine the specific type of domain which is energetically most favorable, say Domain 1 in Fig. 3. 4. Check if Domain 2 can switch to Domain 1 via 180 switching. 5. Check whether domain 3 or domain 4 has the higher Gibbs energy, say Domain 3. 6. Assume 0.5 fraction of domain 3 can switch to Domain 1 via 90° switching, then update the global traction matrix caused by 90° domain switching. Re-analyze the electric field/stress state of the specific grain and determine whether 0.5 fraction of 90° switching from Domain 2 to Domain 1 can be accomplished.
a
If it can, then the next iterative fraction of 90° switching is 0.75, otherwise, the next iterative switching-fraction is 0.25, and so on. The iterations are carried out for 10 steps and achievable fraction of 90 switching is obtained. 7. Analyze the next grain using the original electric field/stress, repeat step 2–5. 8. When the analysis of all the grains is finished, calculate the polarization and strain of the whole ceramics and update the global stiffness matrix for next step of loading. 9. Increase (or decrease) loading, repeat step 2–8. As a 10-step dichotomy approximation is used to obtain the fraction of achievable 90° domain switching at each step of loading (unloading), the computational complexity is about one order higher than the conventional FE models for ferroelectrics [20]. In the present study for a 15 10 finite element mesh with 160 loading steps, the required time on a P4 3.0G PC is about 73 min. 5. Results and discussion Using the 2D FE model, domain switching in polycrystalline BaTiO3 ceramics is studies under uni-axial electromechanical loading. The material constants used in the simulations are: P0 ¼ 0:26 C=m2 , c=a 1 ¼ 3c0 =2 ¼ 0:01, EC ¼ 0:4 MV=m, k=k0 ¼ 1900, sE ¼ 8:93 1012 m2 =N, d33 ¼ 191 pC=N, d31 ¼ 79 pC=N, d15 ¼ 270 pC=N.
0.15
0.10
0.10
0.10
0.05
0.05
0.05
0.00
0.00
-0.05
-0.05
-0.05
-0.10
-0.10
-0.10
-0.15
-0.15
-0.15
2
D3 (C/m )
0.15
2
D3 (C/m )
c 0.15
-1.0
-0.5
0.0
0.5
0.00
-1.0
1.0
-0.5
b
0.0
0.5
1.0
E3 (MV/m)
E3 (MV/m)
d
0.030
0.030
0.030
0.025
0.025
0.020
0.020
0.015
0.015
0.010
0.010
0.005
0.005
0.000
0.000
-0.005
-0.005
-0.005
-0.010
-0.010
0.025
-0.010 -1.0
-0.5
0.0
E3 (MV/m)
0.5
1.0
0.015
ε33 (%)
ε33 (%)
0.020
0.010 0.005 0.000
-1.0
-0.5
0.0
0.5
1.0
E3 (MV/m)
Fig. 6. Response of BaTiO3 ceramics under pure electric loading (a) D–E hysteresis loops, (b) longitudinal strain curves, (c) D–E hysteresis loops without 90° switching, and (d) longitudinal strain curves without 90° domain switching.
328
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
Fig. 6 shows the D–E hysteresis loops and longitudinal strain curves under a cyclic electric loading obtained from the present model with 90° domain switching. For comparison, the simulated hysteresis loops and strain curves in the absence of 90° domain
0.030
0.030
Switching strain Total strain
0.025
0.025 0.020
0.015
0.015
0.010
0.010
0.005
0.005
0.000
0.000
-0.005
-0.005
ε33 (%)
0.020
-0.010
-0.010 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m) Fig. 7. Stabilized strain-electric field of BaTiO3 ceramics under electric loading.
a 0.15
0.10
2
D3 (C/m )
0.05
0 MPa 25MPa 50MPa 100MPa
0.00
-0.05
-0.10
-0.15 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m)
b
0.010
0.010
0 MPa 0.005
Switching strain (%)
0.005
0.000
0.000
25 MPa -0.005
-0.005
-0.010
-0.010
50 MPa -0.015
-0.015
-0.020
-0.020
100 MPa -0.025
-0.025 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m) Fig. 8. Computed polarization and switching strain curves of BaTiO3 under combined electric and uniaxial compression loading.
switching are also shown. It can be seen from Fig. 6 that only after several cycles of electric loading the hysteresis loops and butterfly curves stabilize. This is due to the 10-step dichotomy approximation. In the case of complete 180° switching without 90° switching (see Fig. 6c), the resultant polarization is Pr ¼
4
R p=4 0
P0 cos hdh pffiffiffi ¼ 2P0 =p ¼ 0:45P 0 ¼ 0:117 C=m2 2p
ð26Þ
The remnant polarization in Fig. 6a is about 0.133 C/m2, thus it can be deduced from Fig. 6a and c that some 90° domain switching occurs in BaTiO3 ceramics during electric poling. To evaluate the fraction of 90° switching, the stabilized strainelectric field curve is shown in Fig. 7 where both the switching strain and total strain are present. The switchable strain in Fig. 7 is about 0.007%. As in the case of complete 90° switching, the resultant switching strain is: R p=4 2 c0 ð2 cos2 h sin hÞ=2 dh c0 =4 ¼ 3c0 =ð2pÞ cmax ¼ 0 p=4 ¼ ðc=a 1Þ=p ¼ 0:32%
ð27Þ
The fraction of achievable 90° switching is then estimated to be 0.007/0.32 = 2.2%. Fig. 8 shows the simulated polarization and strain curves of BaTiO3 ceramics under combined electric loading and uni-axial compression. With an increase in applied compressive stress, both the switchable polarization and remnant strain decrease. This is because some domains in Zone I and Zone III (Fig. 5) switch to and get locked in Zone II and Zone IV due to the applied compression. However, the compressive stress has little effect on the switchable strain even at 100 MPa, which may imply that the internal stress caused by 90° switching is much higher than this level. The predicted fraction of 90° switching in BaTiO3 ceramics is 2.2%, which is smaller than the experimental values (12% in [21] or 9% in [22]). This implies that plane strain assumption and the use of rectangular finite elements in the modeling may considerably overestimate the constraint by 90° domain switching and thus underestimate the achievable fractions of 90° switching. On the other hand, in a real ceramic, higher internal stresses may cause plastic deformations which in turn relax the stresses. Domain switching in tetragonal PZT-45/55 ceramics [43] is also considered. Fig. 9 shows the stabilized strain-electric field behaviour. The material constants used in the simulations are c=a 1 ¼ 3c0 =2 ¼ 0:0277, EC ¼ 1 MV=m, P0 ¼ 0:52 C=m2 , sE ¼ 11:7 1012 m2 =N. The piezoelectric constants are identical to the values corresponding to BaTiO3 ceramics as there are no piezoelectric properties reported for tetragonal PZT ceramics in the literature. Note that the piezoelectric constants have only a minor effect on the fractions of 90° switching. The switchable strain in Fig. 9 is about 0.017% and in the case of complete 90° switching it is 0.882% according to Eq. (27). Therefore, the fraction of achievable 90° switching in tetragonal PZT-45/55 ceramics is about 1.9% under electric loading up to 3 MV/m. Note that the simulation results agree qualitatively with the experimental results in [43]. It is clear from the above numerical simulations that the FE formulation of the present constrained domain switching model, which is based on a multi-domain grain assumption and estimates the fraction of 90° switching via a dichotomy approximation, is free from domain switching related numerical instability problems that may be encountered in conventional FE models for ferroelectrics. The present model can predict the fraction of 90° switching in polycrystalline ferroelectrics without using fitting parameters. The current 2D FE model can only be used for tetragonal ferroelectrics, but it can be extended to the 3D cases covering both tetragonal and rhombohedral ferroelectrics. However, the dichotomy approximation in the 2D case cannot be extended in a straightforward
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
0.08
0.08
Switching strain Total strain
Strain (%)
0.06
0.06
0.04
0.04
0.02
0.02
0.00
0.00
-0.02
-0.02
329
nal stress caused by non-180° switching is naturally considered by the FE model. The fraction of non-180° switching in each grain is obtained through a dichotomy approximation. The switching instability problem is thus well resolved in the proposed constrained domain switching model. Numerical simulations show that during electrical poling, 180° switching is usually complete while 90° switching is only about 2% for both BaTiO3 and tetragonal PZT ceramics, which agrees quite well with the existing experimental results. The FE model can be extended to the 3D cases to cover both tetragonal and rhombohedral ferroelectrics by using further assumptions. Acknowledgement
-3
-2
-1
0
1
2
3
This work was supported by a grant from Natural Science and Engineering Research Council of Canada (NSERC).
E3 (MV/m)
References Fig. 9. Computed strain-electric field response of tetragonal PZT-45/55 ceramics.
z 1
2 5
y 3 4
6
x Fig. 10. Six types of domains in 3D tetragonal ferroelectric crystals.
manner to the 3D case because two or more types of non-180° switching may occur simultaneously. As shown in Fig. 10, when an electric field is applied close to z axis of the crystallite coordinates in tetragonal ferroelectrics, there may exists two types of 90° switching (say from domain 6 and domain 4 to domain 1) simultaneously. But if it is further assumed that domain 4 and domain 6 always switch together at equal fractions to domain 1, then the dichotomy approximation can still be applied to obtain the fraction of non-180° switching. Obviously, the model can also be extended to rhombohedral materials if such an assumption is made. 6. Conclusions A nonlinear finite element (FE) formulation for ferroelectric materials is briefly reviewed and it is illustrated that for FE models based on a phenomenological material model with the electric potential as a nodal variable, the Newton–Raphson iteration is usually unstable. While for FE models based on a micromechanical domain switching criterion, the large depolarization field/stress caused by domain switching makes the domains switch back to its original state, thus switching instability arises. To overcome these domain switching related instabilities and to predict the experimentally observed incomplete non-180° switching in ferroelectric ceramics, a constrained domain switching model is proposed and implemented in a FE code. A multi-domain grain assumption is used and charge screening effect is also taken into account in the present model. The depolarization field by polarization switching is thus compensated by free charges and the inter-
[1] Y. Xu, Ferroelectric Materials and their Applications, Amsterdam Press, NorthHolland, 1991. [2] B. Jaffe, W.R. Cook, H. Jaffe, Piezoelectric Ceramics, Academic Press, London and New York, 1971. [3] H. Allik, T.J.R. Hughes, Int. J. Num. Meth. Eng. 2 (1970) 151–157. [4] C.L. Hom, N. Shankar, Int. J. Solids Struct. 33 (1996) 1757–1779. [5] X. Gong, Z. Suo, J. Mech. Phys. Solids 44 (1996) 751–769. [6] K. Ghandi, N.W. Hagood, Smart structures and materials: mathematics and control in smart structures, in: V.V. Varadan (Ed.), Proceedings of the SPIE 2715 (1996) 121–140. [7] K. Ghandi, N.W. Hagood, Smart structures and materials: mathematics and control in smart structures, in: V.V. Varadan, J. Chandra (Eds.), Proceedings of the SPIE 3039 (1997) 97–112. [8] M. Kamlah, Ü. Böhle, Int. J. Solids. Struct. 38 (2001) 605–633. [9] M. Kamlah, A.C. Liskowsky, R.M. McMeeking, H. Balke, Int. J. Solids Struct. 42 (2005) 2949–2964. [10] J.E. Huber, N.A. Fleck, J. Mech. Phys. Solids 49 (2001) 785–811. [11] W. Chen, C.S. Lynch, J. Eng. Mater. Tech. 123 (2001) 169–175. [12] T. Felt, D. Munz, G. Thun, J. Am. Ceram. Soc. 86 (2003) 1427–1429. [13] F.X. Li, D.N. Fang, J. Intell. Mater. Syst. Struct. 16 (2005) 583–588. [14] S.C. Hwang, C.S. Lynch, R.M. McMeeking, Acta Metall. Mater. 43 (1995) 2073– 2084. [15] S.C. Hwang, R.M. McMeeking, Ferroelectrics 211 (1998) 177–194. [16] S.C. Hwang, R.M. McMeeking, Int. J. Solids Struct. 36 (1999) 1541–1556. [17] W. Chen, C.S. Lynch, Eng. Fract. Mech. 64 (1999) 539–562. [18] T. Steinkopff, J. Eur. Ceram. Soc. 19 (1999) 1247–1249. [19] R.M. McMeeking, S.C. Hwang, Ferroelectrics 200 (1997) 151–173. [20] F.X. Li, D.N. Fang, Mech. Mater. 36 (2004) 959–973. [21] D. Berlincourt, H.H.A. Krueger, J. Appl. Phys. 30 (1959) 1804–1810. [22] E.C. Subbarao, M.C. McQuarrie, W.R. Buessem, J. Appl. Phys. 28 (1957) 1194– 1200. [23] O. Axelsson, Iterative Solution Methods, Cambridge University Press, Cambridge, 1994. [24] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, vol. 1, Wiley, Chichester, 1991. [25] C.M. Landis, Int. J. Num. Meth. Eng. 55 (2002) 613–628. [26] M. Kamlah, Z. Wang, Comput. Mater. Sci. 28 (2003) 409–418. [27] C.M. Landis, J. Appl. Mech.-T ASME 70 (2003) 470–478. [28] C.M. Landis, Current Opin. Solid State Mater. Sci. 8 (2004) 59–69. [29] V. Mehling, Ch. Tsakmakis, D. Gross, J. Mech. Phys. Solid 55 (2007) 2106–2141. [30] S.C. Hwang, J.E. Huber, R.M. McMeeking, N.A. Fleck, J. Appl. Phys. 83 (1998) 1530–1540. [31] H. Kessler, H. Balke, J. Mech. Phys. Solids. 49 (2001) 953–978. [32] H. Kessler Jr., E.R. Fuller, H. Balke, Smart structures and materials: behavior and mechanics of active materials, in: C.S. Lynch (Ed.), Proceedings of the SPIE 3992 (2000) 234–244. [33] X.W. Zhang, C. Lei, K.P. Chen, J. Am. Ceram. Soc. 88 (2005) 335–338. [34] M.E. Lines, A.M. Glass, Principles and Applications of Ferroelectrics and Related Materials, Clarendon Press, Oxford, 1977. [35] D.C. Lupascu, Fatigue in Ferroelectric Ceramics and Related Issues, SpringerVerlag Berlin Heidelberg, Berlin-New York, 2004. [36] G.I. Taylor, J. Inst. Metals 62 (1938) 307–324. [37] K. Bhattacharya, R.V. Kohn, Acta Mater. 44 (1996) 529–542. [38] J.Y. Li, R.C. Rogan, E. Üstündag, K. Bhattacharya, Nature Mater. 4 (2005) 776– 781. [39] S. Wan, K.J. Bowman, J. Mater. Res. 16 (2001) 2306–2313. [40] F.X. Li, R.K.N.D. Rajapakse, J. Appl. Phys. 101 (2007) 054110. [41] Q.M. Zhang, H. Wang, N. Kim, L.E. Cross, J. Appl. Phys. 75 (1994) 454–459. [42] W. Lu, D.N. Fang, K.C. Hwang, Acta Mater. 47 (1999) 2913–2926. [43] M.J. Hoffmann, M. Hammer, A. Endriss, D.C. Lupascu, Acta Mater. 49 (2001) 1301–1310.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Nonlinear finite element modeling of polycrystalline ferroelectrics based on constrained domain switching F.X. Li, R.K.N.D. Rajapakse * Department of Mechanical Engineering, The University of British Columbia, 2054-6250 Applied Science Lane, Vancouver, Canada V6T 1Z4
a r t i c l e
i n f o
Article history: Received 9 December 2007 Received in revised form 16 March 2008 Accepted 25 March 2008 Available online 21 May 2008 PACS: 46.15.Cc 77.80.Fm 77.84.s Keywords: Constitutive models Finite element method Ferroelectrics Domain switching Charge screening
a b s t r a c t A review of the nonlinear finite element formulations for ferroelectrics shows that instability related to domain switching may arise in formulations using a phenomenological constitutive model or micromechanical domain switching criterion. To overcome the instability problems and to predict the experimentally observed incomplete non-180° switching in ferroelectrics, a constrained domain switching model for polycrystalline ferroelectrics is proposed in this paper and implemented in a finite element code. In the proposed model, a ferroelectric polycrystalline is made up of numerous random oriented grains, each of which contains multi-domains and is represented by a finite element. Charge screening effect in real ferroelectrics is also taken into account in the model, thus 180° domain switching is not constrained. While the internal stress induced by non-180° switching cannot be compensated in a similar manner and is considered naturally by the finite element method. The fraction of non-180° switching in each grain is obtained through a dichotomy approximation. Selected numerical results show that during electrical poling, 180° switching is usually complete while 90° switching is only about 2% for both BaTiO3 and tetragonal PZT ceramics, which agree well with the experimental results. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction Ferroelectric materials such as Barium Titanate and PZT ceramics are widely used in smart structures and devices as actuators, sensors, ultrasonic motors, etc. [1]. Domain switching is now accepted as the predominant source of the nonlinear response of these ceramics [2]. Ferroelectric ceramics are also very brittle and are subjected to complex electromechanical loading in practical applications that may cause domain switching and cracking. It is therefore important to develop nonlinear finite element models to analyze ferroelectric materials under complex loading. Since the original work of Allik and Hughes [3], linear finite element analysis of piezoelectric materials has received considerable attention. Hom and Shankar [4], and Gong and Suo [5] independently presented nonlinear finite element formulations to study multi-layer actuators made of relaxor ferroelectric ceramics. A fully coupled nonlinear finite element model for phase transition of electromechanically coupled materials was proposed by Ghandi and Hagood [6] by using nodal displacements and electric potential as degrees of freedoms. They later presented a finite element model based on the Hu–Washizu functional, taking the electric displace* Corresponding author. Present Address: Institute for Computing, Information and Cognitive Systems, The University of British Columbia, Canada. E-mail address: [email protected] (R.K.N.D. Rajapakse). 0927-0256/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2008.03.040
ment as nodal variables [7]. Recently, Kamlah and Bohle [8] presented a nonlinear finite element model based on a bilinear phenomenological constitutive model. They first solved a pure dielectric boundary-value problem for the history of the electric potential and thereafter prescribing this electric potential, an electromechanical stress analysis with the mechanical boundary conditions was conducted. Kamlah et al. [9] also conducted finite element simulations by using a multi-domain single crystal switching model. The above mentioned nonlinear finite element (FE) models are based on explicit or implicit phenomenological constitutive models. There are only few experimental studies for ferroelectric materials under multi-axial loading [10–13], and the FE models are therefore often used to simulate one-dimensional problems for the purpose of comparison with experiments [8]. Since the development of the first domain switching criterion for ferroelectric materials by Hwang et al. [14], several nonlinear FE models accounting for domain switching have appeared in the literature [15–18]. Hwang et al. [15,16] formulated two nonlinear FE models for ferroelectric and ferroelastic materials by using an improved domain switching criterion based on potential energy [19]. To simplify the modeling, they solved the mechanical and dielectric boundary-value problems separately, neglecting the piezoelectric effect. Recently, Li and Fang [20] presented an improved domain switching criterion based on the Gibbs free energy and
323
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
developed a fully coupled, three-dimensional (3D) nonlinear FE model to study the effects of intergranular interactions during domain switching. The nonlinear FE models accounting for domain switching usually allow only one domain to switch in each iteration step [15,16,20], thus leading to quite low numerical efficiency. Another problem with the existing FE models is that the energy barrier for 90° (or non-180°) domain switching is not well known, though the energy barrier for 180° switching is widely accepted to be equal to 2P0EC (P0 is the magnitude of the spontaneous polarization and EC the magnitude of coercive field). This paper first examines the so called ‘‘switching-related instability” problem which is likely to be encountered in the FE models accounting for domain switching. A general nonlinear finite element formulation is summarized in Section 2. In Section 3, it is shown that for the nonlinear FE models based on a phenomenological constitutive law, the Newton–Raphson iteration is usually unstable if the electric potential is taken as a nodal variable. In the case of FE models using a micromechanical domain switching criterion, the eigen changes in spontaneous polarization/strain caused by domain switching are found to be equivalent to a rather large depolarization field/stress on the switched domain, which makes the domain switch back to its original state, thus causing switching instability. To overcome these problems and to predict the observed incomplete non-180° domain switching in ferroelectrics [21,22], a constrained domain switching model and its 2D finite element implementation are presented in Section 4. Selected numerical results are presented in Section 5 and the conclusions are in Section 6.
To examine the causes of domain switching related instability, a general nonlinear finite element formulation is summarized in this section [20]. For a ferroelectric material with a remnant polarization and remnant strains, the constitutive equations can be written as
Dj ¼ ejmn ðemn ermn Þ þ
e kjm Em
þ
ð1Þ Prj
ð2Þ
Di;i ¼ 0
ð3Þ
Let ui and / denote the components of the displacement vector and electric potential, respectively. The strains and electric field can be expressed as eij ¼
1 ðui;j þ uj;i Þ; 2
Ei ¼
o/ oxi
uz
/ T
ð9Þ
The corresponding generalized stress, strain and prescribed traction vectors are defined as R ¼ ½rx ; ry ; rz ; rxy ; rxz ; ryz ; Dx ; Dy ; Dz T
ð10Þ
C ¼ ½ex ; ey ; ez ; 2exy ; 2exz ; 2eyz ; Ex ; Ey ; Ez
T
ð11Þ
T ¼ ½tx ; ty ; tz ; pe T
ð12Þ
The constitutive law can be written in the following incremental form. dR ¼ C dC þ DCðC Cr Þ þ dRr
ð13Þ
where C is the generalized constitutive matrix which includes elastic, piezoelectric and dielectric coefficients. Following incremental form of displacement interpolation can be defined for a typical finite element. Du ¼ N Dae
ð14Þ e
where Du; Da and N denote the incremental generalized displacement vector of an element, incremental generalized nodal displacement vector and the interpolation matrix, respectively. The incremental generalized strain vector of an element is given by ð15Þ
where L is the linear differential matrix operator and B = LN. The incremental global finite element equations can be expressed as K Dae ¼ DP
ð16Þ
where X T e K¼ G K G
ð17Þ
e
where rij , emn are the components of the stress tensor and the strain tensor, respectively; Em and Dj are the components of the electric field vector and the electric displacement vector, respectively; e C Eijmn , ejmn and kjm are the elements of the elastic stiffness tensor, piezoelectric tensor and the dielectric permittivity tensor, respectively; and ermn and Prj are the components of the remnant strain tensor and remnant polarization vector, respectively. Consider a piezoelectric body that occupies a volume, V, with a closed surface S. If body forces and body charges can be neglected, the equilibrium equations can be written as, rij;j ¼ 0;
uy
u ¼ ½ ux
DC ¼ LDu ¼ LNDae ¼ BDae
2. Nonlinear finite element formulation
rij ¼ C Eijmn ðemn ermn Þ emij Em
kk denotes the jump across a boundary; pe is the surface charge per unit area; ui denotes the components of the prescribed displacement; and / is the prescribed electric potential. The generalized displacement vector of a finite element is taken as
ð4Þ
The boundary conditions can be expressed as rij nj ¼ ti on St
ð5Þ
kni Di k ¼ pe on SD
ð6Þ
ui ¼ ui on Su
ð7Þ
/ ¼ / on S/
ð8Þ
where nj are the components of the outward unit normal vector of the surface S and t i denotes the components of the traction vector;
and Ke ¼
Z
BT CBdve
ð18Þ
ve
DP ¼
Z
NT DTds þ St þSD
Z
X e
GT
Z
BT CDCr dve þ ve
BT DCðC Cr Þdve ve
X e
GT
Z
X
GT
e
BT DRr dve
ð19Þ
ve
where G is a coordinate transformation matrix that relates generalized element stiffness in local coordinates to the global coordinates and the generalized element nodal vector in local coordinates to the global coordinates. It should be noted that for ferroelectric materials during domain switching, the second term on the right side of Eq. (19) corresponds to the contribution due to the change of the remnant strain and the third term corresponds to the contribution due to the change of the constitutive matrix C, i.e., change of materials properties. The last term of Eq. (19) represents the contribution due to the change of remnant polarization. 3. Domain switching related instabilities in FE modeling of ferroelectrics In this section, the iteration instabilities which might be encountered in nonlinear FE analysis based on phenomenological
324
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
constitutive models and the ‘‘back-switching” phenomenon in nonlinear FE models using a micromechanics based domain switching criterion are discussed. 3.1. Switching induced iteration instability in FE models based on phenomenological constitutive laws For FE models based on phenomenological material models, the material response to applied fields, such as the D–E hysteresis loop, stress–strain curve, can be obtained directly from the material models, explicitly or implicitly, without the help of the FE method. Domain switching (or phase transitions) is not involved in the nonlinear FE formulation, i.e., the second and the last terms on the right hand side of Eq. (19) vanish. The FE formulation of such models is similar to the nonlinear FE formulation of conventional elastoplastic materials except for the additional terms involving dielectric and piezoelectric effects. An incremental nonlinear iterative process is needed in the FE analysis due to the history-dependent behavior of ferroelectric materials. Among the nonlinear iterative methods, the Newton– Raphson method is common and convenient with good accuracy and convergence [23]. Note that in conventional elastoplastic analysis, the Newton–Raphson method is usually stable [24]. However, for ferroelectric materials due to the characteristics of domain switching, the slope of D–E curves near the coercive field is much larger than that of the linear dielectric part. The Newton–Raphson iteration can become unstable due to such abrupt changes of slope. To illustrate this iterative instability, the poling curve of a typical PZT ceramic is shown in Fig. 1. For 1D dielectric problems with the electric potential as the nodal unknowns with prescribed electric displacement D0 , it can be seen from Fig. 1a that Newton– Raphson iteration to solve for the exact solution Ex is unstable, i.e., Ei ! Eiþ1 ! divergence, where Ei and Eiþ1 are the electric field computed at the ith and i + 1th iteration. Landis [25] indicated that by taking the vector electric displacement instead of the scalar electric potential as nodal unknowns, a standard iteration process is stable for ferroelectric materials with a simple bilinear D–E response. However, if the saturation part of the poling curve (BC in Fig. 1b) is also taken into account, even using the electric displacement as nodal unknowns, the Newton–Raphson iterative process for a prescribed electric field may still be unstable. As shown in Fig. 1b, Di is the ith iteration and it is obvious that Di ! Diþ1 ! divergence. This can be true even if the ith iteration, such as Di1 in Fig. 1b, is very close to the exact solution Dx , that ! divergence, and Di ! Diþ1 ! divergence. is Di1 ! Diþ1 1 The above iteration instabilities are caused by the steep slope of the curve AB in Fig. 1a and the curve BC in Fig. 1b associated with domain switching. The D–E curve of a non-polar dielectric material without switching is exactly like the stress–strain curve of a con-
a
ventional elastoplastic material, and the Newton–Raphson iteration is always stable [24]. Therefore, this type of iteration instability is called ‘‘domain switching related iteration instability” to distinguish it from the direct instability of the solution due to domain switching that is discussed in the ensuing section. It should be noted that such ‘‘iteration instabilities” are usually encountered in incremental FE models using phenomenological constitutive laws. It arises from the iterative process when solving for the unknown electric field (electric displacement) with the prescribed electric displacement (electric field). In phenomenological models using numerical methods other than FEM [26–28], the iterations may not be involved in the solution process and above mentioned ‘‘iteration instabilities” may not exist. Moreover, even in FE models using the iterative processes, if the nonlinear items, such as the remnant polarization and strain, are separately taken into account as internal variables [29], ‘‘iteration instabilities” can be avoided by assuming that the linear properties do not change with the internal variables or by using a simplified bilinear D–E assumption [8]. 3.2. Back-switching in FE models utilizing micromechanics based domain switching criteria In the case of nonlinear FE models utilizing micromechanics based domain switching criteria [15–18,20], the equivalent nodal load induced by the change of spontaneous strain and spontaneous polarization during domain switching, i.e., the second and the last terms of the right hand side of Eq. (19), can be rather large compared to the contribution due to moderate electromechanical loading. Consider the case of 180° domain switching without a change in spontaneous strain (Fig. 2a). When the polarization of a domain reverses, a large depolarization field Ed is generated. In the absence of free charge compensation, the depolarization field Ed is on the order of DP=k (10–50 kV/mm), where DP is the magnitude of the change of spontaneous polarization during switching and k is the dielectric constant. Such a large depolarization field (usually higher than 10 times the coercive field) makes the domain switches back to its original state, thus domain switching instability arises. In the case of 90° domain switching (Fig. 2b), in addition to the large depolarization field Ed caused by the spontaneous polarization change, a large mechanical field (tensile stress r1 and compressive stress r2 ) is also generated due to the spontaneous strain change during domain switching. This mechanical field also tends to make the switched domain reverts to its original strain state. r1 and r2 are both of the order of Dc=s (1 GPa, where Dc is the magnitude of the spontaneous strain change and s is the elastic compliance), so that ‘‘back-switching” takes place. Therefore in an unmodified nonlinear FE formulation, domain switching is always unstable due to the terms involved in Eq. (19).
b
D
E
C
C
E
B
D
0
0
B A
A O
i
E
Ex
E
i+1
E
O Di
D
i+1
i D1 Dx D1i+1
D
Fig. 1. Illustration of Newton–Raphson iteration instability: (a) electric potential as nodal unknowns; (b) electric displacement as nodal unknowns.
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
Fig. 2. Depolarization electric field and mechanical field induced by (a) 180° switching or ferroelectric switching and (b) 90° switching or ferroelastic switching.
Domain switching instability associated with the finite element formulation has been noticed in previous modeling studies [16,17,20,25]. To weaken the high electric and mechanical fields induced by ferroelectric and ferroelastic switching in an inclusion model which is equivalent to FE models, Hwang et al. [30] employed a much smaller Young’s modulus (about 20% of the true value) and a rather larger dielectric constant (about 14 times the true value) in their simulations. Chen and Lynch [17] noticed in their modeling, domain switching instability caused by the high local electric field due to the jump of polarization during switching. To reduce this effect, they introduced a factor kð0 < k < 1Þ to the last term of the right hand side of Eq. (19). Kessler et al. [31,32] discussed the unstable switching at the initiation of domain switching in their numerical simulations (not FE model). They suggested that switching instability can be removed by introducing ferroelectric ‘‘hardening” in the material model, which is equivalent to the simplified bilinear D–E curve assumption of Kamlah and Bohle [8]. However, the switching instability addressed in their work is associated with a ‘‘jump” in polarization or strain which is different from the ‘‘back-switching” addressed in this paper. 4. FE model based on constrained domain switching In this section, a constrained domain switching model free from the above mentioned instabilities in the finite element analysis is proposed and implemented in a finite element code. Furthermore, the present model can also predict the experimentally observed incomplete non-180° domain switching in ferroelectric ceramics [21,22,33]. As discussed above, the cause of instability in FE models based on domain switching criteria is the rather large depolarization field/stress induced by the domain switching associated terms. In real ferroelectric ceramics, the large depolarization field by mismatched polarization is fully or partially compensated by free charges inside the material, which is called ‘‘charge screening” [34,35]. The internal stress caused by non-180° switching, however, cannot be compensated in a similar way. The non-180° switching induced constraint in polycrystalline ferroelectrics are actually crystal symmetry dependent. According to the Taylor rule of plasticity [36], a crystal must have at least five slip systems for a
325
polycrystalline to be ductile. In terms of the deformation modes, a ferroelectric crystal is similar to a ferroelastic crystal or a shape memory alloy [37] which has two independent slip systems for the tetragonal and three for the rhombohedral case. Thus in tetragonal or rhombohedral polycrystalline ferroelectric ceramics, non180° domain switching generates large internal stresses, i.e., they are constrained by neighboring grains and may not occur. While in orthorhombic ferroelectric ceramics or PZT ceramics near the MPB, where tetragonal and rhombohedral phases coexist, there are five or more independent deformation modes, which make the materials ductile. Then non-180° domain switching does not generate very large internal stresses from neighboring grains and such switching can almost be accomplished completely in these types of polycrystallines [37,38]. In the present study, we do not study these ‘‘ductile” polycrystallines but focus only on the constrained domain switching phenomena in tetragonal and rhombohedral ferroelectric ceramics. So far, only a few theoretical studies have addressed the experimentally observed incomplete non-180° domain switching [21,22,33] except the work by Li et al. [38] and Wan et al. [39]. Li et al. [38] deduced through a rigorous mathematical analysis that the lower limit of switching strains that can be achieved in tetragonal and rhombohedral ferroelectric ceramics. In the phenomenological model of Wan and Bowman [39], they used an adjustable parameter to address the depolarization field by 90° switching and fit the model to experimental results. As to the FE models, the models using phenomenological material models as well as those based on domain switching criteria with the single-domain grain assumption cannot address the incomplete non-180° switching. The FE simulation by Kamlah et al. [9], which is based on a multi-domain single crystal switching model, could technically predict the fraction of non-180° switching, but they did not specifically address this issue. In the constrained domain switching model proposed in this paper, a ferroelectric polycrystalline is made up of numerous random oriented grains, each of which contains N types of domains (N = 6 for the tetragonal ferroelectrics and N = 8 for the rhombohedral) and is represented by a finite element. In an unpoled ceramic, the fraction of each type of domains in a grain is equal to 1/N, thus both the remnant polarization and strain of each grain is zero. Inside each grain, the electromechanical field is assumed to be uniform. Charge screening effect is taken into account in this model, i.e., the mismatched polarization by domain switching is to be compensated by free charges. Thus, an 180° switching does not generate a depolarization field and not get constrained. The large depolarization stress caused by the spontaneous strain change during non-180° switching cannot be compensated in a similar way and is taken into account by the FE method naturally. As the depolarization stress increases with the fraction of non-180° switching, the achievable fraction of non-180° switching increases with the applied field/stress and is approached by a mathematical approximation. In this model, it is also assumed that a ferroelectric polycrystalline is linear dielectric, elastic and piezoelectric unless domain switching occurs, i.e., all nonlinear polarization and strains are caused by domain switching.
4.1. A simplified 2D finite element model In this paper, the constrained domain switching model is implemented in a 2D finite element code. In the 2D case of tetragonal materials, there only exist four types of domains inside each grain, as shown in Fig. 3a. For a textureless polycrystalline ferroelectrics, the orientations of the grains, which can be denoted by an orientation angle h (0 6 h < 2p) in Fig. 3b, distribute randomly in all directions, i.e., h distributes uniformly in ½0; 2p. Note that in the 2D case,
326
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
a
b
y a
y
a c
1 4
ing cases are considered without loss of any generality. At the bottom face, both the vertical displacement u2 and the electric potential / are fixed for all the nodes. The horizontal displacement of the left node on the bottom face is also fixed to prevent rigid body motion. At the top face, both the electric potential and tractions are prescribed according to the applied electromechanical loading. Both lateral boundaries are traction and charge free.
X2
4
1
θ
x 2
X1 4.4. Domain switching criterion
3
2
3
x
Fig. 3. Domain and grain orientations for tetragonal polycrystalline ferroelectrics. (a) domain orientation in local grain coordinates x–y; (b) grain orientation in sample coordinates X1–X2.
the remnant strain for an unpoled polycrystalline ferroelectric is not zero, but exx ¼ eyy ¼ ðc=a 1Þ=6 ¼ c0 =4
ð21Þ
4.2. Material constants during domain switching Theoretically, all the material constants of a ferroelectric polycrystalline change either substantially or slightly with domain switching. The intrinsic electromechanical properties can be obtained by using the single crystal properties and a domain orientation distribution function (ODF) [40]. However, as the intrinsic effect only contributes about 30% to the dielectric and piezoelectric constants at room temperature [41], it is far from accuracy to simply take into account the intrinsic effect. In this 2D FE model, both the dielectric and elastic constants are assumed to be isotropic and not to vary with domain switching [20]. As to the piezoelectric constants, it is assumed that the properties are proportional to the remnant polarization [14,17]: jPr j s d d¼ P max
ð22Þ s
where Pr is the remnant polarization vector, P max and d are maximum polarization and piezoelectric tensor for the fully poled ceramics, respectively. 4.3. Boundary conditions Fig. 4 shows the boundary conditions and domain orientations in a 2D 4 5 rectangular FE mesh of a polycrystalline ferroelectric. Note that in this paper, only the uni-axial electromechanical load-
traction and charge free on lateral boundaries
In this paper, the change of Gibbs free energy of a specific domain before and after switching is considered as the driving force for domain switching. The Gibbs free energy (per unit volume) of a ferroelectric can be defined as [42]: 1 1 r gðr; EÞ ¼ r : es þ E Ps þ E d : r þ E k E þ r : SE : r 2 2 ð23Þ Bearing in mind the charge screening effect, the electromechanical field inside a grain (element) is assumed to be unaffected by the 180° switching within the grain. The 180° domain switching criterion can then be written as g bef g aft P W 180 bef
ð24Þ
aft
where g and g are Gibbs free energy of the specific domain before and after switching, respectively; and W 180 is the energy barrier (per unit volume) for 180° switching and is taken equal to 2P0 EC [14–20]. Since 180° switching is not constrained by neighboring domains, complete 180° switching occurs inside a grain once Eq. (24) is satisfied. As 90° switching is constrained by the neighboring domains and the degree of constraint is switching-fraction dependent, the 90° switching criterion is also switching-fraction dependent and is expressed as: f90 g bef ðE; rÞ f90 g aft ðE; r; f90 Þ P f90 W 90
ð25Þ
where f90 is the fraction of the completed 90° switching; and W 90 is the energy barrier (per unit volume) for 90° switching. pffiffiffi Referring to Huber et al. [10], W 90 ¼ 2P0 EC to activate 90° switching at coercive field during electric poling. As the constraint by 90° switching increases monotonously with f90 , it is reasonable to assume that f90 increases monotonously with the applied field/ stress. Thus for 90° switching, the problem is to solve the maximum f90 which can satisfy Eq. (25) under a specific electromechanical loading. In this study, f90 , is estimated by a dichotomy approximation. It should be noted that in this paper, the domain switching criterion is based on Berlincourt and Krueger’s definition [21] of
T = T0 and Φ = Φ 0 at the top
X2
y
θ
X1
u1 = 0
x u 2 = 0 and Φ = 0 at the bottom Fig. 4. Boundary conditions and domain orientations in the 2D FE model.
Fig. 5. 2D illustration of complete 180° domain switching based on Berlincourt and Krueger’s definition [21] in tetragonal ferroelectrics during electric poling.
327
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
complete 180° switching during electric poling. As shown in Fig. 5, during electric poling, complete 180° switching means that all the domains in Zone III switch to Zone I via 180° switching. Domains in Zone II and Zone IV can only switch to Zone I via 90° switching. No 180° switching can occur between domains in Zone II and Zone IV in this case. 4.5. FE Algorithm The analysis process of the 2D FE model is described as below. 1. Assemble the global stiffness matrix using the material constants. 2. At each step of loading, calculate the global traction matrix to get the electric field and stress inside each grain. 3. For a specific grain, calculate the Gibbs free energies for all the 4 types of domains to determine the specific type of domain which is energetically most favorable, say Domain 1 in Fig. 3. 4. Check if Domain 2 can switch to Domain 1 via 180 switching. 5. Check whether domain 3 or domain 4 has the higher Gibbs energy, say Domain 3. 6. Assume 0.5 fraction of domain 3 can switch to Domain 1 via 90° switching, then update the global traction matrix caused by 90° domain switching. Re-analyze the electric field/stress state of the specific grain and determine whether 0.5 fraction of 90° switching from Domain 2 to Domain 1 can be accomplished.
a
If it can, then the next iterative fraction of 90° switching is 0.75, otherwise, the next iterative switching-fraction is 0.25, and so on. The iterations are carried out for 10 steps and achievable fraction of 90 switching is obtained. 7. Analyze the next grain using the original electric field/stress, repeat step 2–5. 8. When the analysis of all the grains is finished, calculate the polarization and strain of the whole ceramics and update the global stiffness matrix for next step of loading. 9. Increase (or decrease) loading, repeat step 2–8. As a 10-step dichotomy approximation is used to obtain the fraction of achievable 90° domain switching at each step of loading (unloading), the computational complexity is about one order higher than the conventional FE models for ferroelectrics [20]. In the present study for a 15 10 finite element mesh with 160 loading steps, the required time on a P4 3.0G PC is about 73 min. 5. Results and discussion Using the 2D FE model, domain switching in polycrystalline BaTiO3 ceramics is studies under uni-axial electromechanical loading. The material constants used in the simulations are: P0 ¼ 0:26 C=m2 , c=a 1 ¼ 3c0 =2 ¼ 0:01, EC ¼ 0:4 MV=m, k=k0 ¼ 1900, sE ¼ 8:93 1012 m2 =N, d33 ¼ 191 pC=N, d31 ¼ 79 pC=N, d15 ¼ 270 pC=N.
0.15
0.10
0.10
0.10
0.05
0.05
0.05
0.00
0.00
-0.05
-0.05
-0.05
-0.10
-0.10
-0.10
-0.15
-0.15
-0.15
2
D3 (C/m )
0.15
2
D3 (C/m )
c 0.15
-1.0
-0.5
0.0
0.5
0.00
-1.0
1.0
-0.5
b
0.0
0.5
1.0
E3 (MV/m)
E3 (MV/m)
d
0.030
0.030
0.030
0.025
0.025
0.020
0.020
0.015
0.015
0.010
0.010
0.005
0.005
0.000
0.000
-0.005
-0.005
-0.005
-0.010
-0.010
0.025
-0.010 -1.0
-0.5
0.0
E3 (MV/m)
0.5
1.0
0.015
ε33 (%)
ε33 (%)
0.020
0.010 0.005 0.000
-1.0
-0.5
0.0
0.5
1.0
E3 (MV/m)
Fig. 6. Response of BaTiO3 ceramics under pure electric loading (a) D–E hysteresis loops, (b) longitudinal strain curves, (c) D–E hysteresis loops without 90° switching, and (d) longitudinal strain curves without 90° domain switching.
328
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
Fig. 6 shows the D–E hysteresis loops and longitudinal strain curves under a cyclic electric loading obtained from the present model with 90° domain switching. For comparison, the simulated hysteresis loops and strain curves in the absence of 90° domain
0.030
0.030
Switching strain Total strain
0.025
0.025 0.020
0.015
0.015
0.010
0.010
0.005
0.005
0.000
0.000
-0.005
-0.005
ε33 (%)
0.020
-0.010
-0.010 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m) Fig. 7. Stabilized strain-electric field of BaTiO3 ceramics under electric loading.
a 0.15
0.10
2
D3 (C/m )
0.05
0 MPa 25MPa 50MPa 100MPa
0.00
-0.05
-0.10
-0.15 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m)
b
0.010
0.010
0 MPa 0.005
Switching strain (%)
0.005
0.000
0.000
25 MPa -0.005
-0.005
-0.010
-0.010
50 MPa -0.015
-0.015
-0.020
-0.020
100 MPa -0.025
-0.025 -1.0
-0.5
0.0
0.5
1.0
E3 (MV/m) Fig. 8. Computed polarization and switching strain curves of BaTiO3 under combined electric and uniaxial compression loading.
switching are also shown. It can be seen from Fig. 6 that only after several cycles of electric loading the hysteresis loops and butterfly curves stabilize. This is due to the 10-step dichotomy approximation. In the case of complete 180° switching without 90° switching (see Fig. 6c), the resultant polarization is Pr ¼
4
R p=4 0
P0 cos hdh pffiffiffi ¼ 2P0 =p ¼ 0:45P 0 ¼ 0:117 C=m2 2p
ð26Þ
The remnant polarization in Fig. 6a is about 0.133 C/m2, thus it can be deduced from Fig. 6a and c that some 90° domain switching occurs in BaTiO3 ceramics during electric poling. To evaluate the fraction of 90° switching, the stabilized strainelectric field curve is shown in Fig. 7 where both the switching strain and total strain are present. The switchable strain in Fig. 7 is about 0.007%. As in the case of complete 90° switching, the resultant switching strain is: R p=4 2 c0 ð2 cos2 h sin hÞ=2 dh c0 =4 ¼ 3c0 =ð2pÞ cmax ¼ 0 p=4 ¼ ðc=a 1Þ=p ¼ 0:32%
ð27Þ
The fraction of achievable 90° switching is then estimated to be 0.007/0.32 = 2.2%. Fig. 8 shows the simulated polarization and strain curves of BaTiO3 ceramics under combined electric loading and uni-axial compression. With an increase in applied compressive stress, both the switchable polarization and remnant strain decrease. This is because some domains in Zone I and Zone III (Fig. 5) switch to and get locked in Zone II and Zone IV due to the applied compression. However, the compressive stress has little effect on the switchable strain even at 100 MPa, which may imply that the internal stress caused by 90° switching is much higher than this level. The predicted fraction of 90° switching in BaTiO3 ceramics is 2.2%, which is smaller than the experimental values (12% in [21] or 9% in [22]). This implies that plane strain assumption and the use of rectangular finite elements in the modeling may considerably overestimate the constraint by 90° domain switching and thus underestimate the achievable fractions of 90° switching. On the other hand, in a real ceramic, higher internal stresses may cause plastic deformations which in turn relax the stresses. Domain switching in tetragonal PZT-45/55 ceramics [43] is also considered. Fig. 9 shows the stabilized strain-electric field behaviour. The material constants used in the simulations are c=a 1 ¼ 3c0 =2 ¼ 0:0277, EC ¼ 1 MV=m, P0 ¼ 0:52 C=m2 , sE ¼ 11:7 1012 m2 =N. The piezoelectric constants are identical to the values corresponding to BaTiO3 ceramics as there are no piezoelectric properties reported for tetragonal PZT ceramics in the literature. Note that the piezoelectric constants have only a minor effect on the fractions of 90° switching. The switchable strain in Fig. 9 is about 0.017% and in the case of complete 90° switching it is 0.882% according to Eq. (27). Therefore, the fraction of achievable 90° switching in tetragonal PZT-45/55 ceramics is about 1.9% under electric loading up to 3 MV/m. Note that the simulation results agree qualitatively with the experimental results in [43]. It is clear from the above numerical simulations that the FE formulation of the present constrained domain switching model, which is based on a multi-domain grain assumption and estimates the fraction of 90° switching via a dichotomy approximation, is free from domain switching related numerical instability problems that may be encountered in conventional FE models for ferroelectrics. The present model can predict the fraction of 90° switching in polycrystalline ferroelectrics without using fitting parameters. The current 2D FE model can only be used for tetragonal ferroelectrics, but it can be extended to the 3D cases covering both tetragonal and rhombohedral ferroelectrics. However, the dichotomy approximation in the 2D case cannot be extended in a straightforward
F.X. Li, R.K.N.D. Rajapakse / Computational Materials Science 44 (2008) 322–329
0.08
0.08
Switching strain Total strain
Strain (%)
0.06
0.06
0.04
0.04
0.02
0.02
0.00
0.00
-0.02
-0.02
329
nal stress caused by non-180° switching is naturally considered by the FE model. The fraction of non-180° switching in each grain is obtained through a dichotomy approximation. The switching instability problem is thus well resolved in the proposed constrained domain switching model. Numerical simulations show that during electrical poling, 180° switching is usually complete while 90° switching is only about 2% for both BaTiO3 and tetragonal PZT ceramics, which agrees quite well with the existing experimental results. The FE model can be extended to the 3D cases to cover both tetragonal and rhombohedral ferroelectrics by using further assumptions. Acknowledgement
-3
-2
-1
0
1
2
3
This work was supported by a grant from Natural Science and Engineering Research Council of Canada (NSERC).
E3 (MV/m)
References Fig. 9. Computed strain-electric field response of tetragonal PZT-45/55 ceramics.
z 1
2 5
y 3 4
6
x Fig. 10. Six types of domains in 3D tetragonal ferroelectric crystals.
manner to the 3D case because two or more types of non-180° switching may occur simultaneously. As shown in Fig. 10, when an electric field is applied close to z axis of the crystallite coordinates in tetragonal ferroelectrics, there may exists two types of 90° switching (say from domain 6 and domain 4 to domain 1) simultaneously. But if it is further assumed that domain 4 and domain 6 always switch together at equal fractions to domain 1, then the dichotomy approximation can still be applied to obtain the fraction of non-180° switching. Obviously, the model can also be extended to rhombohedral materials if such an assumption is made. 6. Conclusions A nonlinear finite element (FE) formulation for ferroelectric materials is briefly reviewed and it is illustrated that for FE models based on a phenomenological material model with the electric potential as a nodal variable, the Newton–Raphson iteration is usually unstable. While for FE models based on a micromechanical domain switching criterion, the large depolarization field/stress caused by domain switching makes the domains switch back to its original state, thus switching instability arises. To overcome these domain switching related instabilities and to predict the experimentally observed incomplete non-180° switching in ferroelectric ceramics, a constrained domain switching model is proposed and implemented in a FE code. A multi-domain grain assumption is used and charge screening effect is also taken into account in the present model. The depolarization field by polarization switching is thus compensated by free charges and the inter-
[1] Y. Xu, Ferroelectric Materials and their Applications, Amsterdam Press, NorthHolland, 1991. [2] B. Jaffe, W.R. Cook, H. Jaffe, Piezoelectric Ceramics, Academic Press, London and New York, 1971. [3] H. Allik, T.J.R. Hughes, Int. J. Num. Meth. Eng. 2 (1970) 151–157. [4] C.L. Hom, N. Shankar, Int. J. Solids Struct. 33 (1996) 1757–1779. [5] X. Gong, Z. Suo, J. Mech. Phys. Solids 44 (1996) 751–769. [6] K. Ghandi, N.W. Hagood, Smart structures and materials: mathematics and control in smart structures, in: V.V. Varadan (Ed.), Proceedings of the SPIE 2715 (1996) 121–140. [7] K. Ghandi, N.W. Hagood, Smart structures and materials: mathematics and control in smart structures, in: V.V. Varadan, J. Chandra (Eds.), Proceedings of the SPIE 3039 (1997) 97–112. [8] M. Kamlah, Ü. Böhle, Int. J. Solids. Struct. 38 (2001) 605–633. [9] M. Kamlah, A.C. Liskowsky, R.M. McMeeking, H. Balke, Int. J. Solids Struct. 42 (2005) 2949–2964. [10] J.E. Huber, N.A. Fleck, J. Mech. Phys. Solids 49 (2001) 785–811. [11] W. Chen, C.S. Lynch, J. Eng. Mater. Tech. 123 (2001) 169–175. [12] T. Felt, D. Munz, G. Thun, J. Am. Ceram. Soc. 86 (2003) 1427–1429. [13] F.X. Li, D.N. Fang, J. Intell. Mater. Syst. Struct. 16 (2005) 583–588. [14] S.C. Hwang, C.S. Lynch, R.M. McMeeking, Acta Metall. Mater. 43 (1995) 2073– 2084. [15] S.C. Hwang, R.M. McMeeking, Ferroelectrics 211 (1998) 177–194. [16] S.C. Hwang, R.M. McMeeking, Int. J. Solids Struct. 36 (1999) 1541–1556. [17] W. Chen, C.S. Lynch, Eng. Fract. Mech. 64 (1999) 539–562. [18] T. Steinkopff, J. Eur. Ceram. Soc. 19 (1999) 1247–1249. [19] R.M. McMeeking, S.C. Hwang, Ferroelectrics 200 (1997) 151–173. [20] F.X. Li, D.N. Fang, Mech. Mater. 36 (2004) 959–973. [21] D. Berlincourt, H.H.A. Krueger, J. Appl. Phys. 30 (1959) 1804–1810. [22] E.C. Subbarao, M.C. McQuarrie, W.R. Buessem, J. Appl. Phys. 28 (1957) 1194– 1200. [23] O. Axelsson, Iterative Solution Methods, Cambridge University Press, Cambridge, 1994. [24] M.A. Crisfield, Non-linear Finite Element Analysis of Solids and Structures, vol. 1, Wiley, Chichester, 1991. [25] C.M. Landis, Int. J. Num. Meth. Eng. 55 (2002) 613–628. [26] M. Kamlah, Z. Wang, Comput. Mater. Sci. 28 (2003) 409–418. [27] C.M. Landis, J. Appl. Mech.-T ASME 70 (2003) 470–478. [28] C.M. Landis, Current Opin. Solid State Mater. Sci. 8 (2004) 59–69. [29] V. Mehling, Ch. Tsakmakis, D. Gross, J. Mech. Phys. Solid 55 (2007) 2106–2141. [30] S.C. Hwang, J.E. Huber, R.M. McMeeking, N.A. Fleck, J. Appl. Phys. 83 (1998) 1530–1540. [31] H. Kessler, H. Balke, J. Mech. Phys. Solids. 49 (2001) 953–978. [32] H. Kessler Jr., E.R. Fuller, H. Balke, Smart structures and materials: behavior and mechanics of active materials, in: C.S. Lynch (Ed.), Proceedings of the SPIE 3992 (2000) 234–244. [33] X.W. Zhang, C. Lei, K.P. Chen, J. Am. Ceram. Soc. 88 (2005) 335–338. [34] M.E. Lines, A.M. Glass, Principles and Applications of Ferroelectrics and Related Materials, Clarendon Press, Oxford, 1977. [35] D.C. Lupascu, Fatigue in Ferroelectric Ceramics and Related Issues, SpringerVerlag Berlin Heidelberg, Berlin-New York, 2004. [36] G.I. Taylor, J. Inst. Metals 62 (1938) 307–324. [37] K. Bhattacharya, R.V. Kohn, Acta Mater. 44 (1996) 529–542. [38] J.Y. Li, R.C. Rogan, E. Üstündag, K. Bhattacharya, Nature Mater. 4 (2005) 776– 781. [39] S. Wan, K.J. Bowman, J. Mater. Res. 16 (2001) 2306–2313. [40] F.X. Li, R.K.N.D. Rajapakse, J. Appl. Phys. 101 (2007) 054110. [41] Q.M. Zhang, H. Wang, N. Kim, L.E. Cross, J. Appl. Phys. 75 (1994) 454–459. [42] W. Lu, D.N. Fang, K.C. Hwang, Acta Mater. 47 (1999) 2913–2926. [43] M.J. Hoffmann, M. Hammer, A. Endriss, D.C. Lupascu, Acta Mater. 49 (2001) 1301–1310.