About Mathematical Models of Irreversible Polarization Processes of a Ferroelectric and Ferroelastic Polycrystals

Abstract

This chapter presents the prevalent mathematical models of irreversible processes in polycrystalline ferroelectric materials when they are subjected to intense electrical and mechanical influences. The main purpose of such models is to describe the dielectric hysteresis loops, with which the models of Rayleigh and Preisach coped well, though they were developed almost a 100 years ago. Nevertheless, in order to describe the whole gamut of material properties in irreversible polarization-depolarization processes, it was required in the last three decades to develop new approaches and methods that take into account the material structure and the physics of the process. In this chapter, we attempted to collect the most common one-dimensional models, with a view to give a brief description of the basics and approaches with the application of working formulas, algorithms and graphs of numerical calculations. On one-dimensional models, the basics of three-dimensional models are worked out, such as evolutionary laws, domains switching criteria, generalizations from “hysteron” to the polarization surface, and so on, so they are a necessary step in modeling. However, some of them proved to be so effective that they obtained the right to independent existence, as happened with the Preisach model, which found application in dynamic systems. This research is based on published articles, monographs, proceedings of conferences, and scientific reports of individual collectives published over the past 20–25 years.

1. Introduction

Polycrystalline ferroelectric materials or ceramics are active materials that, by virtue of their internal structure, have the ability to convert mechanical energy into electrical energy, and vice versa. This means that by acting on the sample an electric field or mechanical stress, we observe a response in the form of an electric displacement and the strain [1]. When the intensity of external electric fields or a mechanical stresses are small, the deformations and electrical displacements caused by them are also small. Such processes are called reversible; they say that there is a linear response. The modeling of such response leads to the construction of constitutive relations in the form of linear algebraic equations, which, like the generalized Hooke law, connect external and internal parameters. We can say that the mathematical model of the behavior of such materials is described by the linear algebraic operators in which the elements of elastic, piezoelectric and dielectric constant tensors are found experimentally. The overwhelming majority of problems concerning the calculation of the physical characteristics of transducers with polarized before saturation piezoceramic elements are solved within the linear response of active materials. In this case, the complete formulation of the problem includes the equations of motion, the equations of electrostatics, geometric relationships, and the constitutive relations. The general solutions obey the corresponding initial and boundary conditions. Such problems in the mathematical plan are linear. In simple cases, they can be solved analytically, and in more complex cases, we have to use numerical methods, for example, the finite element method.

The situation changes dramatically as soon as external loads reach thresholds, and their intensity continues to increase. In this case, irreversible processes begin, expressed in the fact that the response of the material will already be nonlinear. The consequence of this is that under increasing loads we have some nonlinear equations, while for decreasing ones we have other nonlinear equations. The constitutive relations become nonlinear and ambiguous. Mathematically, they can no longer be described by simple algebraic relations, but it is necessary to use operator relations of hysteresis type.

Due to the small volume of the article, we confine ourselves of modeling irreversible polarization-depolarization processes by an electric field and mechanical stresses under isothermal processes. Irreversible processes associated with relaxation properties, with the influence of temperature, with the features of the influence of size the ferroelectric granules, the dynamics of processes, and some others will not be considered here. The main circle of questions will be connected with the analysis of existing mathematical models describing the response of the material to external influences of high intensity for the isothermal process: we will consider the principles of constructing the constitutive relations, analyze them, and formulate some conclusions.

First, we note an interesting regularity: many irreversible processes have a similar response in the sense that the relationships between external and internal parameters are mathematically described by similar relationships. For example, in plasticity media, stresses cause elastic and residual deformations; in ferromagnets, the magnetic field leads to induced and remnant magnetization; in ferroelectrics, the electric field generates induced and residual polarization, etc. In the case of cyclic processes, the responses are described by hysteresis dependences, as shown in Figures 1, 2, 3. Therefore, it is not surprising that the mathematical methods being developed for the study of certain processes are often used to describe others. Mathematical modeling of nonlinear responses of polycrystalline ferroelectric materials plays a significant role [2], therefore, the creation and use of such models is based on experimental data.

Figure 1.

Deformation hysteresis loop.

Figure 2.

Magnetic hysteresis loop.

Figure 3.

Dielectric hysteresis loop.

2. Some experimental data

The criterion for the correctness and adequacy of the work of any model is the acceptable coincidence of the predicted data with experimental data. It should be noted that a qualitative experiment is a very complex study, so most of them reflect only certain properties with simple effects. We note only those works that reflect the electric and elastic response due to the action of the electric field and mechanical stresses. Basically, these are the works where the properties of ferroelectric ceramics of the perovskite type are investigated: for example, BaTiO₃, or a ceramics containing lead: PZT, PLZT 8/65/35. Interesting results [3] on the response of PLZT 8/65/35 on the effect of electric and mechanical fields and similar results [4, 5, 6, 7, 8, 9, 10] for complex loads show that the loops of electrical and deformation hysteresis depend significantly on the intensity of the operating fields. Uniaxial mechanical compressive stresses along the electric field axis affect the ability of domains to rotate: the more intense the mechanical stresses, the fewer domains are able switching along the electric field as it is shown in Figures 4 and 5.

Figure 4.

The effect of compressive mechanical stresses on the dielectric hysteresis loops.

Figure 5.

The effect of compressive mechanical stresses on “butterfly” dielectric hysteresis loops.

The types of considered ceramic materials are full ferroelectrics-ferroelastics, so the response of such materials to the effect of purely mechanical loads is of considerable interest. Some results of such tests can be seen in Figures 6 and 7. From this it is easy to see that under purely mechanical impacts, the “solid body”—“solid body” phase transition takes place, the material from isotropic (anisotropic) becomes anisotropic, the elastic modules of the material (tangents to the curves) change, residual strains appear that satisfy the property of incompressibility of the material (compare the values of longitudinal and transverse strains).

Figure 6.

Half-loop “stress-strain” in the case of axial compression.

Figure 7.

Half-loop “stress-transverse strain” in axial compression.

We note one more interesting property in the response of the material, when small loops of dielectric and strain hysteresis are investigated [9]. Small loops of dielectric and strain hysteresis due to the action of an electric field are shown in Figures 8 and 9. Small loops of strain hysteresis due to the action of mechanical stresses are shown in Figures 10, 11. Figure 10 reflects the stretching-compressing process, and Figure 11 reflects only pure compression followed by an increase in intensity.

Figure 8.

Small loops of dielectric hysteresis.

Figure 9.

Small loops of deformation hysteresis.

Figure 10.

Small loops “stress-strain” at axial stretching-compressing loads.

Figure 11.

Small loops “stress-strain” at axial compressing loads.

By behavior of small hysteresis loops, one can judge the change in the elastic, dielectric, and piezoelectric modules of the material. In addition, one can see the property of incompressibility of a material, by estimating the numerical values of strains at zero stresses. Summing up, we can say that the main task of mathematical modeling of irreversible processes is the construction of hysteresis operators taking into account the changing anisotropy of material properties.

3. Simple models

To compare the efficiency of a particular model, it is needed to have graphic images in the form of hysteresis loops. All the graphs given later were obtained in [2], where the algorithms are also described.

3.1. The Rayleigh model

One of the first models describing hysteresis relationships was the Rayleigh model [11], proposed in 1887 for magnetization processes of iron. Applying it to the polarization processes of polycrystalline ferroelectric continuum, it is necessary to replace the magnetic field by electric field and magnetization by polarization. Mathematically, the branches of the dielectric hysteresis are described by parabolic relationships:

$$P=\left\{\begin{array}{l}\alpha E^2,\mathit{\text{initial}}\mathit{\text{polarization}}\mathit{\text{curve}},\\ \alpha E_{\max }E\pm \frac{\alpha }{2}\left(E_{\max }^2-E^2\right),\mathit{\text{descending}}\mathit{\text{and}}\mathit{\text{ascending}}\mathit{\text{branches}},\end{array}\right.$$

where $\alpha =p_s/E_{\max }^2$, $p_s$—spontaneous polarization, $E_{\max }$—maximum value of the electric field, and the hysteresis loop is shown in Figure 12.

Figure 12.

The Rayleigh model.

Conclusion: a quadratic dependence between the electric field and polarization approximates the dielectric hysteresis for small and medium intensity values of the electric field.

3.2. Evolutionary models

Investigating the nonlinear behavior of ferroelectric materials, the authors [12] assumed that the dielectric displacement $D$ is a function of the electric field $E$, the dipole moment $\mu$, and the number of dipoles per unit volume $N$ aligned in the direction of the electric displacement. Thus, irreversible parameters $\mu ,N$ have been introduced into the model, and evolutionary laws are formulated in order to identify them:

$$D=\tilde{D}\left(E,\mu ,N\right);\dot{N}=g\left(E,\mu ,N\right),g\left(0,\mu ,N\right)=0.$$

Under the assumption of linear functions entering into the previous relations and the linear dependence of the dipole moment on the electric field

$$\begin{array}{l}D={\left(\tilde{D}\right)}_0+{\left({\tilde{D}}_E\right)}_{0}E+{\left({\tilde{D}}_{\mu }\right)}_0\left(\mu -{\mu }_0\right)+{\left({\tilde{D}}_N\right)}_0\left(N-N_0\right);\\ \dot{N}={\left(g_E\right)}_{0}E+{\left(g_{\mu }\right)}_0\left(\mu -{\mu }_0\right)+{\left(g_N\right)}_0\left(N-N_0\right);\mu \left(\tau \right)-{\mu }_0\mathit{\alpha E}\left(\tau \right),\end{array}$$

a formula for the electric displacement was obtained

$$D={\left(\tilde{D}\right)}_0+\tilde{\epsilon }\left(0\right)E+\underset{0}{\overset{1}{\int }}\frac{d}{\mathit{dt}}\tilde{\epsilon }\left(t-\tau \right)E\left(\tau \right)d\tau ;\tilde{\epsilon }\left(t\right)=\tilde{\tilde{\epsilon }}\left(t\right)+\alpha \tilde{\tilde{\epsilon }}\left(t\right).$$

The creep functions $\tilde{\tilde{\epsilon }}\left(t\right),\tilde{\tilde{\epsilon }}\left(t\right)$ have an exponential form, and because of their cumbersomeness, they are not given here. Similar arguments are also made for the related problem, where the mechanical stress and the electric displacement are presented in the form of similar integral dependences with the same exponential creep functions. To determine the elastic and dielectric properties [13] of the ceramic material PZT 65/35, a connection between the elastic and dielectric constants and the speed of sound in polarized and unpolarized ceramics was established. Similar studies can be found in [14, 15, 16, 17, 18, 19].

Conclusion: irreversible parameters were introduced in implicit form, for the find of which evolutionary laws were applied. The constitutive relations were obtained in the form of creep integrals. The creep functions are constructed under the assumption of linear dependences for functions entering into evolutionary dependencies. For the formulation of increasing and falling branches of hysteresis, inequalities for functions describing the load and response of the material were formulated.

3.3. Models of the theory of plasticity

In view of the similarity of the plasticity and polarization phenomena noted earlier, for qualitative description of the polarization effects, use the rheological models of the theory of plasticity. At the basis of our subsequent actions lie analogues of mechanical and electrical quantities: the generalized coordinate—the electric charge; generalized speed—current; coefficient of elastic compliance—capacity; the generalized force is the electromotive force. In the transition to continuous media, forces are replaced by mechanical stresses, displacements by strains, etc. As a result, one can write the following correspondence: $\sigma \leftrightarrow E$, $\epsilon \leftrightarrow P$, where $E$ is the electric field; $P$ is polarization; $\sigma$ is mechanical stress; $\epsilon$ is strain. The elastic element is associated with a capacitor and the element of dry friction is a bipolar zener diode (Figure 13). The rheological formulas for a capacitor and a zener diode can conveniently be described using differential inclusions [20] which can be represented by the following expressions:

$$E_e=\frac{1}{c}{P}_e;E_0\in \mathit{rS}\left({\dot{P}}_0\right)\left(c,r>0-\mathit{\text{const}}\right),$$

Figure 13.

Capacitor and bipolar zener diode elements.

where the indices “e” and “0” indicate the induced and residual components, respectively, and $S\left(v\right)$ is a function of sets determined by the rule:

$$S\left(v\right)≔\left\{\begin{array}{l}\left\{-1\right\},v<0;\\ \left[-1,+1\right],v=0;\\ \left\{+1\right\},v>0.\end{array}\right.$$

By connecting these elements in series (Figure 14a), according to the conditions $E=E_e+E_0,P=P_e=P_0$, we obtain a differential inclusion $E\in \frac{1}{c}P+\mathit{rS}\left(\dot{P}\right)$ that determines the “backlash” operator. If these elements are connected in parallel (Figure 14b), then from the conditions $E=E_e=E_0,P=P_e+P_0$ it is easy to derive a differential inclusion $\dot{P}\in c\dot{E}+S^{-1}\left(\frac{E}{r}\right)$ that determines the “stop” operator [20]. Next, one can determine the Prager polarization models by adding capacitors to the considered chains, as shown in Figures 15a, 15b. Not stopping the detailed description of the Prager model, we note only the case shown in Figure 15b, which is described by the following rheological formula [2]:

$$\left(1+\frac{c_1}{c_2}\right)\dot{P}-c_1\dot{E}\in S^{-1}\left(\frac{E-\left(1/c_2\right)P}{r}\right).$$

Figure 14.

Connected elements: a) in series; b) in parallel.

Figure 15.

Prager polarization models: a) capacitor in parallel; b) capacitor in series.

If we collect a battery of one of the chains in parallel, we see that the resulting polarization for such a connection is determined by summing the individual polarization components of each chain, and the electric field is the same for any of them. Let the chains be distributed continuously [2, 21, 22], then there is a distribution density function $f\left(x\right):d\mu \left(x\right)=f\left(x\right)\mathit{dx}$ and three families of positive parameters $\alpha \left(x\right),\beta \left(x\right),r\left(x\right)$. For any $x\in X$ exist a generalized Prager operator ${\left.(E,\{E\left(0\right)-\alpha P^0\left(x\right)\in [-r,+r\left]\right\}\right|}_{x\in X})\mapsto P$, which is defined by the following system:

$$\begin{array}{l}\dot{P}\left(x\right)-\beta \left(x\right)\dot{E}\in S^{-1}\left(\frac{E-\alpha \left(x\right)P\left(x\right)}{r\left(x\right)}\right);P=\underset{X}{\int }P\left(x\right)d\mu \left(x\right);\\ {\left.E\left(t\right)-\alpha P(x,t)\right|}_{t=0}=E\left(0\right)-\alpha P^0\left(x\right)\left(\forall x\in X\right).\end{array}$$

Taking a battery of continuously distributed chains for a representative volume, we can determine it polarization for any value of the electric field. It is possible to perform calculations using the generalized Prager model if, in addition to the probability function of the distribution density, three functions $\alpha \left(x\right),\beta \left(x\right),r\left(x\right)$ are also given. We give one numerical experiment [2], by selecting appropriately incoming in it functions. Let the electric field vary in time according to the law of the sine. The probability density function is chosen as a function of the normal Gaussian distribution:

$$f\left(x\right)=\frac{1}{\sqrt{2\pi }\sigma }\text{exp}\left\{-\frac{{\left(x-b\right)}^2}{2{\sigma }^2}\right\}.$$

We will assume that the generalized Prager operator is generated by a battery of parallel connected chains with functions of the following form:

$$\begin{array}{l}\alpha \left(x\right)=\frac{1}{c_1\left(x\right)+c_2\left(x\right)},\beta \left(x\right)=c_2\left(x\right),c_1\left(x\right)=a\frac{2-\mathrm{th}\mid x\mid }{2};\\ c_2\left(x\right)=b\frac{2-\mathrm{th}\mid x\mid }{2};r\left(x\right)=c\frac{1+\mathrm{th}\mid x\mid }{2}\left(a,b,c-\mathit{\text{const}}\right).\end{array}$$

The result of the calculations is shown in Figure 16.

Figure 16.

Generalized Prager model.

Conclusion: varying the parameters of the Gaussian distribution, as well as the functions included in the differential inclusions, it is possible to substantially change the shape of the hysteresis loop. This circumstance gives confidence that by the selection of these parameters a good coincidence of the modeled and experimental hysteresis loop in the one-dimensional case can be achieved.

3.4. The Preisach model

This model was proposed in 1935 by Preisach [23]. The model is based on the concept of hysteron [24, 25], which approximates the switching of 180° domain. Simulation of the polarization of ceramics is associated with its representation as a set of a very large number of 180° hysterons with the probability density function $\mu \left(x,y\right):$

$$\underset{x\ge 0}{\iint }\mu \left(x,y\right)\mathit{\text{dxdy}}=1\left(x=E_i,y=E_c\right),$$

defined on the plane $\left\{|x|\le \infty ,y>\infty \right\}$. The inhomogeneity of physical conditions in different parts of the ceramics generates a large scatter of hysterons at coercive and internal fields, and the switching of each domain is described by a rectangular hysteresis loop with its coercive ($E_c$) and internal ($E_i$) fields (Figure 17). For each state in the positive half-plane, there is a boundary separating the domains of two opposite directions. For an unpolarized state, it is the abscissa axis (Figure 18). When an electric field of definite sign is applied, this boundary moves due to the involvement of new hysterons in the switching process (Figure 19). If the electric field changes the growth direction, then the direction of the boundary movement also changes (Figure 20).

Figure 17.

Hysteron.

Figure 18.

Areas with nonzero density of hysterons for opposite directions of polarization.

Figure 19.

The motion of the separation boundary with increasing field.

Figure 20.

The motion of the separation boundary with decreasing field.

The distribution function in the vicinity of the coercive fields had to have a pronounced peak, which allows approximating its using known distributions, for example Gauss, with subsequent determination of the parameters entering into it. There are also distribution functions in the form of polynomials degree of intense and coercive fields [26, 27, 28, 29]. For a convenient mathematical representation, we introduce the concept of an elementary dipole hysteresis operator (Figure 21) or the “relay” operator [22] using the functions of sets. Suppose that for any fixed pair $\alpha ,\beta \in R\left(\alpha <\beta \right)$, any function $E\in C^0\left[0,T\right]$, and any of the values $\xi =-1$, or $\xi =+1$ a function is defined $p:\left[0,T\right]\to \left\{-1,+1\right\}$:

$$p\in \left\{\begin{array}{lll}p\left(0\right),& \mathit{\text{if}}& G_t=0;\\ -1,& \mathit{\text{if}}& G_t\ne 0,E\left(\text{max}{G}_t\right)=\alpha ;\\ +1,& \mathit{\text{if}}& G_t\ne 0,E\left(\text{max}{G}_t\right)=\beta ,\end{array}\right.p\left(0\right)\in \left\{\begin{array}{lll}-1,& \mathit{\text{if}}& E\left(0\right)\le \alpha ;\\ \xi ,& \mathit{\text{if}}& \alpha \le E\left(0\right)\le \beta ;\\ +1,& \mathit{\text{if}}& E\left(0\right)\ge \beta ,\end{array}\right.$$

Figure 21.

The definition of elementary dipole hysteresis operator.

and $G_t$—the set of singular time points, such that for $t\in \left(0,T\right]$

$$G_t≔\left\{\tau \in \left(0,t\right]:E\left(\tau \right)=\alpha \mathit{\text{or}}\beta \right\}.$$

In view of the finiteness of the oscillations of any continuous function between $\alpha$ and $\beta$, it follows that $p\left(t\right)$ can have a limited number of jumps between $\left\{+1\right\},\left\{-1\right\}$. This allows us to determine the “relay” rate-independent operator

$${\gamma }_{\alpha \beta }:C^0\left[o,T\right]\times \left\{-1,+1\right\}\to \left\{-1,+1\right\}.$$

The parameters $\alpha$, $\beta$ and $x$, $y$ for the hysteron are related by the following linear relations: $\alpha =y-x;\beta =y+x$. Using the “relay” operator and the distribution function of the hysterons, the irreversible polarization can be determined by the integral:

$$P\left(t\right)=p_s^{\ast }\underset{x\ge 0}{\iiint }\mu \left(x,y\right)\left({\widehat{\gamma }}_{\alpha \beta }E\right)\left(t\right)\mathit{\text{dxdy}},$$

where $p_s^{\ast }$—the maximum polarization value achievable in the polarization process of ceramic by a homogeneous electric field. In addition, a notation ${\gamma }_{x-y,x+y}={\widehat{\gamma }}_{\mathit{xy}}$ is introduced here. It is noteworthy that if we choose an equable distribution function

$$\mu \left(x,y\right)=\left\{\begin{array}{l}1/E_{\max }^2\mathit{at}\left(x,y\right)\in S_{\mathrm{\Delta }};\\ 0\mathit{at}\left(x,y\right)\notin S_{\mathrm{\Delta }};\end{array}\right.$$

where $S_{\mathrm{\Delta }}$—the region on the half-plane of the variables $x$ and $y$, indicated in Figure 18 by a triangle, then by simple calculations of the integrals we easily find hysteresis dependences of the Rayleigh method.

Various forms of loops of dielectric hysteresis can be found in [2]. For example, if we neglect the correlation coefficient in the Gaussian distribution and use the symmetry condition for the function with respect to the variable $y$, we obtain

$$\mu \left(x,y\right)=\frac{1}{\sqrt{2\pi }{\sigma }_1}\text{exp}\left[\frac{-{\left(x-a_1\right)}^2}{2{\sigma }_1^2}\right]\frac{1}{\sqrt{2\pi }{\sigma }_2}\text{exp}\left(\frac{-y^2}{2{\sigma }_2^2}\right).$$

To estimate the influence of the parameter $a_1$ on the behavior of the curves of the hysteresis loop, we set $E_{\max }=2\cdot {10}^{6}V/m;{\sigma }_1={\sigma }_2=2\cdot {10}^{2}V/m;$ and we will increase the parameter $a_1$ sequentially, assuming $a_1=0;1\cdot {10}^6;2\cdot {10}^6;5\cdot {10}^6;7\cdot {10}^{6}V/m$. The resulting hysteresis loops are shown in Figure 22, where the cases a, b, c, d, and e correspond to the indicated values of the parameter.

Figure 22.

Preisach model: effect of increasing the parameter $a_1$ on the form of hysteresis loops.

A Preisach model has become widely used not only in the description of magnetic and ferroelectric hysteresis, but it is intensively used in calculating the damping coefficients of many dynamical systems, including the dynamic of dipole switching [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48]. Two-dimensional models were also developed [49, 50, 51, 52, 53, 54, 55, 56].

Conclusion: Preisach model in the simplest case “includes” the domain structure, but operates only 180° domains. In the mathematical plan, the approximation of the real loop is carried out by elementary rectangular hysteresis loops. A model is intended only for finding the residual polarization. This model does not include mechanical stresses, which significantly reduces its practical application.

3.5. Model of orientation switching

A representative volume is considered as a set of $N$ domains oriented in space in an arbitrary manner [57]. Here, not only 180°, but also 90° switching are considered. Let the applied electric field has a tension $\mathbf{E}$ and each domain is characterized by vector of spontaneous polarization ${\mathbf{p}}_S$ and crystallographic axes $\mathbf{a},{\mathbf{a}}^{\text{'}},\mathbf{c}$ (Figure 23). By simple averaging one can determine the vector of residual polarization and the tensor of strain:

$${\mathbf{P}}_{\mathbf{0}}\equiv <{\mathbf{p}}_S>=\frac{p_S}{N}\sum _k{\left(\mathbf{с}\right)}_k;{\mathbf{\epsilon }}_{\mathbf{0}}\equiv <{\mathbf{\epsilon }}_S>=\frac{{\epsilon }_S}{N}\sum _k{\left(\mathbf{с}\otimes \mathbf{с}-\frac{1}{2}\mathbf{a}\otimes \mathbf{a}-\frac{1}{2}\mathbf{b}\otimes \mathbf{b}\right)}_k.$$

Figure 23.

Determination of angles in the model of orientation switching.

Denote the plane perpendicular to the vector $\mathbf{c}$ through B, and the plane passing through the vectors $\mathbf{E}$ and $\mathbf{c}$ by A. We introduce the following angles: the angle ${\gamma }_c$ between the direction of the field and the axis $\mathbf{c}$; the angle ${\gamma }_a$ between the vector $\mathbf{a}$ and the field $\mathbf{E}$; angle $\omega$ between the axis $\mathbf{a}$ and the line of intersection $\mathit{OK}$. It is obvious that the angles introduced are within

$$0\le {\gamma }_c\le \pi ;0\le \omega \le \pi /4;\omega \le {\gamma }_a\le \pi /2.$$

Let $E_{\mathit{cc}},E_{\mathit{ca}}$ be the coercive fields of 180° and 90° switching, respectively. The main conditions of 180° domains switching (turn) are

$$E\text{cos}{\gamma }_c\ge E_{\mathit{cc}};\frac{E\text{cos}{\gamma }_c}{E_{\mathit{cc}}}-\frac{E\text{cos}{\gamma }_a}{E_{\mathit{ca}}}\ge 0.$$

For 90° domains switching, we have $\frac{E\text{cos}{\gamma }_a}{E_{\mathit{ca}}}-\frac{E\text{cos}{\gamma }_c}{E_{\mathit{cc}}}\ge 1$, if ${\gamma }_c\le \pi /2$ and $E\text{cos}{\gamma }_a\ge E_{\mathit{ca}};\frac{E\text{cos}{\gamma }_a}{E_{\mathit{ca}}}-\frac{E\text{cos}{\gamma }_c}{E_{\mathit{cc}}}\ge 0$, if ${\gamma }_c>\pi /2$.

A schematic representation of all domains axes $\mathbf{c}$ before and after polarization can be seen in Figure 24. Dielectric hysteresis loop are constructed for quasi-static processes, that is, for a sequence of equilibrium states $\left\{E_i\right\}$. With this goal, the process of loading by an electric field $E=E\left(t\right)$ is replaced by a sequence of values of the electric field $E\in \left\{E_i\right\}:E_i=E\left(t_i\right)$, and for each state $E=E_i$, the domain switching conditions are checked, after which the residual polarization and the residual strain are calculated by simple averaging.

Figure 24.

Schematic representation of domains axes $\mathbf{c}$ before and after polarization.

As an example [2], a hysteresis loop (Figure 25) was calculated for the next set of parameters: $N=1273248;E_{\mathit{cc}}=2\cdot {10}^{6}V/m;E_{\mathit{ca}}=3\cdot {10}^{6}V/m;E_{\max }=6\cdot {10}^{6}V/m$, where $E_{\max }$—maximum value of the electric field. It is interesting to note that a model makes it possible to determine the value of the angle of cone in which the directions of all vectors of spontaneous polarization are distributed after the removal of the electric field. With the abovementioned numerical values of the parameters, we obtained: $\alpha ={129}^0;P_{01}=-0.18\cdot {10}^{-4}{p}_s;P_{02}=-0.79\cdot {10}^{-5}{p}_s;P_{03}=0.81p_s$.

Figure 25.

Model of orientation switching.

Conclusion: in physical terms, this model is closer to the physics of the phenomenon of polarization of polycrystalline ferroelectrics, but does not take into account the effect of neighboring domains on each other during the polarization process. Because of this, the loop acquires angular shapes and is not suitable for describing the differential properties of the material. On the other hand, it allows us to find a value of the angle in which the directions of all the spontaneous polarization vectors are located after reaching the saturation polarization. Another drawback of the model is that it does not include mechanical stresses, which does not allow us to investigate ferroelastic phenomena.

3.6. Energy model of switching

In this model, a representative volume is also considered, including a set of domains oriented in space in an arbitrary manner [2] and a residual polarization vector is determined by averaging. The main difference is in describing the orientation of domains and determining the criteria for switching domains. For example, for a ferroelectric of the perovskite type, the axes of the local system are introduced $\mathbf{a}\mathbf{,}\mathbf{b}\mathbf{,}\mathbf{c}$, which are the axes of the crystallographic system directed along the ribs of a rectangular parallelepiped, as shown in Figure 26. The local system $\mathbf{a}\mathbf{,}\mathbf{b}\mathbf{,}\mathbf{c}$ with respect to a fixed system ${\mathit{Ox}}_1{x}_2{x}_3$ is determined by three angles $\phi ,\psi ,\omega$. The application of an electric field or mechanical stresses of high intensity causes a process of domain switching. However, the switching is consistent with the crystallographic axes of the ferroelectrics and is possible only in certain directions (Figure 27). Each domain in the field of external loads (electric field and mechanical stresses) possesses some energy. When the loading factors change, the energy of each of the domains will change. If a domain has to switch, then that switching will be only in such direction where its energy minimal. The switching process starts only if a difference of domain energy of current state and the state with the minimum energy exceeds the threshold value [1, 58, 59]:

$$-{\mathbf{p}}_{\mathbf{S}}\cdot \mathbf{E}+{\mathbf{p}}_{\mathbf{S}}^{\mathbf{min}}\cdot \mathbf{E}-{\mathbf{\epsilon }}_{\mathbf{S}}:\mathbf{\sigma }+{\mathbf{\epsilon }}_{\mathbf{S}}^{\mathbf{min}}:\mathbf{\sigma }\ge U_c,$$

Figure 26.

Determination of the angles of the transition from the local to the global coordinate system.

Figure 27.

Possible positions of the spontaneous polarization vector and deformation of the unit cell.

where $\mathbf{E}$ is the electric field vector; $\mathbf{\sigma }$ is tensor of mechanical stresses; ${\mathbf{p}}_{\mathbf{S}}$ is vector of spontaneous polarization, ${\mathbf{\epsilon }}_{\mathbf{S}}$ is tensor of spontaneous strain, and $U_c$ is the threshold value of energy per unit volume. The form of the dielectric hysteresis loop practically coincides with the previous case; therefore, it is not given here.

Conclusion: the energy model of switching at physical essence takes into account the physical phenomena of polarization of polycrystalline ferroelectrics but also does not take into account the effect of neighboring domains on each other during the switching process. The dielectric hysteresis loop has angular shapes and is not suitable for describing the differential properties of a material. In addition, this model operates only with the residual parameters of polarization and strain and does not include induced components. However, this model has an unquestionable advantage because it already includes a mechanical stresses, which allows us to describe not only ferroelectric, but also ferroelastic phenomena.

3.7. The Giles-Atherton model

The model is based on the so-called “limiting” (or anhysteretic) curve, derived analytically on the basis of Weiss theory and Boltzmann statistics [60, 61, 62, 63, 64, 65, 66, 67]. If in the polycrystalline ferroelectric material were no mechanisms for locking the walls of the domains, then after the removal of the electric field, the polarization of the representative volume would be zero. With great assumptions, one could say that such a material would behave like polar liquid, where the electric field tries to align dipoles along one direction, and the thermal fields impede it: impacts try to destroy the alignment pattern, and after removing the electric field, the thermal fields neutralize the polarization. Of course, a different domain switching mechanism is observed in polycrystalline ferroelectrics, but the switching model is borrowed from polar liquids.

To determine the dependence connecting the polarization of a representative volume with the strength of the effective electric field [60] let the representative volume of the ferroelectric material contain $N$ domains with a spontaneous vector ${\mathbf{p}}_{\mathbf{s}}$. For the averaging operation, a point of reduction is chosen in space, and to each vector of spontaneous polarization is placed in conformity a unit vector parallel to it. The resulting vector of polarization is found by determining the area of the unit sphere onto which the ends of the unit vectors are exit. In an unpolarized state, all unit vectors are distributed uniformly over the entire surface, so that the resulting polarization is zero.

Using the positions of the Weiss theory, we assume that the field of forces acting on the ferroelectric domain is a sum of the electric field $\mathbf{E}$ and some “molecular field” proportional to its polarization: ${\mathbf{E}}^{\mathit{ef}}=\mathbf{E}+\alpha {\mathbf{P}}_0$, where $\alpha$ is a certain constant.

The potential energy of a domain that is similar to a dipole depends on the direction of its polarization vector and is expressed by the formula: $U=-{\mathbf{p}}_s\cdot {\mathbf{E}}^{\mathit{ef}}V$, where $V$—domain volume measure.

To determine the distribution of the domain axes in the presence of an effective field orienting them, it is necessary to use the Boltzmann theorem of statistical mechanics. Under the conditions of thermodynamic equilibrium, the distribution law for domains in the presence of a conservative electrostatic field differs from the law of their distribution in the absence of this field by a factor $\text{exp}\left(-U/k_{\ast }T\right)$, where $T$ is the absolute temperature, $k_{\ast }=1.3807\cdot {10}^{-23}\mathrm{J}/\mathrm{K}$ is the Boltzmann constant. Therefore, the averaging operation yields the maximum possible polarization

$${\mathbf{P}}_{\infty }=\frac{{\int }_0^{2\pi }d\phi {\int }_0^{\pi }\text{exp}\left(\frac{{\mathbf{E}}^{\mathit{ef}}\cdot {\mathbf{p}}_s}{k_{\ast }T}\right){\mathbf{p}}_s\text{sin}\psi d\psi }{{\int }_0^{2\pi }d\phi {\int }_0^{\pi }\text{exp}\left(\frac{{\mathbf{E}}^{\mathit{ef}}\cdot {\mathbf{p}}_s}{k_{\ast }T}\right)\text{sin}\psi d\psi }.$$

In the one-dimensional theory, the direction of the electric field coincides with the direction of the axis $\mathit{Oz}$ of the Cartesian fixed coordinate system, and the residual polarization in this direction is estimated. Taking into account this projection and calculating the integral, we determine the maximum residual polarization $P_{\infty }$ of the representative volume in the form of a Langevin distribution:

$$P_{\infty }=p_s\left[\mathrm{cth}\left(\frac{E+\alpha P_0}{a}\right)-\frac{a}{E+\alpha P_0}\right],a=\frac{k_{\ast }T}{p_{s}V}.$$

In this model, a reversible (induced) part $P_e$ as well as the residual part $P_0$ of the polarization is taken into account:

$$P=P_e+P_0.$$

The reversible part is the state parameter. It can be defined as:

$$P_e=с\left(P_{\infty }-P_0\right),$$

where $с$ is still an indeterminate factor. The irreversible part of the polarization is the parameter of the process. Therefore, to determine it, we can use the evolutionary law, which can be written in differentials:

$${\mathit{dP}}_0=\frac{1}{k\delta }\left(P_{\infty }-P_0\right){\mathit{dE}}^{\mathit{ef}},$$

where $\delta =\text{sign}\left(\mathit{dE}\right)$, and $k$ is a positive constant, also subjecting to determination. Passing from the effective to the applied field, we obtain an ordinary differential equation:

$$\frac{{\mathit{dP}}_0}{\mathit{dE}}=\frac{P_{\infty }-P_0}{k\delta -\alpha \left(P_{\infty }-P_0\right)}.$$

To solve the equation, we need use numerical methods, for example, the fourth-order Runge-Kutta method. When the electric field varies according to the harmonic law, we obtain dielectric hysteresis loops. The model includes five parameters: $p_s,\alpha ,c,a,k$, which are found from the condition of coincidence of calculated and experimental data. Numerous experiments were carried out in [2], and the effect of the model parameters on the shape of loops was investigated. It is shown that the coefficient $\alpha$ is responsible for the amplitude of the loop, the coefficient $a$—for its slope, the coefficient $k$ for the loop area, the coefficient $c$ for the flatness of the loop. Varying the values of the coefficients, one can achieve not only a qualitative but also a quantitative coincidence with the experimental data. It is shown that it is possible to start up the procedure for determining the model coefficients by the coincidence of the calculated values with the experimental data in the set of points.

Some examples of numerical calculations in which the influence of the coefficient k on the form of the hysteresis curves is studied are shown in Figures 28, 29, 30.

Figure 28.

The Giles-Atherton model: the smaller the coefficient $k$, the narrower the hysteresis loop.

Figure 29.

Model Giles-Atherton: with increasing coefficient $k$, the slope of the hysteresis loop changes.

Figure 30.

The Giles-Atherton model: with increasing coefficient $k$, the amplitude of the hysteresis loop changes.

Conclusion: Giles-Atherton model deals with a huge number of distributed domains in a representative volume. Using of Boltzmann statistics made it possible, first, to construct an extremely possible polarization in the form of a Langevin distribution; second, to obtain the density of distribution of domains as a function of exponential type; third, to use the switching criterion. In this case, the switching is reduced to the energy rotations of the domains in dependency of the intensity of the electric field. The model contains both a reversible and a residual part of the polarization. The model contains a set of five parameters, which are chosen from the condition of coincidence of the calculated large hysteresis loop and experimental data. Due to shortcomings of the model, one can be attributed include a very rough approximation in the description of domain rotations in ferroelectrics, where this process is replaced by the process of rotation of dipoles in polar liquids in the construction of limiting polarization. This immediately affects when trying to build small hysteresis loops, where the results of calculations lead to large discrepancies.

3.8. Modeling of ferroelastics

It is interesting to note that, despite its shortcomings, the Giles-Atherton model can be used to describe deformation of ferroelastics [68]. Repeating the main contents of the previous model, we can introduce the limiting strain, which in the one-dimensional case of stretching-compression reduces to calculating the probability integral of the imaginary argument:

$${\mathrm{\epsilon }}_{\infty }=\frac{{\mathrm{\epsilon }}_s}{2}\frac{\underset{-1}{\overset{1}{\int }}\text{exp}\left[\frac{3\left(\mathrm{\sigma }+\beta \mathrm{\epsilon }\right)}{2a}{t}^2\right]t^{2}dt}{\underset{-1}{\overset{1}{\int }}\text{exp}\left[\frac{3\left(\mathrm{\sigma }+\beta \mathrm{\epsilon }\right)}{2a}{t}^2\right]dt},a=k_{\ast }T/{\mathrm{\epsilon }}_{s}V.$$

The total deformation consists of the reversible ${\epsilon }_e$ and irreversible ${\epsilon }_0$ parts: $\epsilon ={\epsilon }_e+{\epsilon }_0$, and to determine the irreversible part, we obtain an ordinary differential equation

$$\frac{d{\epsilon }_0}{d\sigma }=\frac{{\epsilon }_{\infty }-{\epsilon }_0}{k\delta -\alpha \left({\epsilon }_{\infty }-{\epsilon }_0\right)}.$$

The large and small loops of strain hysteresis, calculated from this model, can be seen in Figures 31, 32.

Figure 31.

Modeling of ferroelastics: large hysteresis loop.

Figure 32.

Modeling of ferroelastics: small hysteresis loops.

In [68, 69], two modifications were made for these models in order to obtain small hysteresis loops. First, to the induced component of polarization and strain was added a summand proportional to the loading field. Second, restrictions on the process of domain switching during the construction of small cycles are introduced. The results of this modernization are shown in Figures 33, 34. The obtained results are in good agreement with the experimental data shown in Figures 8 and 11, not only qualitatively, but also quantitatively.

Figure 33.

Improved Giles-Atherton model: small loops of dielectric hysteresis.

Figure 34.

Improved Giles-Atherton model: small loops of deformation hysteresis.

Conclusion: the dignities and drawbacks of the model for ferroelastics are the same as those of the Giles-Atherton model.

4. Discussion

Each of the presented models performs the basic function related to the description of hysteresis dependencies. However, each of them is based on different prerequisites. Therefore, the results of the work of a particular model differ from the corresponding experimental data. The Rayleigh model gives a coincidence with the experimental data for small and medium fields. The Preisach model is widely used in areas where the integral properties of hysteresis are important. The model of orientation switching and model of the energy switching allow us to determine the cone of the angles all domains for which the directions of the vectors of the spontaneous polarization lined up inside this cone after the polarization to saturation, but the hysteresis loops have angular shapes. The Giles-Atherton model, taking into account the modernization, allows describing not only large but also small hysteresis loops.

Advantages and disadvantages of each of the considered models were presented earlier. Each of the models has a number of parameters that influence on the shape of the hysteresis loops. To determine these parameters, it is necessary to compare the calculated and experimental data. However, in the mathematical plan, this model belongs to the class of inverse problems and has a number of specific difficulties.

To sum up, it can be noted that the one-dimensional models allow describing hysteresis dependences of irreversible polarization processes. The more accurate the model is based on the physics of the phenomenon, the more accurate the results of its work. For example, those models that include lot of domains are more accurate in describing quantitative dependencies. However, analysis shows that at the present moment a universal mathematical model of the polarization and deformation of polycrystalline ferroelectric materials has not yet been developed, which would accurately reflect the data of the experimental dependences. Many problems that are importance remain unresolved. Unresolved problems include the following:

• small loops of dielectric and strain hysteresis have been poorly investigated;

• existing models do not reflect functional dependencies on the effects of any loading components with simultaneous exposure to electric fields and mechanical loads;

• the question of the validity of the application of the developed models to thin films remains open;

• there are no studies of how the anisotropy of the material changes with the simultaneous action of an electric field and mechanical stresses;

• very poorly represented models, which take into account the effect of neighboring domains on current switching;

• there are practically no studies related to the analysis of the similarity and differences in the existing models of Preisach, Rayleigh, Giles-Atherton, plasticity materials, except for [70, 71];

• in fact there is no mathematical analysis of the effect of model parameters on the final result, there is no formulation of a minimum set of parameters;

• very poorly represented models, which use the functions of the density of the distribution of domains.

5. Conclusions

One-dimensional models that reflect the hysteretic properties of polycrystalline ferroelectric media on the external effects of high-intensity electric and mechanical fields are considered. The bases of construction of each model are disassembled. The fundamentals of the construction of each of the well-known Rayleigh models, evolution models, models of plasticity theory, Preisach models, models of orientation switching, energy switching models, the Giles-Atherton model are analyzed and the results of their work in the form of hysteresis loops are presented. The main advantages and disadvantages of each model are noted; in particular, the differential and integral properties of hysteresis dependences are marked. It is noted that each of the models has a number of parameters, the choice of which affects the shape of the hysteresis loop. However, in the mathematical sense, such a choice generates additional difficulties, which are connected with the solution of inverse problems. Unresolved problems in the general problem of modeling the polarization of polycrystalline ferroelectric materials are noted. A list of works is given in which the main ideas and results of each of the earlier models are reflected.

Acknowledgments

This work was supported by the Russian Foundation for Basic Research (Grant 17-08-00860-a).