DTIC ADA556961: The Homogenized Energy Model (HEM) for Characterizing Polarization and Strains in Hysteretic Ferroelectric Materials: Implementation Algorithms and Data-Driven Parameter Estimation Techniques

1.6 MB·English

by  Defense Technical Information Center

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