An evaluation of a class of phenomenological theories of ferroelectricity in polycrystalline ceramics

Various phenomenological theories of ferroelectricity in polycrystalline ceramics have been proposed in recent years. A particularly attractive class of multiaxial theories with a reduced number of internal variables hinges upon an additive decomposition of the strain and the dipole density into reversible internal variables associated with elasticity and dipole perturbations, and irreversible internal variables associated with dipole switching. It has, however, been recently recognized that these theories can provide unexpected predictions for certain—yet unexceptional—loading histories. The source of the problem was pinned down to the nonconvex dependence of the internal energy on the irreversible variables. The purpose of the present study is to evaluate this more thoroughly. It is found that predictions become unstable above a certain level of mechanical stress, which can be on the order of a few megapascals or even lower for typical sets of material parameters employed in the literature. It is argued that this class of theories should be used with caution, even within their presumed range of validity.


Introduction
The increasing need for viable computational tools to assess the performance of electro-deformable devices has motivated the recent development of various phenomenological theories of ferroelectricity in polycrystalline ceramics. A particularly attractive class of multiaxial theories with a reduced number of internal variables hinges upon an additive decomposition of the strain and the electric dipole density into reversible internal variables associated with elasticity and dipole perturbations, and irreversible internal variables associated with dipole alterations or switching [1,2]. Constitutive relations and evolution laws for these variables then follow from specific forms assumed for the internal energy and dissipation of the solid in accordance with the thermodynamic framework of Bassiouny and Maugin [3]. These theories are able to reproduce essential features of ferroelectric material responses such as nonlinear stress-strain curves, hysteresis and butterfly loops, and dipole rotation, and are amenable to numerical implementation into efficient finite-element codes [4][5][6].
More recently, the above class of phenomenological theories was put in incremental form by Miehe and Rosato [4] following ideas originally advocated in the context of viscoplasticity (e.g., [7]). The resulting incremental theories allowed for the recast of boundary-value problems in variational form. This variational character was exploited by Bottero and Idiart [8] to derive homogenization estimates for the piezoelectric properties of composite ferroceramics.
In the course of that study, however, it was recognized that this class of theories can provide unexpected predictions for certain-yet unexceptional-loading histories such as electric cycling under moderate uniaxial traction. The source of the problem was pinned down to the nonconvex dependence of the internal energy on the irreversible dipole density. The purpose of the present study is to evaluate this more thoroughly. A brief overview of the theories is provided in Sect. 2, and sample predictions for typical ferroelectric solids are reported in Sect. 3. Finally, a concluding discussion is provided in Sect. 4.

Précis of the phenomenological theories
The class of phenomenological theories considered in this work is based on an additive decomposition of the infinitesimal strain and the electric dipole density into reversible internal variables associated with elasticity and dipole perturbations, and irreversible internal variables associated with dipole alterations or switching, with no account of thermal and free-charge effects. These theories are relevant to polycrystalline ferroceramics subjected to quasi-static electromechanical loadings. We restrict attention to a subclass of simple theories proposed by McMeeking and Landis [2] and refined by Miehe and Rosato [4]. These theories describe the electromechanical state of the solid by a strain tensor ε, a polarization vector P, and an irreversible polarization vector characterizing the average amount and orientation of the permanent electric dipoles relative to a stress-free configuration with no net polarization. The fact that permanent dipoles are characterized by a single variable simplifies the mathematical formulation considerably but renders the theory unable to account for dipole switching induced by purely mechanical loadings. Predictions are therefore expected to be reasonably accurate when electric fields are strong and stresses are small; notwithstanding this expectation, however, the theory should be robust for the entire range of physically admissible electric field intensities and stresses.

Thermodynamics and constitutive relations
In view of the above basic assumptions, the free energy of a system comprising a ferroelectric solid occupying a domain Ω and an electric field is given by where e is the free energy density of the solid, E is the intensity of the electric field, and 0 is the electric permittivity of vacuum. The free energy of the solid is assumed to be a convex function of ε and P. In turn, the dissipation of the system is entirely ascribed to dipole switching within the solid, and is assumed to be of the form where the dissipation potential ϕ is a convex, positive function of the irreversible polarization rate˙ such that ϕ(0) = 0. The form (2) guarantees positivity of the dissipation as required by the principles of thermodynamics.
Thermodynamic arguments then imply that the constitutive relations of the solid are given by (see [3]) where the first two expressions relate the electric field intensity E and stress tensor σ with the polarization and strain, respectively, and the last expression provides the evolution law for the irreversible polarization. In the case of nonsmooth potentials, the derivatives in (3) should be understood as subdifferentials.
The polarization can be eliminated from the constitutive description in favor of the electric field intensity by introducing the free enthalpy density where the first term corresponds to a partial Legendre transformation of e with respect to P and ε. The constitutive relations (3) can then be rewritten as where D = 0 E + P is the electric displacement within the solid.

Incremental constitutive relations
Miehe and Rosato [4] recasted the above constitutive relations (5) in incremental form by discretizing them in time following an implicit Euler scheme and introducing a single incremental potential. More specifically, given the value i of the irreversible polarization at the time step i of an electromechanical process, the electric displacement and strain at the time step i + 1 are expressed in terms of the stress and the electric field intensity at that time step as where w is an incremental potential defined by In the above expression, Δt = t i+1 − t i is the time interval and the optimal corresponds to the irreversible polarization i+1 at the time step i + 1 that solves the discretized evolution law (5) 3 for the given i . Given a particular choice of thermodynamic functions in (7), expressions (6) provide an incremental phenomenological description of ferroelectricity which allows for the variational formulation of boundary-value problems.

Specific forms of the thermodynamic functions
The theories considered in this work amount to postulating a free energy of the form where the elasticity and electric polarizability tensors are given in terms of isotropic constants by C = λI ⊗ I + 2μI and κ = κI, (9) and the piezoelectricity tensor is given in terms of three piezoelectric moduli by In the above expressions, I and I are the second-and fourth-order identity tensors with major and minor symmetry, and the symbol ⊗ s refers to symmetrization in the first two indices. The first three terms in (8) represent the elastic, electric polarization, and piezoelectric energies, respectively, and define a quadratic form which must satisfy the inequalities to comply with the convexity requirement on the free energy. In turn, the stored energy e pd due to permanent dipoles is given by where p s is the saturation polarization and h s is a material parameter characterizing the electric hysteresis slope, while the straining associated with the permanent dipoles is given bŷ where the symbol ⊗ d denotes the deviatoric part of the tensor product, ε s is the uniaxial strain at saturated polarization, and the exponent γ is an additional material parameter characterizing the form of the butterfly loop such that γ ≥ 1. Finally, the theories postulate a dissipation potential of the form where e c is the coercive field strength of the solid, i.e., the electric field level above which dipole switching is triggered; e 0 andṗ 0 are, respectively, a reference electric field and a polarization rate characterizing the ratedependence of the switching process; and n is a rate-sensitivity exponent such that 0 ≤ n ≤ 1. Finally, the corresponding free enthalpy density is given by witĥ The values roughly reproduce the rate-dependent behavior of a polycrystalline lead zirconate titanate at low frequencies [4,10] where S = C −1 is the compliance tensor of the unpoled solid,κ = κ − h T Sh, and the superscript T in the case of third-order tensors denotes transposition between the first pair of indices and the last index.
With the thermodynamic functions completely specified, predictions for any electromechanical loading history can be readily obtained by successive evaluations of the incremental relations (6). However, it should be emphasized at this point that, while the dissipation potential (14) is a convex function of the irreversible polarization rate, the free energy (8) is a nonconvex function of the irreversible polarization, and so the constitutive theory does not conform to the generalized standard material model [9]. Consequently, stability and uniqueness of the predicted response are not guaranteed. In fact, it will be seen in the following sections that predictions can indeed be unstable, as anticipated by Bottero and Idiart [8].

Sample predictions
The phenomenological theory of Sect. 2 is used here to generate predictions for ferroelectric solids subjected to the simultaneous action of electric fields and mechanical loads. Table 1 shows the numerical values adopted for the various material parameters. These values roughly reproduce the rate-dependent behavior of a polycrystalline lead zirconate titanate at low frequencies [4,10] and comply with conditions (11). The irreversible strain exponent is commonly set in the literature to either γ = 1 or γ = 2; we consider both values.
In all cases, initially unpoled specimens are subjected to a mechanical stress σ along a given direction n, and are subsequently subjected to a triangular electric signal E(t) with a peak amplitude of 4e c and a frequency of 1 Hz along the same direction, with the applied stress level held fixed. For each loading step, the incremental potential (7) is computed by maximizing the function J with a direct search complex algorithm for nonsmooth functions, and the material response is obtained by evaluating expressions (5) 1,2 at the optimal irreversible polarization. The time step employed in the calculations is Δt = 10 −3 s.
Predictions for specimens with γ = 1 and γ = 2 are shown in Figs. 1 and 2, respectively, in the form of plots for the electric displacement (D ) and axial strain (ε ) in the direction of the applied loading as a function of the electric field intensity. Also shown in these figures are corresponding plots for the peak values of the irreversible polarization components along the parallel ( = · n) and perpendicular ( ⊥ = | ⊥ | = | − n|) directions to the γ = 1 a P M 0 5 a P M 0 a P M 0 5 - Fig. 1 Predictions for specimens with irreversible strain exponent γ = 1 at various stress levels: electric displacement (D ) and longitudinal strain (ε ) along the compression axis, and parallel ( ) and perpendicular ( ⊥ ) irreversible polarization components applied loads versus the number of loading cycles. Three different stress levels corresponding to compressed (− 50 MPa), free-standing (0 MPa), and tensioned (50 MPa) specimens are considered 1 . Initial transients and stabilized responses are indicated in dotted and continuous lines, respectively. We begin by noting that the responses of freestanding and tensioned specimens predicted by both irreversible strain exponents are consistent with experimental observations. These predictions reproduce the hysteretic electric response and butterfly loops typically observed in ferroceramics. The exponent γ = 2 is seen to produce somewhat more rounded responses than the exponent γ = 1, but other than that, the responses are qualitatively similar. The underlying irreversible polarization remains aligned with the direction of applied loads, as it should in view of the initial isotropy of the material, and the reponse stabilizes right after the first loading ramp. The predicted responses for the compressed specimens, on the other hand, exhibit unexpected features. Indeed, these responses exhibit longer transients and stabilized irreversible polarizations misaligned with the direction of applied loads. This means specimens develop polarization charges on the free surfaces perpendicular to the applied loads and distort with misaligned principal axes. The source of these features can be pinned down to the convex dependence of the assumed free energy density (8) on the irreversible polarization. This convexity is introduced via the first and third terms of the free energy, which represent the elastic and piezoelectric contributions. The relevant term is the first one. Indeed, neglecting the piezoelectric effect and the contribution of the irreversible polarization to the complianceŜ and permittivityˆ tensors, the objective function (7) 2 being maximized with respect to at each time step is given by where σ is the applied stress level, Y and are the Young's modulus and electric permittivity, respectively, of the solid,˙ = ( − i )/Δt. The evolution equation (5) 3 for the irreversible polarization in the plane perpendicular to the loading axis, ⊥ , is thus for the exponent γ = 1, and for the exponent γ = 2. The primal solution ⊥ (t) = 0 to these equations is stable, provided the bracketed factors multiplying ⊥ on the left-hand sides are positive. While this is the case for any tensile stress level (σ > 0), it is not the case for sufficiently large compressive stresses (σ < 0). In the latter case, an exponential build-up of small perturbances in ⊥ (t) should be expected. This expectation conforms with the results reported in figures 1 and 2 . Note that ⊥ (t) becomes finite right after the first loading cycle in specimens with γ = 1, whereas it does so only after several cycles in specimens with γ = 2. Of course, the number of loading cycles it takes ⊥ (t) to become finite increases with the decreasing stress level. This instability is by no means restricted to uniaxial loadings. Consider, for instance, specimens subject to the same loading program considered above except that the mechanical stress is now equibiaxial on a plane perpendicular to the applied electric field. In view of the isochoric character ofε( ), the evolution equations for the irreversible polarization perpendicular to the direction of the applied electric field are given by the same expressions (18) and (19) but with the opposite sign in front of σ . Thus, the primal solution ⊥ (t) = 0 is now stable for any compressive stress level (σ < 0), but it is unstable for sufficiently large tensile stresses (σ > 0).
Stability limits for arbitrary loading conditions follow from knowledge of the Hessian of the free energy (8), which is quite involved. In any event, the above analysis suggests instabilities should be expected whenever specimens are cycled under electric field intensities exceeding e c along with stress levels exceeding ∼ h s p s 2 /ε s . Such stress levels are in the order of 5 MPa for the specimens considered here. While it is conceivable that features of this sort may arise in small anisotropic monocrystalline samples (see, for instance, [12]), they have never been reported in highly disordered polycrystalline samples, to the best of our knowledge. In fact, recent micromechanical theories for ferroceramics intended to generate physically-based models [13,14] make use of energy densities that are convex with respect to the volume fractions of variants and are therefore unlikely to predict unstable responses of this sort. Thus, even though these phenomenological theories are intended for low stress levels only, for example, McMeeking and Landis [2], they should be used with caution even within their presumed range of validity.

Concluding discussion
The study indicates that the above class of phenomenological theories generate unexpected predictions above certain stress levels, which were found to be in the order of a few megapascals or even lower for the sets of material parameters considered. At such stress levels, the predicted response can become unstable. This is a consequence of the assumed convex dependence of the elastic energy on the permanent dipole density. Even though alternative functional dependences have not been considered in this study, the convexity seems inevitable if the set of internal variables describing the permanent dipole distribution is to be restricted to a single vectorial variable as in this class of theories. This is because the physics of ferroelectricity then require that the strain tensor associated with permanent dipoles be an even function of such vectorial variable, and the elastic energy be a quadratic form of such strain tensor. Thus, the issue is expected to persist in more sophisticated versions of the theories making use of an irreversible strain tensor as an additional internal variable in order to capture ferroelasticity, as in Kamlah [1] and Landis [15], since stress levels typically required to trigger ferroelasticity are significantly higher than those alluded to above. These versions have also been found to suffer from additional problems arising from the assumed form for the dissipation potential [16].
We note in passing that Miehe et al. [17,18] proposed a phenomenological theory of ferromagnetism in polycrystalline ceramics that is completely analogous to the above theory of ferroelectricity, where the electric field intensity, electric polarization, and electric displacement vectors are identified with the magnetic field strength, magnetization, and magnetic induction vectors, respectively. This correspondence is of course valid, provided free charges and internal currents are neglected. The above observations are also valid for this theory in view of the analogous mathematical structure.
It is thus concluded that robust phenomenological theories require the use of enriched sets of internal variables to describe the distribution of permanent dipoles within the solid. Efforts along these lines have been pursued by Kamlah [1], Mehling et al. [19] and, more recently, by Stark et al. [13] and Tan and Kochmann [14], but the resulting theories are yet to be thoroughly evaluated.