The Homogenized Energy Model (HEM) for Characterizing Polarization and Strains in Hysteretic Ferroelectric Materials: Implementation Algorithms and Data-Driven Parameter Estimation Techniques Zhengzheng Hu and Ralph C. Smith Department of Mathematics Center for Research in Scientific Computation North Carolina State University Raleigh, NC 27695 Jon Ernstberger Department of Mathematics LaGrange College LaGrange, GA 30240 Abstract Ferroelectric materials, such as PZT, PLZT, PMN and BaTiO , provide unique actuator and 3 sensor capabilities for applications including nanopositioning, high speed valves and fuel injectors, camera focusing and shutter mechanisms, ultrasonic devices for biomedical imaging and treatment, and energy harvesting devices. However, to achieve the full potential of the materials, it is necessary to develop and employ models that quantify the creep, rate-dependent hysteresis, and constitutive nonlinearities that are intrinsic to the materials due to their domain structure. The success of mod- els requires that they be highly efficient to implement since real-time applications can require kHz to MHz rates. The calibration of models for specific materials, devices, and applications, requires efficient and robust parameter estimation algorithms. Finally, control designs can be facilitated by models that admit efficient and robust approximate inversion. The homogenized energy model (HEM) is a multiscale, micromechanical framework that quantifies a range of hysteretic phenomena intrinsic to ferroelectric, ferromagnetic and ferroelastic materials. In this paper, we present highly efficient implementation and parameter estimation algorithms for the ferroelectric model. This in- cludes techniques to construct analytic Jacobians and data-driven algorithms to determine initial parameter estimates to facilitate subsequent optimization. The efficiency of these algorithms fa- cilitates material and device characterization and provides the basis for constructing efficient and robust inverse algorithms for model-based control design. The model implementation, calibration, and validation are illustrated using rate-dependent PZT data and single crystal BaTiO data. 3 i Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE 3. DATES COVERED JAN 2012 2. REPORT TYPE 00-00-2012 to 00-00-2012 4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER The Homogenized Energy Model (HEM) for Characterizing Polarization 5b. GRANT NUMBER and Strains in Hysteretic Ferroelectric Materials: Implementation Algorithms and Data-Driven Parameter Estimation Techniques 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION North Carolina State University,Center for Research in Scientific REPORT NUMBER CRSC-TR12-02 Computation,Department of Mathematics,Raleigh,NC,27695-8212 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release; distribution unlimited 13. SUPPLEMENTARY NOTES 14. ABSTRACT Ferroelectric materials, such as PZT, PLZT, PMN and BaTiO3, provide unique actuator and sensor capabilities for applications including nanopositioning, high speed valves and fuel injectors camera focusing and shutter mechanisms, ultrasonic devices for biomedical imaging and treatment and energy harvesting devices. However, to achieve the full potential of the materials, it is necessary to develop and employ models that quantify the creep, rate-dependent hysteresis, and constitutive nonlinearities that are intrinsic to the materials due to their domain structure. The success of models requires that they be highly efficient to implement since real-time applications can require kHz to MHz rates. The calibration of models for specific materials, devices, and applications, requires efficient and robust parameter estimation algorithms. Finally, control designs can be facilitated by models that admit efficient and robust approximate inversion. The homogenized energy model (HEM) is a multiscale, micromechanical framework that quantifies a range of hysteretic phenomena intrinsic to ferroelectric, ferromagnetic and ferroelastic materials. In this paper, we present highly efficient implementation and parameter estimation algorithms for the ferroelectric model. This includes techniques to construct analytic Jacobians and data-driven algorithms to determine initial parameter estimates to facilitate subsequent optimization. The efficiency of these algorithms facilitates material and device characterization and provides the basis for constructing efficient and robust inverse algorithms for model-based control design. The model implementation, calibration and validation are illustrated using rate-dependent PZT data and single crystal BaTiO3 data. 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF ABSTRACT OF PAGES RESPONSIBLE PERSON a. REPORT b. ABSTRACT c. THIS PAGE Same as 35 unclassified unclassified unclassified Report (SAR) Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 Nomenclature d Piezoelectric constant (m/V = C/N) h Piezoelectric constant (V/m = N/C) E Electric field (V/m) g Gibbs energy for dipoles (CV) G Gibbs energy density of α-variant (CV/m3) α Ge Electric Gibbs energy density of α-variant (CV/m3) α P Polarization (C/m2) P Polarization kernel (C/m2) Pα Polarization of α-variant (C/m2) Pα Remanent polarization of α-variant (C/m2) R Pα Minimum polarization of α-variant (C/m2) m sE Elastic compliance of α-variant at constant field (m2/N) α T Temperature (K) x ,x ,x Fraction of positively, negatively and 90o oriented dipoles (Unitless) + − 90 YP Elastic stiffness of α-variant at constant polarization (N/m2) α (cid:15) Permittivity (F/m = C/Vm) ε Strain (Unitless) ε Strain kernel εα Strain of α-variant εα Remanence strain of α-variant R εα Minimum strain of α-variant m ηε Inverse susceptibility at constant strain (m/F = Vm/C) σ Stress (N/m2) τ Relaxation time for 90o switching (s) 90 τ Relaxation time for 180o switching (s) 180 χ Electric susceptibility (Unitless) e χσ Ferroelectric susceptibility of α-variant at constant stress (F/m = C/Vm) α ψ Helmholtz energy density of α-variant (CV/m3) α 1 Introduction Ferroelectric materials, such as lead zirconate titanate (PZT), lanthanum-doped lead zirconate ti- tanate (PLZT), lead manganese niobate (PMN), polyvinylidene flouride (PVDF), and barium ti- tanate (BaTiO ), are being widely considered as actuators, sensors, and structural units for a range 3 ofaerospace,aeronautic,automotive,industrialandbiomedicalapplications. Theiradvantagesderive from the intrinsic properties of the materials. The complementary direct and converse piezoelectric effectsimbuethematerialswithmultipledesignpropertiesincludingsensor,actuator,self-monitoring and nondestructive evaluation, and energy harvesting capabilities. Their functionality is augmented by their solid state nature which promotes miniaturization and simplified design and reduces power requirements and heat generated by the units. Furthermore, ferroelectric actuators can provide nanometer positioning resolution and can operate at frequencies ranging from DC to MHz. As de- tailed in the companion paper [19], ferroelectric transducers are being considered, or are already being employed, for a large number of applications including high speed valves for fuel injectors, fer- roelectric memory technologies (e.g., FeRAM), nanopositioning units (e.g., AFM and STM stages), 1 high speed camera shutters and autofocusing units, ultrasonic transducers for biomedical imaging and treatment, micro and pico air vehicle design, and energy harvesting devices. However, the ferroelectric properties that imbue the materials with unique transducer capabil- ities also produce rate-dependent hysteresis, creep, and constitutive nonlinearities as illustrated in Figure 1. The field-polarization and field-strain data from [26] illustrates that rate-dependent effects aresignificantatfrequenciesas low as1 Hz. PZTdatafrom [24]illustratesthatforfixedfieldinputs, both the strain and polarization exhibit significant creep on timescales of 1 to 20 seconds. Finally, the MFC data from [8] illustrates nested minor loop behavior typical of moderate drive regimes. Whereas the full switching behavior shown in Figure 1(a)-(e) will typically not be encountered in applications, general models must account for the full range of rate-dependent and creep behavior to provide comprehensive device characterization. To be optimally utilized in applications, constitutive models must satisfy the following criteria: (i) they must adequately quantify the range of rate-dependent behavior exhibited by the materials, (ii) they must be in a form that can be readily extended to provided distributed models, for complex structuresanddevices,thatareamenabletofiniteelementimplementation,(iii)theymustbeefficient to implement, (iv) they must be readily calibrated for a specific material or device, and (v) models employed in control designs can prove advantageous if they admit the construction of approximate inverse relations. It is shown in [18–21] that the homogenized energy model (HEM) satisfies (i) and (ii). Fur- thermore, the construction of approximate inverses for previous formulations of the framework are reported in [4,18] and inversion techniques for the present formulation [19] are addressed in [15]. In this paper, we address (iii) and (iv) by presenting highly efficient implementation techniques and data-drivenparameterestimationalgorithmsforthehomogenizedenergymodel. Incombination,this 0.2 0.2 0.% 0.15 0.15 %m1 0.1 3)45i64ti38 ,9/!0.01 Strain (%)0.00.51 Strain (%) 0.00.51 2 0 !0.% 0 !0.05 !1.5 !1 !0.5 0 0.5 1 1.5 !1.5 !1 !0.5 0 0.5 1 1.5 0 100 200 300 400 500 &ie)* ,-./m1 Field (-V/m) Time (s) (a) (b) (c) * $!!) ! !!$!!37> 0.% 4.5 !!%!!37> !% !!)!!37> 0.& 4 !!’!!37> !) 24lari7ati4n (:/m)!000...0121 Strain (%)123...23555 7-89/:;.5.<= 156!!!!$$%!’( 3!0.2 1 0.01 Hz !$) 0.01 ;7 !0.& 0.1 ;7 0.5 0.1 Hz !$( 1 ;7 1 Hz 1?6 !0.%!2 !1 0 1 2 0!2 !1 0 1 2 !$’! !"# $ $"# % %"# & Field (.V/mm) Field (kV/mm) ,-./d 1M3456 (d) (e) (f) Figure 1: (a) Field-polarization, (b) field-strain and (c) time-strain PZT data from [24]. Rate- dependent (d) polarization and (e) strain PZT data from [26]. (f) Field-strain MFC data from [8]. 2 provides the framework with significant flexibility for use in simulation packages, design algorithms, and control designs for systems utilizing ferroelectric transducers. To provide context for the framework, we first summarize the mechanisms that produce the hys- teretic behavior shown in Figure 1. The creepand rate-dependencies are due to the fact that kinetics associated with dipole switching typically differ from mechanical or electrical loading rates. Further- more, it is discussed in [19] that the 180o switching associated with large polarization changes often occurs at a different rate than 90o switching associated with large strain changes. In combination, this establishes that multiple time scales must be incorporated in dynamic models. As illustrated for a PZT-based macro-fiber composite (MFC) in Figure 2, hysteresis also involves multiple spatial scales. Within a domain, strains or changes in polarization are due to stress-induced material deformations or field-induced ion movement and the behavior is often reversible and linear. Field or stress-induced dipole switching at the grain level produces irreversible hysteresis in both the field-polarization and field-strain relations. For single crystal materials comprised of a single grain, the switching is typically rapid thus producing sharp hysteresis and butterfly loops in quasistatic operating regimes. For polycrystalline materials with distributed interaction and coercive fields, hysteresis and butterfly loops are smoothed due to nonuniform grain contributions. The behavior of hysteresis loops is further modified when hysteretic actuator or sensor materials are employed on distributed structures. Figure 2: Multiscale behavior of a PZT-based MFC transducer at the application, device, material and unit cell levels. 3 As noted in Figure 2, material-level constitutive relations for the polarization P and strain ε can generically be expressed as P(E,σ) = d(E,σ)σ+χσE +P (E,σ) irr ε(E,σ) = sEσ+d(E,σ)E +ε (E,σ) irr where E,σ are input fields and stresses, χσ is the ferroelectric susceptibility at constant stress, sE is the elastic compliance at constant field, and d is a piezoelectric coupling coefficient. We note that d(E,σ),P (E,σ) and ε (E,σ) incorporate the nonlinear and irreversible history-dependence due irr irr to dipole switching. Various modeling hierarchies can be defined by the manner in which d,P and ε are con- irr irr structed. Micromechanical, or microscopically-motivated models are based on an energy description ofthematerialatthedomainorgrainlevelincombinationwithvarioushomogenizationtechniquesto provide expressions for the nonlinear, effective components d,P and ε . Phenomenological mod- irr irr elscircumventthedifficultiesassociatedwithquantifyingcomplex, orpoorlyunderstood, micro-scale physics by constructing relations for d,P and ε based on macroscale observations or experimen- irr irr tal measurements, often guided by thermodynamic constraints. Details regarding various models for ferroelectric materials can be found in [2,9,13,16,18,19]. Thehomogenizedenergymodelisamultiscale,micromechanicalapproachthatbeginswithenergy analysis at the domain level to construct the linear relations shown in Figure 2 for Pα and εα where α designates the dipole variant; e.g., α = ±180,90 for tetragonal materials. Switching processes due to domain wall nucleation and movement are incorporated by tracking the evolution of dipole fractions x which serve as internal variables. Differential equations quantifying the dynamics of x α α are driven by likelihood rates constructed using Boltzmann theory to quantify the scaled probability of transitioning between stable equilibria associated with dipole variants. This incorporates the rate-dependence and multiple-time scales exhibited by the data in Figure 1 and yields grain-level or single crystal kernels P and ε. For polycrystalline materials, material-level models are constructed by assuming that properties such as coercive fields, critical driving forces, and interaction fields are manifestations of underlying densities rather than constants. The three primary issues that must be addressed to construct implementation algorithms are: (i) efficient approximation of the integrals, posed on infinite and semi-infinite domains, arising in the homogenizationstepusedtoconstructtheirreversiblecomponentsd(E,σ),P (E,σ)andε (E,σ), irr irr (ii) efficient approximation of the evolution equations for x , and (iii) efficient evaluation of the like- α lihoods which incorporate rate-dependent effects and drive the evolution equations. The first issue is addressed by employing density representations that admit efficient numerical approximation on finite domains whereas implicit Euler discretizations readily address the second issue. The repeated evaluationoflikelihoodsisavoidedbyconstructingarraysthatpermithigh-speedaccessduringmodel implementation. In combination, the evaluation of d(E,σ),P (E,σ) and ε (E,σ) are reduced to irr irr operations solely involving componentwise matrix multiplication and summation. Hence the oper- ations are highly efficient and inherently parallelizable which facilitates subsequent implementation using devices such as field programmable gate arrays (FPGA) [3]. We note that approximate model inversion techniques rely on forward solutions of the model. Hencetheefficiencyofforwardalgorithmsalsodirectlycontributestothespeedofinversealgorithms. Therearetwoprimaryissuesthatmustbeaddressedwhenemployinggradient-basedoptimization algorithmsfortheparameterestimationrequiredformodelcalibration: (i)determinationofaccurate initialparameterestimatesand(ii)efficientandaccurateconstructionofgradientsorJacobians. Due to its physical nature, parameters in the homogenized energy model can be directly correlated with properties of measured data. We utilize that here to construct data-driven algorithms to determine 4 initial parameter estimates. Secondly, we construct analytic Jacobian relations which significantly improves both the speed and accuracy of optimization routines. In combination, this helps mitigate some of the difficulties associated with the shallow slopes and nonconvexity (multiple local minima) inherent to functionals employed in parameter estimation problems of this type. We note that various aspects of the implementation and data-driven parameter estimation algo- rithms can be optimized to provide improved performance for specific data sets. Instead, our goal was to provide algorithms that balance simplicity, robustness, accuracy, and optimality for a wide range of materials, applications and operating regimes. Points were algorithms can be modified or optimized are noted in the discussion. Following a brief summary of the model in Section 2, highly efficient implementation algorithms are presented in Section 3. The parameter estimation problem is discussed in Section 4 where we provide algorithms to construct analytic Jacobians and data-driven techniques to determine initial parameter estimates. Model calibration and validation are illustrated in Section 5 using PZT data from [24] and single crystal BaTiO data from [5]. Additional examples illustrating the capability of 3 themodeltocharacterizecreepandrate-dependentdata, suchastheshowninFigure1, areprovided in the companion paper [19] which details the model development. 2 Polarization and Strain Models As detailed in [19], consideration of 180o dipole switching yields a macroscopic model quantifying the nonlinear and hysteretic field-polarization map for ferroelectric materials whereas the additional incorporation of 90o switching yields a homogenized model quantifying both the polarization and strains due to input fields and stresses. We summarize the latter model and indicate simplifications that occur if considering only the polarization in the absence of applied stresses. Polarization and Strain Model At the lattice level, we consider the Helmholtz and Gibbs energy densities 1 1 ψ (P,ε) = ηε(P −Pα)2+ YP(ε−εα)2+h (P −Pα)(ε−εα) α 2 α R 2 α R α R R and G (E,σ;P,ε) = ψ (P,ε)−EP −σε α α where we indicate ±180o and 90o orientations by α = ±,90. As summarized in the nomenclature table at the beginning of the paper, Pα,εα,ηε,YP and h denote the remanence polarization, rema- R R α α α nence strain, inverse susceptibility at constant strain, elastic stiffness at constant polarization, and piezoelectric constant. For a fixed applied field and stress, the conditions ∂G = 0 and ∂G = 0 yield the relations ∂P ∂ε Pα = Pα+χσE +d σ m R α α (1) εα = εα +d E +sEσ m R α α where YP h ηε χσ = α , d = α , sE = α α YPηε −h2 α h2 −YPηε α YPηε −h2 α α α α α α α α α aretheferroelectricsusceptibilityatconstantstress,thepiezoelectricconstant,andelasticcompliance at constant field. The minimum of the Gibbs energy in each α-well can then be expressed as 1 1 G (E,σ) = − χσE2− sEσ2−d Eσ−EPα−σεα. αm 2 α 2 α α R R 5 The dipole fractions x ,x and x associated with positively, negatively, and 90o degree dipoles + − 90 evolve according to the differential equations x˙ = −(p +p )x +p x +p x − −90 −+ − 90− 90 +− + x˙ = p x −(p +p )x +p x (2) 90 −90 − 90− 90+ 90 +90 + x˙ = p x +p x −(p +p )x + −+ − 90+ 90 +90 +− + where 1 pαβ(E,σ) = e−∆Gaαβ(E,σ)V/kT (3) τ αβ quantifies the likelihood of transitioning from an α-well to a β-well. As noted in Table 1, (2) can be simplified using the relation x +x +x = 1. The activation energy is specified by the relation + − 90 ( ∆G (1−F (E,σ)/F )2 , F (E,σ) ≤ F ∆Ga (E,σ;F ) = 0 αβ c αβ c αβ c 0 , F (E,σ) > F . αβ c Here F (E,σ) = G (E,σ)−G (E,σ) is the thermodynamic driving force. The specific form of αβ αm βm F , based on the physical assumption that αβ χσ = χσ = χσ = χσ + − 90 sE = sE = sE = sE + − 90 P90 = 0,P+ = −P−,ε+ = ε− (4) R R R R R d = 0,d = −d 90 − + τ = τ = τ = τ = τ , τ = τ = τ , 90− −90 90+ +90 90 +− −+ 180 is given in Table 1. Note that ( 1F , 180o Switching ∆G = 4 c 0 1 F , 90o Switching 16 c is the energy barrier at zero driving force. The polarization and strain kernels are given by X X P = x Pα , ε = x εα. α m α m α=±,90 α=±,90 Based on (1) and (4), these relations can be expressed as P(E,σ) = d¯(E,σ)σ+χσE +P (E,σ) irr ε(E,σ) = sEσ+d¯(E,σ)E +ε (E,σ) irr where X d¯(E,σ) = d x (E,σ) α α α=±,90 X P (E,σ) = Pαx (E,σ) irr R α α=±,90 X ε (E,σ) = εαx (E,σ). irr R α α=±,90 6 In the final step in the development of [19], macroscopic models Z ∞Z ∞ P(E(t),σ(t);x0) = P(E(t)+E ,σ(t);F )ν (E )ν (F )dE dF + I c I I c c I c 0 −∞ Z ∞Z ∞ ε(E(t),σ(t);x0) = ε(E(t)+E ,σ(t);F )ν (E )ν (F )dE dF + I c I I c c I c 0 −∞ are constructed by considering interaction fields E and thermodynamic driving forces F to be I c manifestations of underlying densities ν and ν . The final models can thus be expressed as I c P(E,σ) = d(E,σ)σ+χσE +P (E,σ) irr ε(E,σ) = sEσ+d(E,σ)E +ε (E,σ) irr where Z ∞Z ∞ d(E,σ) = d¯(E ;F )ν (E )ν (F )dE dF e c I I c c I c 0 −∞ Z ∞Z ∞ P (E,σ) = P (E ;F )ν (E )ν (F )dE dF irr irr e c I I c c I c 0 −∞ Z ∞Z ∞ ε (E,σ) = ε (E ;F )ν (E )ν (F )dE dF . irr irr e c I I c c I c 0 −∞ Here E (t) = E(t)+E is the effective field. e I Details regarding various choices for the densities are provided in [19]. For the implementation and parameter estimation algorithms detailed in this paper, we employ the representations νc(Fc) = P1‘α‘ kXK=kαααkφk(Fc) , φk(Fc) = σckFc1√2πe−[ln(Fc)−ln(F¯c)]2/2(σck)2 , σck = 2kσc (5) K 1 Xβ νI(EI) = P β βkϕk(EI) , ϕk(EI) = σk√12πe−EI2/2(σIk)2 , σIk = 2kσI ‘ ‘ k=k I β where the preceding sums ensure integration to unity. During model calibration, the parameters {α ,β } are determined through a least squares fit to the data. k k The complete polarization and strain model is summarized in Table 1. Parameters for the Polarization-Strain Model To construct the density basis functions ϕ (E ) and φ (F ), it is necessary to determine values k I k c for the driving force mean µ = ln(F¯ ) and variance σ2 as well as the interaction field variance σ2. c c c I We denote this set of parameters by p¯= [µ ,σ2,σ2]. (6) c c I We provide an algorithm in Section 4.2 to estimate values for these parameters based on measured attributes of E-P data. We note that these parameters are not updated or optimized through a least squares fit to data. Based on the assumption (4), the remaining parameters are denoted by p = [P+,ε+,ε90,χσ,d ,sE,γ,τ ,τ ,α ,β ]. (7) R R R + 90 180 k k This comprises the set that is optimized during model calibration. 7