- Email: [email protected]
Viscoplastic behaviour of perovskite type ferroelectrics
Materials Science and Engineering B 118 (2005) 7–11
Viscoplastic behaviour of perovskite type ferroelectrics A.Yu. Belov ∗ , W.S. Kreher Technische Universit¨at Dresden, Institut f¨ur Werkstoffwissenschaft, 01062 Dresden, Germany
Abstract We present a viscoplastic (rate-dependent) constitutive model describing both hysteretic and sub-coercive behaviour of ferroelectric ceramics. Viscoplastic models employ rate equations for the volume fractions of orientation variants along with a distribution function for grain orientations, providing a statistical macroscopic description for the processes of polarization reversal in polycrystalline ferroelectrics. Here attention is focused on two viscoplastic models for perovskite type ferroelectrics undergoing a cubic-to-tetragonal phase transition. Both models incorporate 90◦ and 180◦ domain switching, but differ in a number of representative domain orientations, which is 6 and 42, respectively. Different parameterisations for the models were tested, with a special emphasis on orientation effects in the behaviour of poled ferroelectrics under multi-axial loading as well as on rate effects. Simulation of the polarization rotation in soft PZT ceramics by electric field applied at an angle to the original poling direction show good agreement between the predictions of the viscoplastic model and experimental data. © 2004 Elsevier B.V. All rights reserved. Keywords: Ferroelectrics; Ceramics; Domains; Viscoplasticity
1. Introduction The ferroelectric perovskite oxides ABO3 are used in a wide range of functional ceramics and form the materials base for many electronic applications. It is widely recognized that elastic, dielectric, and piezoelectric properties of ferroelectric ceramics of this family such as lead zirconate titanate, Pb(Zr,Ti)O3 (PZT), are largely due to irreversible processes associated with the displacement of domain walls. Most of the polarization reversal modelling in ferroelectric ceramics is based on switching criteria formulated in terms of a critical field or work required for domain reorientation [1–3]. Upon application of a field below the coercive field (Ec ) a near-linear reversible polarization is induced and the remanent polarization occurs only as Ec is exceeded. In the case of multi-axial loading this approach is modified by introducing a yield or switching surface (in the combined electric field and mechanical stress space), within which switching does not ∗ Corresponding author. Present address: Institute of Crystallography RAS, Moscow, Russia. Tel.: +49 351 463 31 409; fax: +49 351 463 31 422. E-mail address: [email protected] (A.Yu. Belov).
0921-5107/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.mseb.2004.12.057
occur [3,4]. However, a number of experimental data shows [5–9] that irreversible processes in ferroelectric ceramics are important also at sub-coercive fields E < Ec . Therefore, it is an attractive goal to formulate a constitutive model, which would describe the response of a ferroelectric ceramic not only at E ≥ Ec but also in the sub-coercive region E < Ec . The obvious requirements to such a model are as follows. First, it must capture the basic physical processes underlying the polarization reversal. Secondly, it should remain computationally efficient to be implemented into finite element models of real scale devices such as sensors or actuators. In this study we approach the problem by using a viscoplastic (rate-dependent) model without a switching criterion. Its functional form is similar to that used in [3]. The difference between the two models consists in the number of representative polarization directions used to reproduce the isotropic response of the unpoled polycrystalline ferroelectrics. The viscoplastic models also incorporate the rate dependence of the domain reversal. The polarization reversal in perovskite ferroelectrics was reported [9] for an extremely broad time range (10−8 < t < 102 s), indicating a variety of physical mechanisms limiting the domain switching process. Some experimental data [6–8] for soft PZT ceramics show remarkable rate effects even at low
8
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
field frequencies of 0.1–1.0 Hz. Correspondingly, the time scale for the rate dependent slow processes can approach seconds.
S are the wall velocity and average area, respectively. For Vrs an Arrhenius type dependency is adopted g(E, σ) Vrs = Vc exp − , (5) kT
2. Viscoplastic model of a ferroelectric ceramic
where g is the Gibbs energy required for a domain wall to overcome an obstacle (activation energy) and Vc is the wall velocity at zero activation energy. The activation energy depends on both the nature of pinning defects and their distribution. In principle, it could be found from atomistic calculations. However, for formulation of an empirical constitutive model, it is sufficient to employ some approximate expressions for g, similar to that used in dislocation plasticity. Here, we suggest that whatever loading scheme is applied (uni- or multi-axial), the activation energy depends on the electric field and stress tensor only via the driving force frs = g(r) – g(s) > 0. Following [12], one has p frs g(E, σ) = g0 1 − . (6) fc
Unlike the time-independent models [1–3], a viscoplastic model does not imply a critical field for the domain switching. It takes into account the thermally activated nature of the domain wall motion in the random potential relief of pinning centres. Therefore, at nonzero temperature domain walls can be displaced irreversibly at any electric field, whatever small it may be. We consider a ferroelectric material with n distinct orientation states (variants). For instance, n is equal to six and eight in tetragonal and rhombohedral materials, respectively. Each variant r is characterized by its volume fraction ξ r , as well as its spontaneous strain εs(r) and polarization Ps(r) . The volume averages of the spontaneous strain and polarization in a polydomain crystal (or grain) have the form s n s(r) ε ε = ξ. (1) Ps P s(r) r r=1
The corresponding constitutive relation for the crystal (grain) is derived under an assumption (known as the Reuss method) that the electric field E and stress σ are identical in all variants s n (r) ε S σ ε d (r)T = + ξr . (2) D Ps E d (r) κ(r) r=1
Here S(r) , d(r) , and κ(r) are the tensors of elastic, piezoelectric, and dielectric constants of the r variant. Under zero electric field and stress all n variants possess the same Gibbs free energy g(r) = g0 (r = 1, . . ., n). Application of a field eliminates this energetic degeneracy and the Gibbs energy of the variant r becomes 1 g(r) = g0 − σ : εs(r) − P s(r) E − σ : S (r) 2 1 : σ − Ed (r) : σ − Ek(r) E. (3) 2 The fast domain switching typical for thin films and proceeding at the time scale of 10–100 ns is essentially controlled by nucleation kinetics and is limited by the formation time for a nucleus of a critical size. The formalisms to describe the fast processes are given in [10,11]. In the case of slow processes controlled by dynamics of domain walls, the decrease in the volume fraction of the r variant with time (transformation rate) due to its switching into a state s of a lower Gibbs energy can be expressed by an Orowan type equation v˙ rs = ρrs S Vrs .
(4)
Here ρrs is the density of the mobile domain walls separating the r- and s-type domains (transformation system), Vrs and
Here, g0 is the activation energy at zero driving force, fc and p characterise the obstacle strength and shape, respectively. The parameters g0 , fc , and p depend on the domain wall type and details of its interaction with pinning centres. Tetragonal ferroelectrics with 180◦ and 90◦ domain walls need two sets of their adjustable parameters. In turn, rhombohedral crystals allow for 180◦ , 109◦ , and 71◦ polarization switching and three parameter sets are required. Some ceramics, like PZT with compositions near the morphotropic phase boundary, admit a coexistence of the tetragonal and rhombohedral phases and hence the transformation systems of both phases have to be incorporated into the model. In the limit of frs → fc , the activation energy (6) in combination with (5) and (4) yields to the power (viscoplastic) law for the transformation rate, m frs . (7) v˙ rs = v˙ c fc The viscoplastic index m is proportional to g0 /kT and normally largely exceeds unity. Hence, Eq. (7) provides an estimate for v˙ rs in the whole range 0 ≤ frs ≤ fc . In fact, in the athermal regime (frs > fc ) the viscoplastic law must be replaced by a viscous dependency like v˙ rc ∼ frs . Nevertheless, we adopt it also for the athermal regime to describe high transformation rates in this case. It will be further assumed that the factor v˙ c , which generally depends on the microstructure, is expressed in terms of the volume fraction ξ r of the r-domains. According to [3], saturation of a transformation system rs is controlled by a simple switching function as m α frs ξr v˙ rs = v˙ c . (8) fc ξ0 The re-defined parameter v˙ c in Eq. (8) is now independent of the microstructure. So far, our consideration addressed the case of a single polydomain grain with n orientation variants.
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
9
Fig. 1. (a) A stereographic projection of the polarization directions in the models I and II. Six polarization orientations Ps(r) (r = 1, . . ., 6) of model I are situated along the cubic axes of the paraelectric phase. In model II six additional sets of polarization Ps(r,n) (r, n = 1, . . ., 6) are situated in the diagonal planes of the cube. Only two sets, Ps(r,1) and Ps(r,2) , are shown. The remaining four sets are related to other diagonal planes. (b) The electric field orientation with respect to the coordinate system (x1 , x2 , x3 ).
In ceramics one has a random distribution of grain orientations described by a continuous function and Eqs. (1) and (2) are to be modified to take into account the additional averaging over grain orientations. Here two discrete approximations (models I and II) of this function for tetragonal ceramics will be used. While model I is based on one system of six variants and corresponds to one possible grain orientation in a paraelectric cubic phase, model II employs seven systems of six variants in each. The total number of accessible polarization directions in model II is 42. However, there is no polarization switching between orientation variants, which belong to different systems. The relative arrangement of domains in both models is shown in Fig. 1.
This effect is mainly due to a larger number of transformation systems for 90◦ switching in comparison with only one system for direct 180◦ switching in this model.
4. Frequency dependence and sub-coercive behaviour The simulated polarization hysteresis curves for the field amplitudes of 2.0 and 0.5 MV/m are shown in Figs. 3 and 4,
3. Uni-axial loading The effect of the model parameters on polarization switching under a uni-axial electrical loading E3 = Esin ωt of large amplitude was first simulated for model I. Both viscoplastic models are suitably calibrated by introducing (90) (90) m two parameters ω(90) = v˙ c (P s E0 /fc ) and ω(180) = m (180) (180) v˙ c (2P s E0 /fc ) with a dimension of frequency. These parameters characterize the effective mobilities of the 180◦ and 90◦ domain walls under a reference electric field with an amplitude E0 of 2 MV/m. The dependence of the wall mobilities and the corresponding transformation rates in Eq. (8) on an electric field amplitude E is given by a powerlaw ∼(E/E0 )m . Fig. 2 illustrates the dependence of polarization hysteresis under a quasi-static field with amplitude of 2.0 MV/m and a frequency of 0.01 Hz on the parameter α. The simulation was performed for the frequency ω(90) of 20 Hz and two frequency ω(180) values of 20 and 1 Hz. Fig. 2 shows that the coercive field Ec decreases with increasing α, indicating the gradual onset of domain switching. However, larger values of α impede the polarization saturation and reduce the remanent polarization. Interestingly, the simulation results are nearly independent of 180◦ domain switching.
Fig. 2. Dependence of the polarization hysteresis loops on the parameter α. The sinusoidal electric field has an amplitude of 2.0 MV/m and a frequency of 0.01 Hz. The simulation was performed for two parameter sets: (a) m = 8, ω(90) = 20 Hz, ω(180) = 20 and (b) m = 8, ω(90) = 20 Hz, ω(180) = 1 Hz.
10
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
5. Multi-axial loading
Fig. 3. Polarization vs. sinusoidal electric field hysteresis curves for frequencies of 0.1–100 Hz. The field has an amplitude of 2.0 MV/m. Parameter set: α = 3.0, ω(90) = 10 Hz, ω(180) = 10 Hz.
respectively. The simulation was performed for the frequencies of 0.1, 10 and 100 Hz, using the model I. Its parameterisation corresponds to the experimental data [8] of the commercial soft PZT piezoceramic PIC151 in the quasi-static regime with a frequency of 0.01 Hz. This material is known to possess a pronounced frequency dependence of the polarization hysteresis already at a frequency of 1 Hz. The following parameter values were used: m = 8, α = 3, ω(90) = 10 Hz and ω(180) = 10 Hz. In agreement with experiment [8] the simulation predicts incomplete polarization saturation at the maximum field for the loading rates exceeding 1 Hz, with the coercive field Ec gradually increasing with frequency. In the frequency range of 10–100 Hz the remanent polarization sharply decreases with frequency. Fig. 4 shows the response of a ferroelectric to a weak sub-coercive field when the viscoplastic index m is modified. Although polarization varies with the field mostly linearly, contribution of irreversible processes to the polarization can be enhanced, decreasing the viscoplastic index. This is in accordance with the essentially thermoactivative nature of the ferroelectric switching under weak fields
Fig. 4. Dependence of polarization vs. electric field curves on the viscoplastic index m. The sinusoidal field has an amplitude of 0.5 MV/m and a frequency of 0.1 Hz. Parameter set: ␣ = 3.0, ω(90) = 10 Hz, ω(180) = 10 Hz.
Since the electric field and stress components enter the viscoplastic model via the driving forces, the model can be easily evaluated for the case of multi-axial loading with only minor changes in the parameters. As an example, we simulate the polarization rotation tests [3,4], which were performed for both hard (PZT-4D) and soft (PZT-5H) ferroelectric ceramics. An unpoled ferroelectric ceramic was first poled by a constant electric field of 1.5 MV/m for 100 s. Then the sample was subjected to an electric field at an angle θ to the polarization direction for 30 s. The polarization rotation was measured for fields in the range of 0.0–2.0 MV/m. Fig. 5 shows results of a simulation performed using model II under exactly the same conditions as in the experiments [3,4]. The change D in the dielectric displacement (along E) was calculated as a function of θ and φ, as shown in Fig. 1b. Comparison of the simulation results with experiment shows that the current parameterisation of the viscoplastic model corresponds to the soft PZT-5H ceramics, which exhibit pronounced anisotropy of the ferroelectric switching for field amplitudes above Ec . However, one should mention that the model II nicely reproduces the isotropic response of unpoled ceramics. The simulation results for two values of φ show that the polarization rotation is nearly independent of the azimuth as it would be expected for isotropic polycrystalline materials.
Fig. 5. The simulated change in electric displacement D vs. electric field E inclined at an angle θ to the polarization direction. Parameter set: m = 8, α = 3, ω(90) = 20 Hz, and ω(180) = 1 Hz.
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
11
voltage as well as upon unloading are illustrated in Fig. 6. The poling was performed by a sinusoidal electric field with a frequency of 0.1 Hz. The total loading-unloading time π/ω corresponds to 31.4 s. The simulation shows some additional material poling during uploading, however, for the geometrical parameters used, the thin layer above the pore remains non-polarized.
7. Conclusions
Fig. 6. Distribution of the spontaneous polarization Ps 3 in a double-layer film with a pore. The values of Ps 3 are given in C/m2 . The thick film and thin layer have a thickness of 100 m and 1 m, respectively. The pore diameter is 10 m and its depth below the thick film surface is 0.1 m. The bottom surface was set to zero potential and the upper surface to the sinusoidal potential with amplitude of −202 V, corresponding to the average electric field of 2.0 MV/m in the double-layer film. The polarization distribution is shown at the maximum voltage (left) and upon unloading (right).
6. Finite element implementation of the model Computational efficiency and a possibility of implementation into finite element models of ferroelectric devices possessing complex shapes is an important requirement to a constitutive material model. Here we apply the viscoplastic model to the finite element analysis of the poling process in PZT double-layer films containing microstructural inhomogeneities. Such films consist of a thin upper layer and a thick substrate film and are of interest for application as a ferroelectric printing form. It is expected that non-polarizable defects like pores in the substrate film hinder the poling of the thin layer, deteriorating the printing quality. As an initially unpoled double-layer film with a pore is subjected to an electric loading perpendicular to the film plane, the poling process near the pore is prohibited by the electrical boundary conditions on its surface. Therefore, the resulting inhomogeneous polarization distribution is dependent on the geometrical parameters such as the pore size and depth, and also on the thin layer thickness. The non-linear analysis of the poling process was performed using the viscoplastic model I with the parameter set m = 8, α = 3, ω(90) = 20 Hz and ω(180) = 1 Hz. The volume fractions ξ r (R, Z) of the orientation variants now depend on the spatial coordinates R and Z = x3 due to the axial symmetry of the problem. The poling process is simulated by integrating the rate equations in time numerically. At every time step the electric field and mechanical stress inside the film are found by solving an axially-symmetric electroelastostatic problem. This is done by means of a commercial finite element code FlexPDE [13]. The set of material parameters is taken from [14]. The distributions of the spontaneous polarization Ps 3 in the double-layer film at the maximum
A viscoplastic constitutive model for ferroelectric ceramics has been presented to include irreversible ferroelectric switching at low and high loading levels. Viscoplastic modelling is demonstrated to be a computationally efficient method for the analysis of various processes in ferroelectric ceramics, including rate and temperature dependence of the polarization hysteresis. A viscoplastic model with 42 representative domain orientations is shown to be useful for simulations of ferroelectric ceramics under multi-axial loading conditions. It can reproduce the features of the polarization rotation tests and in addition exhibits isotropic behaviour, which is typical for polycrystals. The computational efficiency of the viscoplastic model was demonstrated by its implementation into a 3D finite element model of the poling process in ferroelectric multilayer films containing microstructural inhomogeneities.
Acknowledgement Support by Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References [1] S.C. Hwang, C.S. Lynch, R.M. McMeeking, Acta Metall. Mater. 43 (1995) 2073. [2] S.C. Hwang, J.E. Huber, R.M. McMeeking, N.A. Fleck, J. Appl. Phys. 84 (1998) 1530. [3] J.E. Huber, N.A. Fleck, J. Mech. Phys. Solids 49 (2001) 785. [4] J. Shieh, J.E. Huber, N.A. Fleck, Acta Mater. 51 (2003) 6123. [5] D. Damjanovic, J. Appl. Phys. 82 (1997) 1788. [6] Y.H. Chen, D. Viehland, Appl. Phys. Lett. 77 (2000) 133. [7] D. Viehland, Y.H. Chen, J. Appl. Phys. 88 (2000) 6696. [8] D. Zhou, M. Kamlah, D. Munz, in: C.S. Lynch (Ed.) Smart Structures and Materials 2001: Active Materials: Behavior and Mechanics, vol. 4333, SPIE, 2001, p. 64. [9] C. Jullian, J.F. Li, D. Viehland, Appl. Phys. Lett. 83 (2003) 1196. [10] E. Fatuzzo, Phys. Rev. 116 (1962) 1999. [11] Y. Ishibashi, I. Takagi, J. Phys. Soc. Jpn. 31 (1971) 506. [12] U.F. Kocks, A.S. Argon, M.F. Ashby, Progr. Mater. Sci. 19 (1975) 1. [13] G. Backstr¨om, Fields of Physics by Finite Element Analysis, GB Publishing, Malm¨o, 2000. [14] J. R¨odel, W. Beckert, W.S. Kreher, J. Eur. Ceram. Soc. 24 (2004) 1647.
Viscoplastic behaviour of perovskite type ferroelectrics A.Yu. Belov ∗ , W.S. Kreher Technische Universit¨at Dresden, Institut f¨ur Werkstoffwissenschaft, 01062 Dresden, Germany
Abstract We present a viscoplastic (rate-dependent) constitutive model describing both hysteretic and sub-coercive behaviour of ferroelectric ceramics. Viscoplastic models employ rate equations for the volume fractions of orientation variants along with a distribution function for grain orientations, providing a statistical macroscopic description for the processes of polarization reversal in polycrystalline ferroelectrics. Here attention is focused on two viscoplastic models for perovskite type ferroelectrics undergoing a cubic-to-tetragonal phase transition. Both models incorporate 90◦ and 180◦ domain switching, but differ in a number of representative domain orientations, which is 6 and 42, respectively. Different parameterisations for the models were tested, with a special emphasis on orientation effects in the behaviour of poled ferroelectrics under multi-axial loading as well as on rate effects. Simulation of the polarization rotation in soft PZT ceramics by electric field applied at an angle to the original poling direction show good agreement between the predictions of the viscoplastic model and experimental data. © 2004 Elsevier B.V. All rights reserved. Keywords: Ferroelectrics; Ceramics; Domains; Viscoplasticity
1. Introduction The ferroelectric perovskite oxides ABO3 are used in a wide range of functional ceramics and form the materials base for many electronic applications. It is widely recognized that elastic, dielectric, and piezoelectric properties of ferroelectric ceramics of this family such as lead zirconate titanate, Pb(Zr,Ti)O3 (PZT), are largely due to irreversible processes associated with the displacement of domain walls. Most of the polarization reversal modelling in ferroelectric ceramics is based on switching criteria formulated in terms of a critical field or work required for domain reorientation [1–3]. Upon application of a field below the coercive field (Ec ) a near-linear reversible polarization is induced and the remanent polarization occurs only as Ec is exceeded. In the case of multi-axial loading this approach is modified by introducing a yield or switching surface (in the combined electric field and mechanical stress space), within which switching does not ∗ Corresponding author. Present address: Institute of Crystallography RAS, Moscow, Russia. Tel.: +49 351 463 31 409; fax: +49 351 463 31 422. E-mail address: [email protected] (A.Yu. Belov).
0921-5107/$ – see front matter © 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.mseb.2004.12.057
occur [3,4]. However, a number of experimental data shows [5–9] that irreversible processes in ferroelectric ceramics are important also at sub-coercive fields E < Ec . Therefore, it is an attractive goal to formulate a constitutive model, which would describe the response of a ferroelectric ceramic not only at E ≥ Ec but also in the sub-coercive region E < Ec . The obvious requirements to such a model are as follows. First, it must capture the basic physical processes underlying the polarization reversal. Secondly, it should remain computationally efficient to be implemented into finite element models of real scale devices such as sensors or actuators. In this study we approach the problem by using a viscoplastic (rate-dependent) model without a switching criterion. Its functional form is similar to that used in [3]. The difference between the two models consists in the number of representative polarization directions used to reproduce the isotropic response of the unpoled polycrystalline ferroelectrics. The viscoplastic models also incorporate the rate dependence of the domain reversal. The polarization reversal in perovskite ferroelectrics was reported [9] for an extremely broad time range (10−8 < t < 102 s), indicating a variety of physical mechanisms limiting the domain switching process. Some experimental data [6–8] for soft PZT ceramics show remarkable rate effects even at low
8
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
field frequencies of 0.1–1.0 Hz. Correspondingly, the time scale for the rate dependent slow processes can approach seconds.
S are the wall velocity and average area, respectively. For Vrs an Arrhenius type dependency is adopted g(E, σ) Vrs = Vc exp − , (5) kT
2. Viscoplastic model of a ferroelectric ceramic
where g is the Gibbs energy required for a domain wall to overcome an obstacle (activation energy) and Vc is the wall velocity at zero activation energy. The activation energy depends on both the nature of pinning defects and their distribution. In principle, it could be found from atomistic calculations. However, for formulation of an empirical constitutive model, it is sufficient to employ some approximate expressions for g, similar to that used in dislocation plasticity. Here, we suggest that whatever loading scheme is applied (uni- or multi-axial), the activation energy depends on the electric field and stress tensor only via the driving force frs = g(r) – g(s) > 0. Following [12], one has p frs g(E, σ) = g0 1 − . (6) fc
Unlike the time-independent models [1–3], a viscoplastic model does not imply a critical field for the domain switching. It takes into account the thermally activated nature of the domain wall motion in the random potential relief of pinning centres. Therefore, at nonzero temperature domain walls can be displaced irreversibly at any electric field, whatever small it may be. We consider a ferroelectric material with n distinct orientation states (variants). For instance, n is equal to six and eight in tetragonal and rhombohedral materials, respectively. Each variant r is characterized by its volume fraction ξ r , as well as its spontaneous strain εs(r) and polarization Ps(r) . The volume averages of the spontaneous strain and polarization in a polydomain crystal (or grain) have the form s n s(r) ε ε = ξ. (1) Ps P s(r) r r=1
The corresponding constitutive relation for the crystal (grain) is derived under an assumption (known as the Reuss method) that the electric field E and stress σ are identical in all variants s n (r) ε S σ ε d (r)T = + ξr . (2) D Ps E d (r) κ(r) r=1
Here S(r) , d(r) , and κ(r) are the tensors of elastic, piezoelectric, and dielectric constants of the r variant. Under zero electric field and stress all n variants possess the same Gibbs free energy g(r) = g0 (r = 1, . . ., n). Application of a field eliminates this energetic degeneracy and the Gibbs energy of the variant r becomes 1 g(r) = g0 − σ : εs(r) − P s(r) E − σ : S (r) 2 1 : σ − Ed (r) : σ − Ek(r) E. (3) 2 The fast domain switching typical for thin films and proceeding at the time scale of 10–100 ns is essentially controlled by nucleation kinetics and is limited by the formation time for a nucleus of a critical size. The formalisms to describe the fast processes are given in [10,11]. In the case of slow processes controlled by dynamics of domain walls, the decrease in the volume fraction of the r variant with time (transformation rate) due to its switching into a state s of a lower Gibbs energy can be expressed by an Orowan type equation v˙ rs = ρrs S Vrs .
(4)
Here ρrs is the density of the mobile domain walls separating the r- and s-type domains (transformation system), Vrs and
Here, g0 is the activation energy at zero driving force, fc and p characterise the obstacle strength and shape, respectively. The parameters g0 , fc , and p depend on the domain wall type and details of its interaction with pinning centres. Tetragonal ferroelectrics with 180◦ and 90◦ domain walls need two sets of their adjustable parameters. In turn, rhombohedral crystals allow for 180◦ , 109◦ , and 71◦ polarization switching and three parameter sets are required. Some ceramics, like PZT with compositions near the morphotropic phase boundary, admit a coexistence of the tetragonal and rhombohedral phases and hence the transformation systems of both phases have to be incorporated into the model. In the limit of frs → fc , the activation energy (6) in combination with (5) and (4) yields to the power (viscoplastic) law for the transformation rate, m frs . (7) v˙ rs = v˙ c fc The viscoplastic index m is proportional to g0 /kT and normally largely exceeds unity. Hence, Eq. (7) provides an estimate for v˙ rs in the whole range 0 ≤ frs ≤ fc . In fact, in the athermal regime (frs > fc ) the viscoplastic law must be replaced by a viscous dependency like v˙ rc ∼ frs . Nevertheless, we adopt it also for the athermal regime to describe high transformation rates in this case. It will be further assumed that the factor v˙ c , which generally depends on the microstructure, is expressed in terms of the volume fraction ξ r of the r-domains. According to [3], saturation of a transformation system rs is controlled by a simple switching function as m α frs ξr v˙ rs = v˙ c . (8) fc ξ0 The re-defined parameter v˙ c in Eq. (8) is now independent of the microstructure. So far, our consideration addressed the case of a single polydomain grain with n orientation variants.
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
9
Fig. 1. (a) A stereographic projection of the polarization directions in the models I and II. Six polarization orientations Ps(r) (r = 1, . . ., 6) of model I are situated along the cubic axes of the paraelectric phase. In model II six additional sets of polarization Ps(r,n) (r, n = 1, . . ., 6) are situated in the diagonal planes of the cube. Only two sets, Ps(r,1) and Ps(r,2) , are shown. The remaining four sets are related to other diagonal planes. (b) The electric field orientation with respect to the coordinate system (x1 , x2 , x3 ).
In ceramics one has a random distribution of grain orientations described by a continuous function and Eqs. (1) and (2) are to be modified to take into account the additional averaging over grain orientations. Here two discrete approximations (models I and II) of this function for tetragonal ceramics will be used. While model I is based on one system of six variants and corresponds to one possible grain orientation in a paraelectric cubic phase, model II employs seven systems of six variants in each. The total number of accessible polarization directions in model II is 42. However, there is no polarization switching between orientation variants, which belong to different systems. The relative arrangement of domains in both models is shown in Fig. 1.
This effect is mainly due to a larger number of transformation systems for 90◦ switching in comparison with only one system for direct 180◦ switching in this model.
4. Frequency dependence and sub-coercive behaviour The simulated polarization hysteresis curves for the field amplitudes of 2.0 and 0.5 MV/m are shown in Figs. 3 and 4,
3. Uni-axial loading The effect of the model parameters on polarization switching under a uni-axial electrical loading E3 = Esin ωt of large amplitude was first simulated for model I. Both viscoplastic models are suitably calibrated by introducing (90) (90) m two parameters ω(90) = v˙ c (P s E0 /fc ) and ω(180) = m (180) (180) v˙ c (2P s E0 /fc ) with a dimension of frequency. These parameters characterize the effective mobilities of the 180◦ and 90◦ domain walls under a reference electric field with an amplitude E0 of 2 MV/m. The dependence of the wall mobilities and the corresponding transformation rates in Eq. (8) on an electric field amplitude E is given by a powerlaw ∼(E/E0 )m . Fig. 2 illustrates the dependence of polarization hysteresis under a quasi-static field with amplitude of 2.0 MV/m and a frequency of 0.01 Hz on the parameter α. The simulation was performed for the frequency ω(90) of 20 Hz and two frequency ω(180) values of 20 and 1 Hz. Fig. 2 shows that the coercive field Ec decreases with increasing α, indicating the gradual onset of domain switching. However, larger values of α impede the polarization saturation and reduce the remanent polarization. Interestingly, the simulation results are nearly independent of 180◦ domain switching.
Fig. 2. Dependence of the polarization hysteresis loops on the parameter α. The sinusoidal electric field has an amplitude of 2.0 MV/m and a frequency of 0.01 Hz. The simulation was performed for two parameter sets: (a) m = 8, ω(90) = 20 Hz, ω(180) = 20 and (b) m = 8, ω(90) = 20 Hz, ω(180) = 1 Hz.
10
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
5. Multi-axial loading
Fig. 3. Polarization vs. sinusoidal electric field hysteresis curves for frequencies of 0.1–100 Hz. The field has an amplitude of 2.0 MV/m. Parameter set: α = 3.0, ω(90) = 10 Hz, ω(180) = 10 Hz.
respectively. The simulation was performed for the frequencies of 0.1, 10 and 100 Hz, using the model I. Its parameterisation corresponds to the experimental data [8] of the commercial soft PZT piezoceramic PIC151 in the quasi-static regime with a frequency of 0.01 Hz. This material is known to possess a pronounced frequency dependence of the polarization hysteresis already at a frequency of 1 Hz. The following parameter values were used: m = 8, α = 3, ω(90) = 10 Hz and ω(180) = 10 Hz. In agreement with experiment [8] the simulation predicts incomplete polarization saturation at the maximum field for the loading rates exceeding 1 Hz, with the coercive field Ec gradually increasing with frequency. In the frequency range of 10–100 Hz the remanent polarization sharply decreases with frequency. Fig. 4 shows the response of a ferroelectric to a weak sub-coercive field when the viscoplastic index m is modified. Although polarization varies with the field mostly linearly, contribution of irreversible processes to the polarization can be enhanced, decreasing the viscoplastic index. This is in accordance with the essentially thermoactivative nature of the ferroelectric switching under weak fields
Fig. 4. Dependence of polarization vs. electric field curves on the viscoplastic index m. The sinusoidal field has an amplitude of 0.5 MV/m and a frequency of 0.1 Hz. Parameter set: ␣ = 3.0, ω(90) = 10 Hz, ω(180) = 10 Hz.
Since the electric field and stress components enter the viscoplastic model via the driving forces, the model can be easily evaluated for the case of multi-axial loading with only minor changes in the parameters. As an example, we simulate the polarization rotation tests [3,4], which were performed for both hard (PZT-4D) and soft (PZT-5H) ferroelectric ceramics. An unpoled ferroelectric ceramic was first poled by a constant electric field of 1.5 MV/m for 100 s. Then the sample was subjected to an electric field at an angle θ to the polarization direction for 30 s. The polarization rotation was measured for fields in the range of 0.0–2.0 MV/m. Fig. 5 shows results of a simulation performed using model II under exactly the same conditions as in the experiments [3,4]. The change D in the dielectric displacement (along E) was calculated as a function of θ and φ, as shown in Fig. 1b. Comparison of the simulation results with experiment shows that the current parameterisation of the viscoplastic model corresponds to the soft PZT-5H ceramics, which exhibit pronounced anisotropy of the ferroelectric switching for field amplitudes above Ec . However, one should mention that the model II nicely reproduces the isotropic response of unpoled ceramics. The simulation results for two values of φ show that the polarization rotation is nearly independent of the azimuth as it would be expected for isotropic polycrystalline materials.
Fig. 5. The simulated change in electric displacement D vs. electric field E inclined at an angle θ to the polarization direction. Parameter set: m = 8, α = 3, ω(90) = 20 Hz, and ω(180) = 1 Hz.
A.Yu. Belov, W.S. Kreher / Materials Science and Engineering B 118 (2005) 7–11
11
voltage as well as upon unloading are illustrated in Fig. 6. The poling was performed by a sinusoidal electric field with a frequency of 0.1 Hz. The total loading-unloading time π/ω corresponds to 31.4 s. The simulation shows some additional material poling during uploading, however, for the geometrical parameters used, the thin layer above the pore remains non-polarized.
7. Conclusions
Fig. 6. Distribution of the spontaneous polarization Ps 3 in a double-layer film with a pore. The values of Ps 3 are given in C/m2 . The thick film and thin layer have a thickness of 100 m and 1 m, respectively. The pore diameter is 10 m and its depth below the thick film surface is 0.1 m. The bottom surface was set to zero potential and the upper surface to the sinusoidal potential with amplitude of −202 V, corresponding to the average electric field of 2.0 MV/m in the double-layer film. The polarization distribution is shown at the maximum voltage (left) and upon unloading (right).
6. Finite element implementation of the model Computational efficiency and a possibility of implementation into finite element models of ferroelectric devices possessing complex shapes is an important requirement to a constitutive material model. Here we apply the viscoplastic model to the finite element analysis of the poling process in PZT double-layer films containing microstructural inhomogeneities. Such films consist of a thin upper layer and a thick substrate film and are of interest for application as a ferroelectric printing form. It is expected that non-polarizable defects like pores in the substrate film hinder the poling of the thin layer, deteriorating the printing quality. As an initially unpoled double-layer film with a pore is subjected to an electric loading perpendicular to the film plane, the poling process near the pore is prohibited by the electrical boundary conditions on its surface. Therefore, the resulting inhomogeneous polarization distribution is dependent on the geometrical parameters such as the pore size and depth, and also on the thin layer thickness. The non-linear analysis of the poling process was performed using the viscoplastic model I with the parameter set m = 8, α = 3, ω(90) = 20 Hz and ω(180) = 1 Hz. The volume fractions ξ r (R, Z) of the orientation variants now depend on the spatial coordinates R and Z = x3 due to the axial symmetry of the problem. The poling process is simulated by integrating the rate equations in time numerically. At every time step the electric field and mechanical stress inside the film are found by solving an axially-symmetric electroelastostatic problem. This is done by means of a commercial finite element code FlexPDE [13]. The set of material parameters is taken from [14]. The distributions of the spontaneous polarization Ps 3 in the double-layer film at the maximum
A viscoplastic constitutive model for ferroelectric ceramics has been presented to include irreversible ferroelectric switching at low and high loading levels. Viscoplastic modelling is demonstrated to be a computationally efficient method for the analysis of various processes in ferroelectric ceramics, including rate and temperature dependence of the polarization hysteresis. A viscoplastic model with 42 representative domain orientations is shown to be useful for simulations of ferroelectric ceramics under multi-axial loading conditions. It can reproduce the features of the polarization rotation tests and in addition exhibits isotropic behaviour, which is typical for polycrystals. The computational efficiency of the viscoplastic model was demonstrated by its implementation into a 3D finite element model of the poling process in ferroelectric multilayer films containing microstructural inhomogeneities.
Acknowledgement Support by Deutsche Forschungsgemeinschaft is gratefully acknowledged.
References [1] S.C. Hwang, C.S. Lynch, R.M. McMeeking, Acta Metall. Mater. 43 (1995) 2073. [2] S.C. Hwang, J.E. Huber, R.M. McMeeking, N.A. Fleck, J. Appl. Phys. 84 (1998) 1530. [3] J.E. Huber, N.A. Fleck, J. Mech. Phys. Solids 49 (2001) 785. [4] J. Shieh, J.E. Huber, N.A. Fleck, Acta Mater. 51 (2003) 6123. [5] D. Damjanovic, J. Appl. Phys. 82 (1997) 1788. [6] Y.H. Chen, D. Viehland, Appl. Phys. Lett. 77 (2000) 133. [7] D. Viehland, Y.H. Chen, J. Appl. Phys. 88 (2000) 6696. [8] D. Zhou, M. Kamlah, D. Munz, in: C.S. Lynch (Ed.) Smart Structures and Materials 2001: Active Materials: Behavior and Mechanics, vol. 4333, SPIE, 2001, p. 64. [9] C. Jullian, J.F. Li, D. Viehland, Appl. Phys. Lett. 83 (2003) 1196. [10] E. Fatuzzo, Phys. Rev. 116 (1962) 1999. [11] Y. Ishibashi, I. Takagi, J. Phys. Soc. Jpn. 31 (1971) 506. [12] U.F. Kocks, A.S. Argon, M.F. Ashby, Progr. Mater. Sci. 19 (1975) 1. [13] G. Backstr¨om, Fields of Physics by Finite Element Analysis, GB Publishing, Malm¨o, 2000. [14] J. R¨odel, W. Beckert, W.S. Kreher, J. Eur. Ceram. Soc. 24 (2004) 1647.