Multi-grain analysis versus self-consistent estimates of ferroelectric polycrystals

ARTICLE IN PRESS

Journal of the Mechanics and Physics of Solids 55 (2007) 648–665 www.elsevier.com/locate/jmps

Multi-grain analysis versus self-consistent estimates of ferroelectric polycrystals Anja Hauga, John E. Huberb, Patrick R. Oncka, Erik Van der Giessena, a

Materials Science Center, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands b Department of Engineering Science, University of Oxford, Parks Rd., Oxford OX1 3PJ, UK Received 6 February 2006; received in revised form 29 May 2006; accepted 26 June 2006

Abstract Ferroelectrics are polycrystalline materials consisting of intragranular regions with different polarization directions, called domains. The domains can be switched into different states by the application of an electric field or mechanical stress. We study the influence of grain-to-grain interactions on the overall and local switching behavior. The behavior inside each grain is represented by the micromechanics model of Huber et al. [1999. A constitutive model for ferroelectric polycrystals. J. Mech. Phys. Solids 47 (8), 1663–1697]. The predictions of a self-consistent model of the polycrystal response are compared with those of a multi-grain model in which grains are represented individually. In one flavor of the multi-grain model, each grain is represented by a single finite element, while in the other the fields inside each grain are captured in more detail through a fine discretization. Different electrical and mechanical loading situations are investigated. It is found that the overall response is only mildly dependent on the accuracy with which grain-to-grain interactions are captured, while the distribution of grain-average stresses is quite sensitive to the resolution of the intragranular fields. r 2006 Elsevier Ltd. All rights reserved. Keywords: Ferroelectricity; Self-consistent model; Polycrystals; Multi-grain model

Corresponding author.

E-mail address: [email protected] (E. Van der Giessen). 0022-5096/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.jmps.2006.06.009

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

649

1. Introduction A range of device applications, such as actuation and sensing, are based on the piezoelectric and inverse piezoelectric effects. Ferroelectric ceramics such as lead zirconate titanate (PZT) are employed in such devices because of their high piezoelectric coupling coefficient. These polycrystalline ceramics are not piezoelectric after sintering. This is because the remanent polarization is zero after processing and the piezoelectric coefficients depend linearly on the remanent polarization (Jaffe et al., 1971). However, the application of stress and electric field can change the state of the remanent polarization and thus the device performance. In order to calculate the device performance a nonlinear constitutive law needs to be employed, for which a number of possibilities have been proposed in the literature. They differ in the level of representation of the material microstructure. The phenomenological models of Landis (2002), McMeeking and Landis (2002), Cocks and McMeeking (1999), Kamlah and Tsakmakis (1999) are completely macroscopic and do not explicitly account for any microstructure. On the other hand, the phase-field model of Wang et al. (2004) represents the motion of the domain walls in detail but does not yet allow to compute the performance of a polycrystalline ferroelectric. The model proposed by Huber et al. (1999) for a single crystal sits in between these two classes of models, and involves a ‘smeared-out’ representation of the domain structure. From this single-crystal model the response of a polycrystal can be computed in various ways. As a refinement of a Taylor-like model, in which interactions between crystallites are ignored, Huber et al. (1999) proposed a self-consistent procedure. In this approach, the interaction between grains is approximated by considering all grains to be embedded as ellipsoidal inclusions in the macroscopic, to be determined, matrix. By analogy with the self-consistent models for polycrystal plasticity, Huber et al. (1999) expected that a self-consistent model would give a rather accurate description of the overall response of the polycrystal. It is unclear however how accurate the self-consistent estimates of the fields inside individual grains are, while, for example, the stress fields inside grains control the susceptibility to damage development and fracture. The objective of the present paper is to study the accuracy of the self-consistent average by comparison with a two-dimensional multi-grain model in which grains inside a polycrystalline aggregate are represented individually. In addition to the predicted overall response, we pay special attention to the way in which grain–grain interactions affect the distributions of strain, remanent strain, and stresses inside the grains. The multi-grain model is analyzed using finite elements, where we adopt two levels of discretization. In the coarser approach, a single element is used to represent a grain, while we also discuss simulations where hexagonal grains are meshed with a relatively fine grid. A similar study, recently presented by Kamlah et al. (2005), used a single element per grain and found that the influence of intergranular constraint on the macroscopic material response is significant close to saturation and during reversal of loading. The emphasis of the present work is on identifying the intragranular stresses and their potential consequences for material reliability; we therefore include a detailed analysis of the distribution of field quantities, in addition to the macroscopic response.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

650

2. Single crystal model The model we use for a single crystal is described in detail by Huber et al. (1999). Each grain is considered to comprise a number of different domains, corresponding to the possible directions of spontaneous polarization. In the case of a tetragonal system (e.g. tetragonal perovskite, BaTiO3 ) six domains (M ¼ 6), or variants, are possible. Assuming that the ratio of grain size to the size of the domains is large, the model does not describe the domains individually but rather adopts a ‘smeared-out’ representation. At each material point in the crystal, the relative size of a domain, I, is characterized by its volume fraction cI , with values between zero and one. It is further assumed that in each domain the strain eij and the electrical displacement Di are composed of a linear part, proportional to the stress skl and electric field E k , and an independent remanent part. The total response is calculated by volume averaging over all domains, both for the remanent (superscript R) and the linear part (superscript L), and is given by eij ¼ eLij þ eR ij ¼

M X

EðIÞ cI ½sijkl skl þ d ðIÞ kij E k  þ

I¼1

Di ¼ DLi þ DR i ¼

M X

cI eRI ij ,

(1)

I¼1

M X

sðIÞ cI ½d ðIÞ ikl skl þ kik E k  þ

I¼1

M X

cI PRI i ,

(2)

I¼1

ðIÞ where, for each domain I, sEðIÞ ijkl is the elastic compliance tensor at constant electric field, d kij

the dielectric tensor, and ksðIÞ ik the dielectric permittivity tensor at constant stress. Note that the latter material parameters are fixed to specific domains, so that the overall elastic, dielectric and piezoelectric properties depend explicitly on the domain volume fractions— which change due to switching. Domains switch through switching systems. This leads to N ¼ MðM  1Þ ¼ 30 switching systems, with separate switching systems for forward and reverse switching. Thus, not all switching systems can be active: denoting the volume fraction of switching a system a ¼ 1; . . . ; N by f a , the rate of change f_ ¼ 0 in the nonactive systems. The rate of a change c_I of domain I’s volume fraction caused by the rate of change f_ is given by c_I ¼

N X

a

AIa f_ ;

I ¼ 1; . . . ; M

(3)

a¼1

in terms of a matrix AIa populated with 0’s and 1’s. The switching increment produces an increment of remanent strain and remanent polarization. The remanent strain of a switching system is described by simple shear in a direction si on a plane of unit normal ni through the Schmid tensor, mij ¼ 12 ðsi nj þ sj ni Þ, as in crystal plasticity (see e.g. Asaro, 1983). The remanent polarization changes in the direction of shear, si . Thus, the rates of remanent R strain, e_R , and polarization, P_ , are obtained by summing the contributions over all ij

i

switching systems: X a f_ maij ga ; e_R ij ¼ a

R P_ i ¼

X

a f_ sai Pa .

(4)

a

Here, ga and Pa are the characteristic strain and polarization of switching system a, respectively.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

651

The rates of changes of these remanent parts are governed by the thermodynamic driving force Ga for switching, appearing in the dissipative work expression X a G a f_ . (5) w_ D ¼ a

From the observation that the latter is the total work rate minus the rate of change of the recoverable stored work, G a can be expressed as (Huber et al., 1999) a

G a ¼ ðmij ga þ 12e~aij Þsij þ ðsai Pa þ 12D~ i ÞE i , where e~ aij ¼

X

ðIÞ AIa ½sEðIÞ ijkl skl þ d kij E k ;

a D~ i ¼

I

(6) X

sðIÞ AIa ½d ðIÞ ikl skl þ kik E k 

(7)

I

are the linear strains and electric displacement, respectively, associated with switching system a. 3. Self-consistent model for polycrystals For use in the self-consistent description, the single crystal model developed in Section 2 a is completed by relating the volume fraction transfer rates f_ to the thermodynamic driving a a force for transformation G . A rate-independent formulation is obtained by setting f_ ¼ 0 a a a unless G ¼ G c , where G c is a threshold value of the driving force required to initiate switching. Linear hardening is incorporated through X b a H ab f_ . (8) G_ c ¼ b

In the present work, the hardening matrix is taken to be H ab ¼ dab G c =100 to approximate ;tan tan nonhardening behavior. Tangent moduli (cD;tan ijkl , hkij , bik ) for the single crystal are defined by _ _kl  htan s_ ij ¼ cD;tan kij Dk , ijkl  E_ i ¼

_kl htan ikl 

þ

_ b;tan ik Dk .

ð9Þ ð10Þ

Following Huber et al. (1999), the tangent moduli for a single crystal are derived in Appendix A. The self-consistent average of the tangent moduli over a set of grains with differing crystallographic orientations provides an estimate of the overall tangent moduli that relate the macroscopic increments of stress and electric field ðs_¯ ij ; E_¯ i Þ to those of strain and electric _¯ Þ: displacement ð_¯ ij ; D i _¯ , _¯kl  h0kij D s_¯ ij ¼ cD;0 k ijkl  ;0 0 _ _¯ . E¯ i ¼ hikl _¯ kl þ bik D k

ð11Þ ð12Þ

Here, an overbar indicates the volume average of a quantity over the entire polycrystal. The key feature of the self-consistent average relevant to the present study is that each grain experiences a constraint imposed upon it by compatibility with the surrounding material. Thus, the macroscopic behavior of the model captures features that result from the interaction between grains. However, in the self-consistent scheme, each grain is

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

652

constrained only by a matrix of material with the overall tangent moduli of the polycrystal. Thus, the details of grain-to-grain interactions due to particular arrangements of grains are lost. Self-consistent calculations are here compared with a multi-grain finite element analysis where the detailed interactions are captured explicitly. In order to compare the self-consistent estimate of the polycrystal response with the multi-grain model, calculations of material response were carried out using identical material parameters to those used in the multi-grain model, which are taken from Table 1 of Huber et al. (1999). Comparable conditions were produced by imposing ¯ 2 ¼ 0), plane strain and plane electric displacement conditions (¯12 ¼ ¯ 22 ¼ ¯ 32 ¼ 0, D and by choosing the crystallographic orientations such that each grain has its local 2-axis aligned with the global 2-axis (the out-of-plane direction). The stress, strain, electric field, electric displacement and corresponding remanent quantities were calculated in each grain at each load step to provide a detailed comparison with the multi-grain simulations.

4. Multi-grain model The self-consistent estimate of the polycrystalline response is compared with that predicted by a model in which a polycrystal is represented by an aggregate of discrete grains. In two dimensions here, two types of aggregates are considered: one containing many triangular grains (see Fig. 1a), the other comprising a smaller set of hexagonal grains, as illustrated in Fig. 1b. In the first type, each grain is represented by single (linear triangular) finite element, while in the hexagonal-grain model each grain is meshed with a relatively fine grid of finite elements. Contrary to the solution method in the self-consistent model, where switching occurs when G a reaches a critical value G ac , we here adopt a power-law relation between the driving force and the incremental volume fraction of the material for numerical

w

w grain R 3

ωR 2

1 h

h

ti=0 Q=0

(a)

t1=0

ti=0

ti=0

ti=0

Q=0

Q=0

Q=0

(b)

t1=0

Fig. 1. Definition of multi-grain model with boundary conditions: (a) one-element triangular grain model; (b) hexagonal grain model with a refined finite element mesh shown in the inset.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

convenience (Huber and Fleck, 2001). Thus, 8     < f_ Ga m cI 1=k for G a 40; a 0 Gac c0 _ f ¼ :0 for G a o0;

653

(13)

where f_0 , m and k are model constants, and I specifies the domain in which transformation a takes place. Gac now is the value of G a at which there is a rapid increase in transformation rate for mb1 and by taking G ac ¼ Gc there is no hardening, as in the self-consistent calculations. The last factor in Eq. (13) has been added to enforce the limits of the switching behavior: when domain type I is depleted, i.e. cI ¼ 0, it cannot contribute to further switching. c0 is the volume fraction of each variant in the unpoled state. Because of the high nonlinearity introduced by (13), use of the incremental form of the expressions (1) and (2) with (4) would lead to a rather poor numerical stability. Therefore a forward gradient method is used following Pierce et al. (1984), the details of which can be found in Appendix B. The finite element formulation adopts material displacements and electric potential as nodal variables (Allik and Hughes, 1970). Therefore, the increments of stress and electric displacement need to be expressed in terms of the increments of strain and electric field. Thus, the rate constitutive equations used in the finite element formulation read _ _ ij , _ s_ ij ¼ cEijkltan e_kl  etan kij ðE k þ E k Þ þ s D_ ¼ etan e_  ktan ðE_ þ E_ Þ,

ð14Þ

ð15Þ k k ikl kl ik _ where E i and s_ ij are defined in Appendix B, Eq. (36). The governing equations are incorporated through the virtual work statement Z Z Z _ Dt ðs_ ij deij  D_ i dE i Þ dV ¼ Dt ðt_i dui  QdfÞ dS  ðsij deij  Di dE i Þ dV V S V  Z  ðti dui  QdfÞ dS 8dui ; df ð16Þ i

S

with virtual strains and virtual electric field defined by deij ¼ 12ðdui;j þ duj;i Þ;

dE i ¼ df;i ,

respectively. Here ui is the displacement, f the electrical potential and Q the surface charge. The terms in square brackets on the right-hand side is an equilibrium correction, included to prevent the solution from drifting off the equilibrium path. 5. Results Although the problem definition in Fig. 1 assumes plane strain and plane electric ¯ 2 ¼ e¯i2 ¼ 0 (i ¼ 1; 2; 3) for the self-consistent and D2 ¼ ei2 ¼ 0 for the displacement, i.e. D finite element calculations, the crystals considered are three-dimensional crystals, but with the normals to the switching systems pointing in the x2 -direction, normal to the plane of view. The crystallographic axes of grain R ¼ 1; . . . ; n are described by a rotation oR about the x2 -axis which is taken to be specified by oR ¼ 2pR=n for the self-consistent calculation and chosen randomly within 0 to 2p for the multi-grain model.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

654

Results for four different loading situations are presented: uniaxial electric cycling, uniaxial mechanical cycling, uniaxial electric cycling under a compressive stress, and mechanical depolarization. The boundary conditions for the finite element calculation are specified in Fig. 1 and Table 1. The different loads are applied on the top and bottom side via the potential and tractions (and indicated by the subscript ‘app’). The lateral sides are always traction and charge free. It is assumed that the bottom of the sample cannot move in the x3 -direction and is shear–stress free t1 ¼ 0. The loading conditions in the selfconsistent calculation are sapp ¼ s¯ 33 and E app ¼ E¯ 3 . Depending on the type of loading situation, the macroscopic response is reported in terms of the overall electric displacement, Dmac , or the macroscopic strain, emac . In the multi-grain computations, the macroscopic electric displacement is calculated from the charges on the top divided with the area of the top surface, while emac is computed from the average displacement u3 along the top of the sample divided with the cell height h. In ¯ 3 and emac ¼ e¯33 . the self-consistent calculation, these are simply given by Dmac ¼ D One thousand eight hundred grains are used in the self-consistent calculation and in the one-element grain calculation, and 187 in the hexagonal-grain mesh (Fig. 1). All hexagonal grains in the multi-grain model of Fig. 1b are meshed with 72 triangular elements. The values of the material constants in the crystallographic axes and the modeling parameters for the switching systems can be found in Table 1 of Huber et al. (1999). In each grain R, ðIÞ sðIÞ the sEðIÞ ijkl , d ijk , kik , mij and si are obtained from the orientation oR . In the initially unpoled state, all six domain volume fractions are taken equal to cI ¼ c0 ¼ 16. All results are normalized by D0 , E 0 , g0 or t0 as in Huber et al. (1999). D0 and E 0 are the values at the onset of switching for purely electrical loading; similarly, g0 and t0 mark the start of switching under mechanical loading. 5.1. Electrical cycling Before comparing the results from the different frameworks, we explore the effect of the exponent m in the power-law (13). The self-consistent model is time independent, corresponding to m ! 1 in the power-law (13). Rate independence is approached by using a high value of m, but this requires very small timesteps. To asses the influence of m, calculations are done for a cyclic electric field of 4E 0 of a grain consisting of one element (so that stress, strain and electric field are uniform, as in the self-consistent calculation). In Fig. 2 the macroscopic response, hysteresis and butterfly loop, are shown for m ¼ 5; 20 and 100. The other parameters in (13) are k ¼ 0:6, f_0 ¼ 0:4 s1 , while the loading rate is E_ ¼ 5:56E 0 f_0 . The rate-dependent calculation predicts the typical features of a hysteresis loop as well as of a butterfly loop. At initial loading only a dielectric response occurs. Table 1 Boundary conditions that differentiate between the different kinds of loadings considered in this paper

Top Bottom

Electrical cycling

Mechanical cycling

Depoling

Electrical cycling under stress

_ ¼ E_ app h f t3 ¼ 0 _ ¼0 f

Q¼0 t_3 ¼ s_ app Q¼0

f¼0 t_3 ¼ s_ app f_ ¼ 0

f_ ¼ E_ app h t_3 ¼ s_ app f_ ¼ 0

The boundary conditions that are common to all cases are shown in Fig. 1.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

εmac/γ0

Dmac/D0

1.5

6 FEM

m 5 20 100

655

4 1 2 0

-4

-2

0 -2

2 Eapplied/E0

4 0. 5

-4 -6

0 -4

-3

-2

-1

0

1

2

3

4

Eapplied /E0 Fig. 2. Hysteresis and butterfly loop for m ¼ 5; 20; 100 in (13) for a single-element grain calculation under electrical cycling with 4E 0 .

At the onset of nonlinearity there is an increase in electric displacement due to the appearance of remanent polarization; similarly, there is an increase in strain due to remanent strain and due to the nonzero piezoelectric tensor. Subsequently, the response becomes linear again. Such a linear response may return once the volume fractions are depleted and no further increase in remanent components can occur. However, this is not the case here; instead, some mild switching still occurs in what looks like the post-poling regime, and is due to the stress that develops in the through-thickness direction, s22 , as a consequence of the plane strain condition (Haug et al., 2006). When the electric field is reversed, the response is again linear at first, until reverse switching starts. First, the material switches mainly due to stresses in the through-thickness direction, and this emerges as a linear reduction in electric displacement and strain. At the onset of the nonlinearity due to the applied electric field, the electric displacement decreases tremendously, whereas the strain only increases due to the change in piezoelectric properties. With further decreasing electric field, once switching in the reverse direction is exhausted, the electric displacement decreases and the strain increases both linearly. Reversing the electric field once more gives the same behavior. We observe that the value of m influences the behavior at the onset of nonlinear behavior. The higher the value of m, the steeper the slope in the hysteresis loop during nonlinear behavior and the sharper and higher is the tip of the butterfly loop. The results for m values of 20 and 100 coincide when no switching occurs, but differ slightly during switching. For m ¼ 5 the response during switching is so different that even the strain in the linear regimes is affected by this. In the remainder of this paper, we use m ¼ 20 as a compromise between computing time and near time-independent behavior. We proceed by comparing the results of the self-consistent model with the multi-grain model. Since there is a spatial distribution of differently oriented grains in the multi-grain model, 24 realizations of triangular grains have been analyzed, one example of which is shown in Fig. 3. According to Fig. 3a, all calculations agree very well in the electrical response taking into consideration the difference in rate dependence between the two classes of computations. In the mechanical response, Fig. 3b, the agreement between the

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

656

Dmac/D0

εmac/γ0

6

P

1.8 1.6 1.4

4

1.2

2

1 0 -4

-3

-2

-1

0

1

-2

2 3 Eapp/E0

P

0.8

4

0.6 0.4

-4 SC FE M FE M

-6

0.2 0 -4

(a)

-2

0 Eapp/E0

(b)

2

4

Fig. 3. Comparison of macroscopic response predicted by the self-consistent model (SCM) and the multi-grain model (FEM) for one particular realization under electrical cycling in terms of (a) electrical response (hysteresis loop) and (b) mechanical response (butterfly loop).

R ε33 /γ0

R ε33 /γ0

2

2 1.5

1.5

1

1

0.5

0.5 8

6

4

2

(a) frequency [%] SCM

0

0

(b)

90 180 270 grain Orientation [deg]

FEM

(24 realizations)

360

0

2

4

6

8

(c) frequency [%] FEM

(1 realization)

Fig. 4. Remanent strain eR 33 =g0 versus grain orientation (b) and its frequency distribution according to SCM (a) and multi-grain model (FEM) (c) at E app ¼ 4E 0 . For an explanation of the plotting procedure, see the text. The two blue curves in (b) indicate the spread among the 24 single-element grain computations analyzed.

hexagonal grains calculation and the self-consistent calculation is very good except for the lower tips. The butterfly loop of the triangular grain calculation has the same size but is slightly shifted to lower strains. In order to gain more insight in the differences between the various models, the R microscopic response is investigated. For this purpose, we analyze the values of eR 33 , P3 , s33 and e33 at the peak applied electric field of 4E 0 . In Figs. 4a–7a, histograms, histograms of the frequency distribution of the different values inside each grain of the SCM are plotted. The same is done in Figs. 4c–7c for one realization of hexagonal grains. In this case, there are several data points per orientation corresponding with the number of integration stations per grain. The blue curve in these plots gives the distribution averaged over all grains in the 24 one-element grain calculation. For more information, the variation of each parameter with respect to the grain orientation is plotted in Figs. 4b–7b where each grain is represented by one symbol. In the case of the self-consistent model, adjacent symbols

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665 P3R /D0

P3R /D0 5

5

4.5

4.5

4

4

3.5

3.5

3

3 12 10 8 6 4 2 (a) frequency [%]

0

657

0

(b)

90 180 270 grain Orientation [deg]

SCM

FEM

360

(c)

0 2 4 6 8 10 12 frequency [%]

FEM

(24 realizations)

(1 realization)

Fig. 5. Remanent polarization PR 3 =D0 versus grain orientation (b) and its frequency distribution of SCM (a) and multi-grain model (FEM) (c) at E app ¼ 4E 0 .

σ33 /τ0

σ33 /τ0 5 4 3 2 1 0 -1 -2 -3

5 4 3 2 1 0 -1 -2 -3 10 8 6 4 2 0 (a) frequency [%]

0

90

180

270

grain Orientation [deg]

(b)

SCM

FEM

0

360 0 2 4 6 8 10 (c) frequency [%]

(24 realizations)

FEM

(1 realization)

Fig. 6. Intra-grain stress in direction of applied field, s33 =t0 , versus grain orientation (b) at and its frequency distribution of SCM (a) and multi-grain model (FEM) (c) at E app ¼ 4E 0 .

ε33/γ0

ε33/γ0

4

SCM

FEM

4

(24 realizations)

3.5 3

2.5

2

2

1.5

1.5

1

1

0.5

(a)

6

4

frequency [%]

2

0

(1 realization)

3

2.5

8

FEM

3.5

0.5 0

(b)

90

180

grain Orientation [deg]

270

360

0

(c)

2

4

6

8

frequency [%]

Fig. 7. Total strain e33 =g0 versus grain orientation (b) and its frequency distribution of SCM (a) and multi-grain model (FEM) (c) at E app ¼ 4E 0 . The dash–dotted horizontal lines in (a) and (c) indicate the median of e33 =g0 .

string together to form a seemingly smooth curve. In contrast, values from the integration stations of every hexagonal grain exhibit a spread due to the nonuniformity over the grain, thus resulting in a wide distribution for all grains. To show the trend of the one-element

ARTICLE IN PRESS 658

A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

grain calculation, the results of 24 realizations are presented in terms of 80 equal intervals of grain orientation between 0 and 2p. For each interval, the mean value and the standard deviation are calculated; in the figures, two curves are plotted: the mean values plus and minus the standard deviation. In all cases, the results show an orientational dependence with a periodicity of 90 . This originates from the fact that the switching system itself has a rotational symmetry of 90 in the plane. Fig. 4 presents the orientation dependence of the remanent strains. It is noted that only 90 switches contribute to the remanent strain; 180 switching is strain free. The contributions of grains whose orientation aligns with the global coordinate system (oR ¼ 0 ) are highest and are lowest for oR ¼ 45 . The frequency plot, Fig. 4a, shows a bimodal distribution. This is a direct consequence of the slope of the sinusoidal curve in Fig. 4b: regions near the peaks and valleys with a very low slope give a relatively large contribution to the distribution function. The single-element grain computations are seen in Fig. 4b to show the same trend. The scatter, indicated by the two blue lines, is observed to be dependent on the grain orientation; around the valleys, where there is very little switching, the standard deviation is smaller than around the peaks, yielding values that are lower than the SCM data. The frequency plot, Fig. 4c, also reveals a bimodal distribution but less pronounced than for the SCM data. The shape of the histogram is also influenced by the scatter behavior; a higher standard deviation broadens and reduces the peak. This explains the broadening of the peak at the upper range of the remanent strain. The hexagonal-grain calculation shows the same tendency in the scatter plot and frequency histogram, but the large-strain peak in the frequency histogram is shifted to higher values.   The remanent polarization PR 3 , analyzed in Fig. 5, is affected by both 90 and 180 switches. In the SCM calculation, three regions (Fig. 5b, horizontal black lines, seen also in the frequency plot), appear inside the sinusoidal-like dependence. The outside regions give rise to the bottom and top peaks in the frequency histogram Fig. 5a as in Fig. 4a. But, now a third weak peak appears in the middle which is caused by the different slope of the curve in between the horizontal lines in the scatter plot. The ensemble-averaged remanent polarization of the one-element grain results follow roughly the same trends but less structure can be obtained due to the scatter of the data. One would again expect a bimodal distribution, but since the standard deviation around the peaks in the orientation dependence is so large, the histogram in Fig. 4c is relatively smooth with just a minor peak at the lower shelf of polarizations. The results for the hexagonal grain are similar. The scatter plot of the self-consistent grain stresses in the polarization direction, s33 , in Fig. 6b as a function of grain orientation reflects a sinusoidal dependence. However, there is a phase shift by 45 compared to the remanent strain and remanent polarization curves: this means that a tensile stress corresponds to a small remanent strain and a compressive stress to a large remanent strain. The single-grain results follow a similar dependence but due to the large standard deviation, the frequency plot shows a normal distribution, Fig. 6c, and not a bimodal distribution as for the self-consistent model, Fig. 6a. The hexagonal grains behave the same way as the triangular grains. In the total strain histograms, Fig. 7a–c, also the median values of self-consistent and hexagonal-grain simulations, plotted in Fig. 3, are shown as dashed–dotted lines. According to Eq. (1), the total strain has contributions from the local stresses due to the compliance, from the electric field due to the piezoelectricity and from the remanent strain. The electric field in the applied direction in the grains fluctuates slightly about the

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

659

σ33/τ0 3.00 2.00 1.00 -1.00 -2.00 -3.00

(a)

(b)

Fig. 8. Distribution of s33 =t0 at E app ¼ 4E 0 in a particular realization of orientations according to the multi-grain models: (a) multiple hexagonal grains and (b) single-element grains.

applied electric field value; thus in comparison with the remanent strain curve, the piezoelectric effect shifts up the total strain. The stresses (Fig. 6) generate elastic strains that decrease the amplitude of the remanent strain curve (Fig. 4). Finally, to gain further insight in this, we analyze the stress distributions in multi-grain aggregates shown in Fig. 8. Note that, on average, s33 is zero, in correspondence with Fig. 6. In the hexagonal grain model, Fig. 8a, high compressive and tensile stresses are found mostly near the grain boundaries, whereas in the single-element calculation the whole grains are affected by the high stresses, Fig. 8b. Thus, in the hexagonal grains mainly the material around the grain boundaries has to accommodate with the surrounding grains thereby generating locally higher stresses than in the one-element grain calculation (see Fig. 6). In summary, the remanent strain and polarization depend predominantly on grain orientation through the orientation factors maij and sai , respectively, in the micromechanics law (4). The enhanced grain-to-grain interactions in the multi-grain models smear out the distributions as compared to the self-consistent results (compare Fig. 4a with c, and Fig. 5a with c). Both the multi-grain models and the SCM predict that compressive stresses develop when the remanent strain in the poling direction is large, and that tensile stresses are induced when the remanent strain is small. This translates into macroscopic predictions that agree well, but at the microscopic grain level the self-consistent predictions deviate substantially from those of the multi-grain model. In particular, the bi-modal distribution functions in Figs. 7a and 6a are unrealistic in comparison with the more Gaussian distributions according to the multi-grain model, Figs. 7c and 6c. This implies that the maximum stresses and strains that can develop inside grains in order to maintain compatibility are severely underestimated by the self-consistent method. 5.2. Mechanical cycling Next, we consider mechanical cycling (cf. Table 1) by loading between 8t0 at a rate of s_ app ¼ 0:52t0 f_0 . The material response, shown in Fig. 9, is analogous to electrical cycling.

ARTICLE IN PRESS 660

A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

8

σapp/τ0

6 4 2 0 -4

-3

-2

-1

0

1 -2 -2.8 -4 -6

2

3

4

εmac/γ0 SC FEM FEM

-8 Fig. 9. Macroscopic strain in x3 -direction due to stress-driven mechanical cycling.

First, the compressive response is linear until switching starts. The subsequent nonlinear response ceases when volume fractions deplete, after which the response is linear again. After reducing the applied stress the behavior repeats itself. It is important to note that this is a purely mechanical response, without a change in remanent polarization and therefore involves no piezoelectricity. Mechanical loading of an unpoled material gives rise to switching where the volume fractions change in such a way that no remanent polarization arises but only remanent strain. The agreement between the different model predictions is very good, with just minor deviations for the single-element grain model. During the nonlinear process the triangular grains behave slightly differently, but the strains at 8sapp are the same for all models. 5.3. Mechanical depoling Until now either a purely electrical or a purely mechanical loading situation has been investigated. Fig. 10 shows the depolarization response: the material is poled first by electrical loading and then a mechanical compressive stress is applied (s_ app o0) and removed. Since the application of a compressive stress in the x3 -direction does not have a preferred switching direction but rather a switching plane—the ðx1 ; x3 Þ-plane—the net polarization gets destroyed. For poling, an electric field of 4E 0 is applied and then reduced to zero with a rate of 5:56E 0 f_0 ; this corresponds to the point P in Fig. 3. Subsequently, a mechanical field of 0:52t0 f_0 up to 8t0 is applied and reduced to zero, as in Fig. 9. The application of the mechanical stress yields a nonlinear response until the volume fractions are depleted or the stresses prevent further switching, after which the response is linear again. As can be seen in Fig. 10, the macroscopic values of strain and electric displacement predicted by the various models after poling differ somewhat (point P, also in Fig. 3). During depoling from this state, the agreement in mechanical response is very good but the electrical responses remain slightly different. However, the total electric displacement difference during depoling is the same for all calculations.

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

εmac/γ0 -3

-2

-1

Dmac/D0 0

1

0

-1

P

2

3

4 P

-1 -2

-3

-3

σapp/τ0

-2

-4

1

0

0

(a)

661

-4

-5

-5

-6

-6

-7

-7

-8

-8

SC FEM FEM

(b)

Fig. 10. Macroscopic strain (a) and macroscopic electric displacement (b) during mechanical depoling from the state denoted by P in Fig. 3.

During unloading, the response is linear and the same for all three calculations. The finite slope in Fig. 10b indicates that the material is still piezoelectric, i.e. not fully depoled. The constant and equal slopes in Figs. 10a and b, during unloading, indicates that the piezoelectric tensor is on average the same in all calculations. 5.4. Electrical cycling under mechanical stress Finally, another mixed loading condition is investigated: electrical cycling under a compressive load, Table 1. Therefore a compressive stress of 2:8t0 is applied with a loading rate of 0:52t0 f_0 and then the electric field is cycled with a rate of 5:56E 0 f_0 . The applied stress already gives rise to extensive switching, as seen previously in Fig. 11, but in the ðx1 ; x3 Þ-plane and therefore without a net polarization. Compressive stresses in throughthickness direction do develop. The effect of the compressive stress can be seen by comparing with Fig. 3. The butterfly loop is shifted to lower values due to the applied compressive stress. The loop also loses symmetry about the strain axis. Quite interestingly, the sign of the strain difference calculated at an applied electric field of þ4E 0 and 4E 0 is opposite for FEM and SCM. Switching starts with the application of the electric field in a complex stress state. The applied stress makes it harder for the material to switch. The total strain difference is about 90% of the uncompressed loop. The tip of the butterfly loop broadens and the values are smaller than the strain at the beginning of electrical cycling indicating that less 180 switches occur. Also due to this, the steep slope around E 0 in Fig. 3 in the hysteresis loop flattens out and thus no sharp distinction between linear and nonlinear response is visible. The macroscopic displacement difference from the top to the bottom wing of the butterfly loop at an applied electric field of 4E 0 reaches only about 80% of the value when no compressive stress is applied. Again the agreement in the hysteresis loop of all the

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

662

εmac/γ0

Dmac/D0

-0.4

6

-0.6

4

-0.8 2 -1 0 -4

-3

-2

-1

0 -2

1

2

3

-1.2

4

-1.4

Eapp/E0

-4 -6

(a)

-1.6 SC FEM FEM

-1.8 -4 (b)

-3

-2

-1

0

1

2

3

4

Eapp/E0

Fig. 11. Hysteresis (a) and butterfly loop (b) during electrical cycling under a compressive stress of 2:8t0 .

modeling results is very good. In the butterfly loop, the top of the wing at E app ¼ 4E 0 varies between the SCM and FEM calculation. However, note that the FEM results compare well with each other, even for the tips of the loops. 6. Summary Results of self-consistent calculations and multi-grain finite element calculations, based on the single crystal law of Huber et al. (1999), have been compared. Four different loading situations were investigated: electrical and mechanical cycling, mechanical depoling and electrical cycling under compression. The general conclusion is that the agreement in macroscopic response is excellent for all loading situations. The ‘microscopic’ response at the (sub)grain scale was investigated for poling by an applied electric field of 4E 0 . It has been found that both the remanent properties in the direction of the electric R field, eR 33 , P3 , and mechanical fields, s33 and e33 , show a clear grain orientation dependence. But, in the multi-grain results there is additional scatter in stress and strain for any given orientation because of fluctuations caused by grain–grain interactions. They lead to a significant enhancement in the peak values of the grain-level stresses and strains as compared to the self-consistent model. These findings are particularly relevant if one is interested in studying the reliability of ferroelectric materials. High stresses and strains can develop inside grains in order to maintain compatibility with neighbouring grains, and they may trigger (inter-)granular failure processes. They are significantly underestimated in the self-consistent model, while also a one-element representation of individual grains is not able to pick-up the field fluctuations with sufficient accuracy. Acknowledgments This research has been carried out within the research programme ‘‘Evolution of the microstructure of materials’’ which is financially supported by the Netherlands foundation

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

663

for Fundamental Research on Matter (FOM) and the Netherlands Institute for Metals Research (NIMR).

Appendix A. Self-consistent tangent modulus In order to construct tangent moduli for the single crystal, first define, following Eqs. (6) and (7): ^ aij ¼ maij ga þ ~aij ,

ð17Þ

a a D^ i ¼ sai Pa þ D~ i ,

ð18Þ

a ^aij  hrpq D^ r , s^ apq ¼ cD pqij  a a E^ k ¼ hkij ^aij þ beki D^ i .

ð19Þ ð20Þ

The elastic compliance tensor at constant electric displacement, cD ijkl , the dielectric tensor hkij and the dielectric impermeability tensor beik are obtained by algebraic manipulation of sEijkl , d kij and keik . Huber et al. (1999) showed that tangent moduli for the single crystal, defined by (9) and (10) are given by XX D ð21Þ s^ apq ðX ab Þ1 s^ bij , cD;tan pqij ¼ cpqij  a

 htan rpq ¼ hrpq 

b

XX a

¼ bki  b;tan ki

XX a

b s^ apq ðX ab Þ1 E^ r ,

ð22Þ

b a b E^ k ðX ab Þ1 E^ i ,

ð23Þ

b

where the elements of matrix X ab are b a X ab ¼ H ab þ s^ bij ^aij þ E^ k D^ k .

(24)

In deriving the tangent moduli, summations over a and b range only over a those transformation systems that satisfy G a ¼ Gac and f_ 40. The further constraint that no volume fraction must fall below zero enforces the saturation of ferroelectric switching. The single crystal tangent moduli are then averaged over a set of grains with distinct crystallographic orientations using a self-consistent procedure (Hill, 1965; Huber et al., 1999) to obtain the overall tangent moduli defined by Eqs. (11) and (12). The self-consistent average treats each grain as a spherical inclusion embedded in a matrix of material which has the (initially unknown) tangent moduli of the composite of grains. The inclusion has uniform stress and electric field. Strain and electric displacement increments in each inclusion are related to the corresponding overall increments in the matrix using the solutions to the ellipsoidal inclusion problem developed by Eshelby (1957). Finally, equating the volume average strain and electric displacement of the grains to the macroscopic values of strain and electric displacement provides a set of overall tangent moduli. The procedure is described in detail in Huber et al. (1999).

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

664

Appendix B. Rate tangent modulus The main goal of the time integration procedure is to determine the correct values of the state variables in a configuration at time instant tnþ1 , starting with the known values of the same variables in the previous configuration, at time tn . Recognizing the fact that instabilities are mainly caused by the highly nonlinear power-law dependence (13) of the a incremental volume fraction, f_ , it can be rewritten as a linear interpolation between the values at time t and t þ Dt: a

a

Df a ¼ ½ð1  yÞf_ ðtÞ þ yf_ ðt þ DtÞDt,

(25)

where y is a parameter in the range ½0; 1 whose value is taken here to be y ¼ 0:9. A Taylor series expansion is then used to estimate the rates of volume fraction transfer at time t þ Dt, a

qf_ a a f_ ðt þ DtÞ ¼ f_ ðtÞ þ a DG a . qG

(26)

For the derivation, the incremental form of Eqs. (1) and (2) needs to be rewritten as X a _ ij  hrpq D_ r  (27) s^ apq f_ , s_ pq ¼ cD pqij e a

E_ k ¼ hkmn e_mn þ beki D_ i 

X

a a E^ k f_ .

(28)

a

The incremental volume fraction is calculated by substitution of Eqs. (27) and (28) into the incremental form of Eq. (6), and subsequent substitution in Eqs. (25) and (26) " #   X X f_b b b b Df a ¼ Dt ymDt M ab s^ ij e_ij þ E^ i D_ i þ M ab f_ , (29) b b G b where N

ab

¼d

ab

a   f_ b a þ Dtym a ðs^ bij mij ga þ e~ aij þ E^ k ðsar Pa þ D~ r ÞÞ G

(30)

and M ab is the inverse of N ab . Eqs. (27) and (28) with Eq. (29) can be rewritten as _ pq , _ _ij  htan s_ pq ¼ cD;tan rpq Dr þ s pqij e

(31)

_ _ _mn þ btan E_ k ¼ htan kmn e ki Di þ E k ,

(32)

where D cD;tan pqij ¼ cpqij  ymDt

X X f_b a

htan rpq ¼ hrpq  ymDt

b

Gb

M ab s^ apq s^ bij ,

X X f_b a

b

Gb

(33)

b

M ab s^ apq E^ r ,

(34)

ARTICLE IN PRESS A. Haug et al. / J. Mech. Phys. Solids 55 (2007) 648–665

btan ki ¼ bki  ymDt

X X f_b a

s_ pq ¼ 

XX a

b

b

G

b f_ M ab s^ apq ;

b

a b M ab E^ k E^ i ,

E_ k ¼ 

(35)

XX a

665

a b f_ M ab E^ k .

(36)

b

References Allik, H., Hughes, T.J.R., 1970. Finite element method for piezoelectric vibration. Int. J. Numer. Methods Eng. 2, 151–157. Asaro, R.J., 1983. Micromechanics of crystals and polycrystals. Adv. Appl. Mech. 23, 1–115. Cocks, A.C.F., McMeeking, R.M., 1999. A phenomenological constitutive law for the behaviour of ferroelectric ceramics. Ferroelectrics 228, 219–228. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. London A 241, 376–396. Haug, A., Onck, P.R., Van der Giessen, E., 2006. Development of inter- and intragranular stresses during switching of ferroelectric polycrystals. Int. J. Solids Struct., in press, doi:10.1016/j.ijsolstr.2006.07.024. Hill, R., 1965. A self-consistent mechanics of composite materials. J. Mech. Phys. Solids 13, 213–222. Huber, J.E., Fleck, N.A., 2001. Multi-axial electrical switching of a ferroelectric: theory versus experiment. J. Mech. Phys. Solids 49 (4), 785–811. Huber, J.E., Fleck, N.A., Landis, C.M., McMeeking, R.M., 1999. A constitutive model for ferroelectric polycrystals. J. Mech. Phys. Solids 47 (8), 1663–1697. Jaffe, B., Cook, W.R., Jaffe, H., 1971. Piezoelectric Ceramics. Academic Press, London, New York. Kamlah, M., Liskowsky, A.C., McMeeking, R.M., Balke, H., 2005. Finite element simulation of a polycrystalline ferroelectric based on a multidomain single crystal switching model. Int. J. Solids Struct. 42 (9–10), 2949–2964. Kamlah, M., Tsakmakis, C., 1999. Phenomenological modeling of the non-linear electro-mechanical coupling in ferroelectrics. Int. J. Solids Struct. 36, 669–695. Landis, C.M., 2002. Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics. J. Mech. Phys. Solids 50 (1), 127–152. McMeeking, R.M., Landis, C.M., 2002. A phenomenological multi-axial constitutive law for switching in polycrystalline ferroelectric ceramics. Int. J. Eng. Sci. 40 (14), 1553–1577. Pierce, D., Shih, C.F., Needleman, A., 1984. A tangent modulus method for rate dependent solids. Comput. Struct. 18 (9–10), 875–887. Wang, J., Shi, S.Q., Chen, L.Q., Li, Y.L., Zhang, T.Y., 2004. Phase field simulations of ferroelectric/ferroelastic polarization switching. Acta Mater. 52 (3), 749–764.