Fatigue comprehensive modeling for ferro- and piezo-electric effect possessed in smart materials

Abstract

Ferroelectric materials have been used in numerous applications as ultrasonic and sonar transducers etc. due to their distinct ferroelectric and piezoelectric properties. Their characteristic loops measured from experiment demonstrate generally nonlinear hysteresis, and the hysteresis is also easily to be deteriorated to offset or asymmetry by some inevitable reasons like electric fatigue. The aforementioned complexities aggravate the modeling difficulty and no doubt impede the utilization potentiality of ferroelectric materials. In view of this, the paper is set out to present an appropriate model that can comprehensively describe a wide class of characteristic loops in consideration the deteriorated circumstance aroused by fatigue factor. Finally experimental data is utilized here for validating the superiority of the proposed method to conventional one.

Keywords

Ferroelectric materials Hysteresis loops Butterfly loops Electric fatigue Comprehensive modeling

1. Introduction

Some smart materials like ferroelectric materials usually possess several field-coupling properties as ferro- and piezo-electricity [1–3], which are usually expressed by electric-field-induced polarization (P–E) or electric-field-induced displacement (D–E) and electric-field-induced strain (ϵ–E) loops. Under harmonic field, these measured loops demonstrate representative nonlinear-hysteretic characteristics that should be theoretical symmetric if domain switching and reverse process was idealized. However, ferroelectric materials are likely to lose partial its ability of switchable polarization resulting in their deterioration. This phenomenon could be mainly attributed to electric fatigue, which is a primary obstacle to the full commercialization of ferroelectric materials and their devices.

Up to date, many proposed models have been used to describe these ferro- and piezo-electric behaviors. One category [4] of the models is obtained by microcosmic perspective that are suitable exploring fundamental material behavior rather than for device design. Another category [4] is macroscopic models that focus on phenomenological descriptions of material behaviors. The latter are arguably the better tools for designing devices featuring ferroelectric materials. But both them still need to be improved because they are not flexible enough to suit for the asymmetric circumstance. With regard to the asymmetric modeling for fatigue, it has remained elusive until now. Some fatigue models have been developed through analyzing the dependence of polarization deterioration on the number of electrical cycles. They can be used to theoretically describe the interfered P–E (D–E) and ϵ–E loops by electric fatigue. Herein it is worth mentioning that most relevant models are inevitably involving sets of unknown or inexplicit coefficients. To evaluate these coefficients is very pivotal before they are substituted into the model, which directly determines whether the model is effective and accurate. So far these coefficients are evaluated by experimental measurement or even just empirical means. These test results are, in realities, often susceptible to some uncertain factors aroused by restricted experimental conditions [5]. Thus models substituted by these evaluations may not be effective for describing the specific characteristic behavior in ferroelectrics.

For such cases, it is especially significant to have a fatigue comprehensive model and calibration methods that can rightly reproduce those ferro- and piezo-electric responses. The related work is carried out as follows.

2. Comprehensive model considering fatigue

2.1. Comprehensive equations

Being different from the constitutive relation of piezoelectric materials, ferroelectric materials possess domain switching to arouse two more extra variables as remanent strain ${\varepsilon}_{ij}^{r}$ and remanent polarization $P_{i}^{r}$ [6]. Hence, the linear relationships fitting ferroelectrics are $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle D_{i}-P_{i}^{r}=d_{ikl}{\sigma}_{kl}+{\kappa}_{ij}E_{j}\end{eqnarray}$ (1) $\begin{eqnarray}\displaystyle & \text{} & \displaystyle \text{}{\varepsilon}_{ij}-{\varepsilon}_{ij}^{r}=s_{ijkl}{\sigma}_{kl}+d_{kij}E_{k}\end{eqnarray}$ (2) where D _i and ϵ_ij are electric displacement and total strain produced by the applied stress field σ_kl and electric field E, d _ikl is a tensor of piezoelectric coefficients, 𝜅_ij, is a dielectric permittivity tensor, and s _ijkl is the elastic compliance.

The remanent strain ${\varepsilon}_{ij}^{r}$ , in the light of Landau-Devonshire theory [6], is approximately linked to the remanent polarization through a quadratic relationship between them such that $\begin{eqnarray}\displaystyle {\varepsilon}_{ij}^{r}=\frac{{\varepsilon}_{0}}{2P_{0}^{2}}\left(3P_{i}^{r}P_{j}^{r}-{\delta}_{ij}P_{k}^{r}P_{k}^{r}\right) & & \displaystyle\end{eqnarray}$ (3) where ϵ₀ is the saturated strain corresponding to the maximum achievable polarization P ₀, and plus δ_ij is the Kronecker delta.

As soon as substituting Eq. (3) into Eq. (2), we found the output variables (D _i and ϵ_ij) are basically dependent on the intermediate variable $P_{i}^{r}$ . As such, it is pivotal to establish the relation of polarization with electric and stress fields. In this paper, we adopted one macroscopic model proposed by Miller and McWhorter [7] that are narrated as follows.

Firstly, a hyperbolic tangent function is utilized to fit this saturated polarization $P_{i}^{sat}$ given by $\begin{eqnarray}\displaystyle P_{i}^{sat}={\xi}P_{s}\tanh \left[\frac{{\xi}E_{j}+E_{c}}{2{\delta}}\right]\quad \text{where }{\delta}=E_{c}\left[\ln \left(\frac{1+P_{r}/P_{s}}{1-P_{r}/P_{s}}\right)\right]^{-1} & & \displaystyle\end{eqnarray}$ (4) where P _s, P _r and E _c are three constants whose names are separately saturated polarization, remnant polarization and coercive field. Here, 𝜉 = 1 for an increasing electric field and 𝜉 = −1 for a reducing field.

In real experiment, the behavior of circuits experiencing arbitrary initial conditions cannot be precisely predicted due to incomplete switching. Thus, minor (un-saturated) polarizations would be witnessed more commonly, which are basically surrounded within the saturated loop. On that basis, it is assumed that the derivative of $P_{i}^{s}$ is always smaller than that of the minor one, $P_{i}^{r}$ . Their relationship is thus expressed by the following equation $\begin{eqnarray}\displaystyle \frac{dP_{i}^{r}}{dE_{j}}={\Omega}\frac{dP_{i}^{s}}{dE_{j}}{\Omega}=1-\tanh \left[\left(\frac{P_{i}^{r}-P_{i}^{s}}{{\xi}P_{s}-P_{i}^{r}}\right)^{1/2}\right] & & \displaystyle\end{eqnarray}$ (5) where 𝛺 is an intermediate function whose value ranges from 0 to 1. This equation, together with an initial condition for the electric field, completely defines the hysteretic relationship between polarization and the electric field.

2.2. Fatigue comprehensive equations

As mentioned before, the induced polarization has been proved demonstrating deterioration of hysteresis responses by some defects such as point defects and space charges [8]. On account of structural damage caused by point defect, Yoo found the reduction rate of actual polarization is proportional to the original polarization along with the increasing of the cycle number of fluctuant electric field. On the other hand, Verdier found the change of center of the polarization hysteresis response shows an approximately linear dependence on the logarithm of the electric cycle number. Based on the relation of fatigue polarization with cycle number of electric field, the current (fatigue) polarization $\bar{P}_{i}^{r}$ can be composed of $\begin{eqnarray}\displaystyle \overline{P}_{i}^{r}=P_{i}^{r}(An+1)^{-{\alpha}}+cq\lg (n+1) & & \displaystyle\end{eqnarray}$ (6) where $P_{i}^{r}$ is the original (non-fatigue) polarization before cyclic loading, n is the cycle number, fatigue parameters A, 𝛼 and c express respectively piling, decay and electrode effect in electric fatigue, q is dependent on the charge status where q = 1 for positive charge, q = −1 for negative charge and q = 0 for regions without space charge. Plugging Eq. (4) into Eqs (1)–(5), the fatigued electric displacement D _i and strain ϵ _ij will be obtained eventually. Since now, we gained a comprehensive model for describing the ferro- and piezo-electric behavior of smart materials by consideration of fatigue factor.

2.3. Parameter analysis

Through observation, this model inevitably involves multiple parameters that affect the hysteretic response by means of different shape loops under specific condition (amplitude, frequency or cyclic number of electric field). It is important to delineate their influences on the overall behavior, and, if possible, reveal their physical/macroscopic meaning in case they do not initially admit a direct interpretation. To this end, we carry out a numerical experiment here, and set the material properties to P _s = 0.3 C/m², P _r = 0.25 C/m², E _c = 0.36 MV/m, 𝜉_r = 5000, d ₃₃₃ = 10⁻³ m/MV, ϵ₀ = 0.001, A = 10⁻⁵, 𝛼 = 0.05, c = 0.005. We then apply a harmonic electric field to the material with amplitude 1 MV/m and frequency 0.2 Hz, by substituting it into comprehensive model. We then integrate these coupled differential equations using a standard fourth-order Runge–Kutta method, and obtain electric-field-induced electric displacement and induced strain loops subsequently. The following observations can be made from the obtained results:

First of all, parameters P _s, P _r and E _c are defined in the saturated dipole polarization curve P ^sat, which are all taken as positive quantities, and named as saturated polarization, remanent polarization and coercive field. Parameter ϵ_r is a dielectric constant that affects the slope of the linear fit of D–E loop correspondent to P ^r–E loop in large-field range.

Fig. 1.

Parameters analysis in the comprehensive model.

Without considering fatigue (n = 0), we found piezoelectric coefficient d ₃₃₃ is proportional to the slope of the ϵ–E loop at zero field point (see Fig. 1(a)). Likewise, saturated strain ϵ₀ is proportional to the average value of the ϵ₃₃ at the whole range. Parameters d ₃₃₃ and ϵ₀, generally speaking, reflect piezoelectric capability of the materials. With considering fatigue (n = 10⁶) the primary indication is the decrease of P _r accompanied by a change in E _c, which directly induces the deformation of ϵ–E loop. A decrease in the electric–induced strain is witnessed to make the loop more flat (by comparison of Fig. 1(a) with Fig. 1(b)). Then parameters A and 𝛼 demonstrate the inactive polarization caused by point defect, which make the electric-induced strain loop more compressed. Furthermore as parameter c increases (while space charge positive i.e., q = 1), it is seen that ϵ₃₃ becomes asymmetric with a decreasing of strain range. The left wing of the butterfly loop degrades more rapidly than the right as result of electric fatigue. Meanwhile, the loop center exhibits a shift to the left. These phenomena have been observed in experiments [9]. Based on the preceding description about parameters in the model, they are organized and explained altogether in Table 1.

Table 1

Characteristic parameters in fatigue comprehensive model

Properties	Parameters	Formulation
P _s		Saturated polarization
P _r	Ferroelectric	Remnant polarization
E _c		Coercive field
𝜉_r	Dielectric	Large-field relative permittivity
d	Piezoelectric	Piezoelectric coefficient
ϵ₀		Remanent strain
_A		Piling coefficient
𝛼	Fatigue	Decay coefficient
c		Electrode coefficient

3. Parameter identification by using GA

As shown in Table 1, characteristic parameters have different certain physical meaning, whose evaluations would determine the effectiveness of the model. Thus, we design a GA procedure to optimize the values of these parameters so as to enhance and guarantee the real practicability of the model. It is introduced as follows.

GAs are heuristic optimization methods that imitate the process of natural selection. A generic GA consists of several basic working parts including chromosome representation, fitness function, genetic operators, and termination criteria. Based on the problem we will resolve, the GA is designed as such way:

(1) Chromosome structure

There are nine variables to be identified in Table 1. Thus the chromosome structure becomes $\begin{eqnarray}\displaystyle {\Pi}=\left\{\left[P_{s},P_{r},E_{c},{\xi}_{r},d_{333},{\varepsilon}_{0},A,{\alpha},c\right]_{i}\right\}\quad i=1,\ldots ,N & & \displaystyle\end{eqnarray}$ (7) where N is maximum number of chromosomes, i is the ith individual in chromosome.

(2) Fitness function

A fitness function measures the status of one chromosome achieving the pre-set objective. Each chromosome needs to be awarded an evaluation that indicates how close it came to meeting the objective. Here, we define the fitness function (“Fit”) by using the Root-Mean-Square Error (RMSE) between the simulated and experimental values, which are: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \text{error}_{1}=\sqrt{\frac{1}{M}\left[\mathop{\sum }_{j=1}^{M}D_{sim,j}-D_{\exp ,j}\right]^{2}}\left/\sqrt{\frac{1}{M}\mathop{\sum }_{j=1}^{M}\left[D_{\exp ,j}\right]^{2}}\right.\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \text{error}_{2}=\sqrt{\frac{1}{M}\left[\mathop{\sum }_{j=1}^{M}{\varepsilon}_{sim,j}-{\varepsilon}_{\exp ,j}\right]^{2}}\left/\sqrt{\frac{1}{M}\mathop{\sum }_{j=1}^{M}\left[{\varepsilon}_{\exp ,j}\right]^{2}}\right.\end{eqnarray}$ (9) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \text{Fit}=\max [\text{error}_{1,\min },\text{error}_{2,\min }]\times 100\%\end{eqnarray}$ (10) where error₁ and error₂ severally denote the displacement and strain discrepancies between simulated and experimental data. Each j-th simulation value and its corresponding experimental data are labeled by the subscripts “sim” and “exp”, and M is the total number of data points. Next we choose the larger one between minimal error₁ and minimal error₂ as the fitness function. Apparently less value of “Fit” indicates the better individual among the population until the best (lowest) individual is figure out.

In addition to the above settings, the rest of the main parameters in GA are also designed as: roulette wheel selection, crossover rate is 0.85 and mutation rate is 0.01, and maximum number of iteration is 50. As soon as the GA procedure is used in the model, the whole fatigue-comprehensive-model method is formed consequently and its flowchart is illustrated in Fig. 2.

Fig. 2.

Flow chart of proposed model method by using GA parameter identification.

Next through numerical simulation for experimental data by using such GA method, those parameters will be optimized or assigned certain values according to specific condition including amplitude, frequency and cyclic number of electric field. Finally, we will use curve-fitting method to establish the relationship of parameters with the aforementioned condition, which is the follow-up study. On the basis, the comprehensive model with identified parameters can be applicable in prediction of ferroelectrics’ characteristic in real engineering.

4. Model validation by experimental data

Nuffer [10] has carried out an experiment for investigating fatigue behavior of commercial bulk Lead Zirconate Titanate (PZT) under bipolar cycling. Polarization and strain hysteresis loops were all monitored as well as acoustic emissions (AE). In the experiment, higher cycling field is found yielding stronger fatigue, more incomplete polarization and asymmetric induced strain hysteresis. Here, we use these experimental data to gauge the ability of the proposed model and its GA-based identification procedure in mimicking the D–E and S–E loops, and in extracting the material properties.

As shown in Fig. 3, partial experiment results of strain butterfly loops with solid dots are simulated by the proposed method for three cases 0, 3 × 10⁶ and 10⁸ cycles field whose highest electric field is set as 1.96 kV/mm the optimized parameters and their simulation accuracies (Fit) for these three cases are listed in Table 2. With regard to the case without fatigue, the GA optimized every parameter in the model. Fatigue parameters (A, 𝛼 and c) are fixed by zero in order to optimize the rest parameter more efficiently. On the contrary, the case considering fatigue (i.e. 3 ⋅ 10⁶ and 10⁸ cycles field), the parameters are identified by different certain values, and the fatigue parameters are not zero anymore.

From practicability perspective, GA-based model method is proved genuine applicable That is because the GA procedure can optimize the whole parameters of the model even though some of them are ambiguous, however experimental method can just predict the parameters with explicit meaning like P _r and E _c. From accuracy perspective, the proposed model method is also testified to achieve satisfactory results. As a example, the proposed method have a less error (9.53%) than Nuffer’s model (24.52%) [10], while electric field fluctuates 3 ⋅ 10⁶ cycles. And for all the cases, the “Fits” of the GA-based model is enough accurate below 10%. Perceptually speaking, the simulated results agree with the experimental data very well no matter of how many cycling number is.

Fig. 3.

Experimental and simulated loops using proposed model method.

Table 2

Optimized parameters and simulation accuracy using proposed model method

Parameters	Optimized by GA / measured from test [10]
	0 cycle	3 ⋅ 10⁶ cycles	10⁸ cycles
P _s	0.3358 / —	0.2541 / —	0.253 / —
P _r	0.3348 / 0.336	0.2245 / 0.167	0.216 / 0.124
E _c	1.037 / 1.019	1.174/ 1.231	1.350 / 1.238
𝜉_r	4810 / —	170 / —	96.41 / —
d ₃₃₃	0.994 / —	0.8204 / —	0.268 / —
ϵ₀	1.798 / —	4.947 / —	3.889 / —
A	— / —	4.6e-5 / —	3.5e-5 / —
𝛼	— / —	0.055 / —	0.033 / —
c	— / —	0.003 / —	0.002 / —
Fit (%)	9.82 / —	9.53 / —	7.86 / —

Footnotes

Acknowledgements

This research is financially supported by the National Natural Science Foundation, People’s Republic of China, grant No. 11502188 and No. 11172226. These supports are gratefully acknowledged.

References

Kobune

Kuriyama

Furotani

, Ferroelectric materials and their applications, Japanese Journal of Applied Physics 54(10S) (2015).

Chen

Yusa

and Miya

, Some advances in numerical analysis techniques for quantitative electromagnetic nondestructive evaluation, Nondestructive Testing and Evaluation 24(1 & 2) (2009), 69–102.

Chen

Takagi

Uchimoto

Chigusa

and Yoshida

, Dependence of deformation-induced magnetic field on plastic deformation for SUS304 stainless steel, Int. J. Appl. Eelectromagn. Mech. 38(1) (2012), 17–26.

Hwang

S.C.

Huber

J.E.

McMeeking

R.M.

, The simulation of switching in polycrystalline ferroelectric ceramics, Journal of Applied Physics 84(3) (1998), 1530–1540.

Zhao

, Ferroelectric polarization and defect-dipole switching in an epitaxial (111) BiFeO₃ thin film, Journal of Applied Physics 117(20) (2015), 204103.

Lines

M.E.

and Glass

A.M.

, Principles and Applications of Ferroelectrics and Related Materials, Oxford University Press, 1977.

Miller

S.L.

and Mcwhorter

P.J.

, Physics of the ferroelectric nonvolatile memory field-effect transistor, Journal of Applied Physics 72(12) (1992), 5999–6010.

Shou-Wen

and Xi-Qiao

, Effects of electric fatigue on the butterfly curves of ferroelectric ceramics, Materials Science and Engineering: A 459.1 (2007), 273–277.

Ozgul

Trolier-McKinstry

and Randall

C.A.

, Influence of electrical cycling on polarization reversal processes in Pb (Zn_1∕3Nb_2∕3) O₃-PbTiO₃ ferroelectric single crystals as a function of orientation, Journal of Applied Physics 95(8) (2004), 4296–4302.

10.

Nuffer

Lupascu

D.C.

and Rödel

, Damage evolution in ferroelectric PZT induced by bipolar electric cycling, Acta Materialia 48(14) (2000), 3783–3794.

11.

Setter

and Waser

, Electroceramic materials, Acta Materialia 48(1) (2000), 151–178.