Open access peer-reviewed chapter

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

Written By

Alexander Skaliukh

Submitted: 12 January 2018 Reviewed: 02 May 2018 Published: 03 October 2018

DOI: 10.5772/intechopen.78262

Chapter metrics overview

1,093 Chapter Downloads

View Full Metrics


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.


  • mathematical model
  • hysteresis
  • domains
  • polarization
  • strain
  • ferroelectric materials
  • ceramics
  • constitutive relations

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, BaTiO3, 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:


where α=ps/Emax2, ps—spontaneous polarization, Emax—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 μ, and the number of dipoles per unit volume N aligned in the direction of the electric displacement. Thus, irreversible parameters μ,N have been introduced into the model, and evolutionary laws are formulated in order to identify them:


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


a formula for the electric displacement was obtained


The creep functions ε˜˜t,ε˜˜t 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: σE, εP, where E is the electric field; P is polarization; σ is mechanical stress; ε 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:


Figure 13.

Capacitor and bipolar zener diode elements.

where the indices “e” and “0” indicate the induced and residual components, respectively, and Sv is a function of sets determined by the rule:


By connecting these elements in series (Figure 14a), according to the conditions E=Ee+E0,P=Pe=P0, we obtain a differential inclusion E1cP+rSṖ that determines the “backlash” operator. If these elements are connected in parallel (Figure 14b), then from the conditions E=Ee=E0,P=Pe+P0 it is easy to derive a differential inclusion ṖcĖ+S1Er 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]:


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 fx:dμx=fxdx and three families of positive parameters αx,βx,rx. For any xX exist a generalized Prager operator (E{E0αP0x[r+r]}xX)P, which is defined by the following system:


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 αx,βx,rx 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:


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


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 μxy:


defined on the plane xy>. 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 (Ec) and internal (Ei) 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.


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 α,βRα<β, any function EC00T, and any of the values ξ=1, or ξ=+1 a function is defined p:0T1+1:


Figure 21.

The definition of elementary dipole hysteresis operator.

and Gt—the set of singular time points, such that for t0T


In view of the finiteness of the oscillations of any continuous function between α and β, it follows that pt can have a limited number of jumps between +1,1. This allows us to determine the “relay” rate-independent operator


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


where ps—the maximum polarization value achievable in the polarization process of ceramic by a homogeneous electric field. In addition, a notation γxy,x+y=γ̂xy is introduced here. It is noteworthy that if we choose an equable distribution function


where SΔ—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


To estimate the influence of the parameter a1 on the behavior of the curves of the hysteresis loop, we set Emax=2106V/m;σ1=σ2=2102V/m; and we will increase the parameter a1 sequentially, assuming a1=0;1106;2106;5106;7106V/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 a1 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 E and each domain is characterized by vector of spontaneous polarization pS and crystallographic axes a,a',c (Figure 23). By simple averaging one can determine the vector of residual polarization and the tensor of strain:


Figure 23.

Determination of angles in the model of orientation switching.

Denote the plane perpendicular to the vector c through B, and the plane passing through the vectors E and c by A. We introduce the following angles: the angle γc between the direction of the field and the axis c; the angle γa between the vector a and the field E; angle ω between the axis a and the line of intersection OK. It is obvious that the angles introduced are within


Let Ecc,Eca be the coercive fields of 180° and 90° switching, respectively. The main conditions of 180° domains switching (turn) are


For 90° domains switching, we have EcosγaEcaEcosγcEcc1, if γcπ/2 and EcosγaEca;EcosγaEcaEcosγcEcc0, if γc>π/2.

A schematic representation of all domains axes 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 Ei. With this goal, the process of loading by an electric field E=Et is replaced by a sequence of values of the electric field EEi:Ei=Eti, and for each state E=Ei, 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 c before and after polarization.

As an example [2], a hysteresis loop (Figure 25) was calculated for the next set of parameters: N=1273248;Ecc=2106V/m;Eca=3106V/m;Emax=6106V/m, where Emax—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: α=1290;P01=0.18104ps;P02=0.79105ps;P03=0.81ps.

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 a,b,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 a,b,c with respect to a fixed system Ox1x2x3 is determined by three angles φ,ψ,ω. 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]:


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 E is the electric field vector; σ is tensor of mechanical stresses; pS is vector of spontaneous polarization, εS is tensor of spontaneous strain, and Uc 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 ps. 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 E and some “molecular field” proportional to its polarization: Eef=E+αP0, where α 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=psEefV, 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 expU/kT, where T is the absolute temperature, k=1.38071023J/K is the Boltzmann constant. Therefore, the averaging operation yields the maximum possible polarization


In the one-dimensional theory, the direction of the electric field coincides with the direction of the axis 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 of the representative volume in the form of a Langevin distribution:


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


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


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:


where δ=signdE, and k is a positive constant, also subjecting to determination. Passing from the effective to the applied field, we obtain an ordinary differential equation:


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: ps,α,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 α 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:


The total deformation consists of the reversible εe and irreversible ε0 parts: ε=εe+ε0, and to determine the irreversible part, we obtain an ordinary differential equation


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.



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


  1. 1. Hall DA. Review. Nonlinearity in piezoelectric ceramics. Journal of Materials Science. 2001;36:4575-4601
  2. 2. Belokon AV, Skaliukh AS. Mathematical Modeling of Irreversible Processes of Polarization. Moscow:Fizmatlit; 2010. 328 p. (In Russian)
  3. 3. Lynch CS. The effect of uniaxial stress on the electro-mechanical response of 8/65/35 PLZT. Acta Metallurgica et Materialia. 1996;44(10):4137-4148
  4. 4. Zhou D. Experimental Investigation of Non-linear Constitutive Behavior of PZT Piezoceramics. Dissertation, Institut für Materialforschung II Von der Fakultät für Maschinenbau der Universität Karlsruhe; 2003 (TH). 141 p
  5. 5. Kamlah M, Wang Z. A Thermodynamically and Microscopically Motivated Constitutive Model for Piezoceramics. Wissenschaftliche Berichte FZKA 6880, Forschungszentrum Karlsruhe in der Helmgoltz-Gemeinschaft; 2003. 99 p
  6. 6. Di Chen. Experimental Investigation of Strain Distributions and Polarization in Lead-based and Lead-free Ferroelectrics. Dissertation Doktors der Ingenieurwissenschaften, des Karlsruher Instituts für Technologie; 2017. 96 p
  7. 7. Zareian Jahromi SA. Nonlinear Constitutive Modeling of Piezoelectric Materials. A thesis submitted to the faculty of graduate studies in partial fulfilment of the requirements for the degree of doctor of Philosophy. Department of Mechanical and Manufacturing Engineering, Calgary; 2013. 165 p
  8. 8. Zhou D, Kamlah M, Munz D. Effects of bias electric fields on the non-linear ferro-elastic behavior of soft lead zirconate titanate piezoceramics. Journal of the American Ceramic Society. 2005;88(4):867-874
  9. 9. Zhou DY, Kamlah M, Munz D. Uniaxial compressive stress dependence of the high-field dielectric and piezoelectric performance of soft PZT piezoceramics. Journal of Materials Research and Technology. 2004;19(3). pp. 834-842
  10. 10. Zhou D, Kamlah M, Munz D. Effects of uniaxial prestress on the ferroelectric hysteretic response of soft PZT. Journal of the European Ceramic Society. 2005;25:425-432
  11. 11. Rayleigh L. The behaviour of iron and steel under the operation of feeble magnetic forces. Philosophical Magazine. 1887;23:225-245
  12. 12. Chen PJ, Peercy PS. One dimensional dynamic electromechanical constitutive relations of ferroelectric materials. Acta Mechanica. 1979;31:231-241
  13. 13. Chen PJ, Tucker TJ. Determination of the polar equilibrium properties of the ferroelectric ceramic PZT 65/35. Acta Mechanica. 1981;38:209-218
  14. 14. Chen PJ, Montgomery ST. A macroscopic theory for the existence of the hysteresis and butterfly loops in ferroelectricity. Ferroelectrics. 1980;23:199-208
  15. 15. Chen PJ, Tucker TJ. One dimensional polar mechanical and dielectric responses of the ferroelectric ceramic PZT 65/35 due to domain switching. International Journal of Engineering Science. 1981;19:147-158
  16. 16. Bailey PB, Chen PJ. Transient electromechanical responses of ferroelectric ceramics to impulsive electric fields. International Journal of Solids and Structures. 1981;17:471-478
  17. 17. Chen PJ, Madsen MM. One dimensional polar responses of the electrooptic ceramic PLZT 7/65/35 due to domain switching. Acta Mechanica. 1981;41:255-264
  18. 18. Chen PJ. Hysteresis effects in deformable ceramics. In: Maugin GA, editor. The Mechanical Behavior of Electromagnetic Solid Continua (IUTAM). North-Holland, Amsterdam: Elsevier Science Publishers; 1984. pp. 137-143
  19. 19. Chen PJ. Three dimensional dynamic electromechanical constitutive relations of ferroelectric materials. International Journal of Solids and Structures. 1980;31:1059-1067
  20. 20. Filatov OP, Khapaev MM. Averaging of Systems of Differential Inclusions: Textbook. Moscow: Publishing House of Moscow University; 1998. 160 p. (In Russian)
  21. 21. Palmov VA. Oscillations of Elastic-plastic Bodies. Moscow: Nauka; 1976. 328 p. (in Russia)
  22. 22. Visintin А. On hysteresis in elasto-plasticity and in ferromagnetism. International Journal of Non-Linear Mechanics. 2002;37:1283-1298
  23. 23. Preisach F. Uber die magnetische nachwirkung. Zeitschrift fur Physik. 1935;94:277-302
  24. 24. Visintin А. Quasilinear first-order PDEs with hysteresis. Journal of Mathematical Analysis and Applications. 2005;312:401-419
  25. 25. Visintin А. Quasilinear hyperbolic equations with hysteresis. Annales de l'Institut Henri Poincaré–AN. 2002;19(4):451-476
  26. 26. Robert G, Damjanovic D, Setter N. Preisach distribution function approach to piezoelectric nonlinearity and hysteresis. Journal of Applied Physics. 2001;90(5):2459-2464
  27. 27. Robert G, Damjanovic D, Setter N. Preisach modelling of ferroelectric pinched loops. Applied Physics Letters. 2000;77(26):4413-4415
  28. 28. Robert G, Damjanovic D, Setter N, Turik AV. Preisach modeling of piezoelectric nonlinearity in ferroelectric ceramics. Journal of Applied Physics. 2001;89(9):5067-5074
  29. 29. Turik AV. To the theory of polarization and hysteresis of ferroelectrics. Solid State Physics. 1963;5(4):1213-1215
  30. 30. Zs S, Tugyi I, Gy K, Fuzi J. Identification procedures for scalar Preisach model. Physica, B. 2004;343:142-147
  31. 31. Azzerboni B, Cardelli E, Finocchio G. A comparative study of Preisach scalar hysteresis models. Physica B. 2004;343:164-170
  32. 32. Ragusa C. An analytical method for the identification of the Preisach distribution function. Journal of Magnetism and Magnetic Materials. 2003;254-255:259-261
  33. 33. Fuzi J. Strong coupling in electromechanical computation. Journal of Magnetism and Magnetic Materials. 2000;215-216:746-748
  34. 34. Kadar G. On the product Preisach model of hysteresis. Physica B. 2000;275:40-44
  35. 35. Hu H, Mrad RB. On the classical Preisach model for hysteresis in piezoceramic actuators. Mechatronics. 2003;13:85-94
  36. 36. Liu J, Sh Z, Chen F, Ch Y, Zhou X. Measurements and simulation of hysteresis loops of donor-doped strontium bismuth tantalate ceramics. Physics Letters A. 2004;321:199-204
  37. 37. Meyer V, Sallese J-M, Fazan P, Bard D, Pecheux F. Modeling the polarization in ferroelectric materials: a novel analytical approach. Solid-State Electronics. 2003;47:1479-1486
  38. 38. Yu Y, Naganathan N, Dukkipati R. Preisach modelling of hysteresis for piezoceramic actuator system. Mechanism and machin theory. 2002;37:49-59
  39. 39. Chuntao L, Yonghong T. A neural networks model for hysteresis nonlinearity. Sensors and Actuators A. 2004;112:49-54
  40. 40. Sjostrom M. Differentiation and power loss computation of classical Preisach model. Physica B. 2004;343:96-100
  41. 41. Pasco Y, Berry A. Consideration of piezoceramic actuator nonlinearity in the active isolation of deterministic vibration. Journal of Sound and Vibration. 2006;289:481-508
  42. 42. Kadar G, Szabo G. Hysteresis modeling. Journal of Magnetism and Magnetic Materials. 2000;215-216:592-596
  43. 43. Fuzi J. Experimental verification of a dynamic hysteresis model. Physica B. 2004;343:80-84
  44. 44. Zirka SE, Moroz YI, Marketos P, Moses AJ. Properties of dynamic Preisach models. Physica B. 2004;343:85-89
  45. 45. Morentin FJ, Alejos O, Francisco C, Munoz JM, Hernandez-Gomez P, Torres C. Simple standard problem for the Preisach moving model. Physica B. 2004;343:107-111
  46. 46. Cardelli E, Bertoncini F, Di Fraia S, Tellini B. Implementation of the modified Preisach scalar model in the finite difference - time - domain numerical modeling. Physica B. 2001;306:126-131
  47. 47. Belbas SA, Mayergoyz ID. Optimal control of dynamical systems with Preisach hysteresis. International Journal of Non-Linear Mechanics. 2002;37:1351-1361
  48. 48. Cross R, Krasnoselskii AM, Pokrovskii AV. A time-dependent Preisach model. Physica B. 2001;306:206-210
  49. 49. Visintin A. Vector Preisach model and Maxwell’s equations. Physica B. 2001;306:21-25
  50. 50. d'Aquino M, Serpico C. A new Preisach-type vector model of hysteresis. Journal of Magnetism and Magnetic Materials. 2004;272-276:731-733
  51. 51. Fuzi J. Two Preisach type vector hysteresis models. Physica B. 2004;343:159-163
  52. 52. Kedous-Lebouc A, Vernescu C, Cornut B. A two-dimensional Preisach particle for vectorial hysteresis modeling. Journal of Magnetism and Magnetic Materials. 2003;254-255:321-323
  53. 53. Fuzi J. Anisotropic vector Preisach particle. Journal of Magnetism and Magnetic Materials. 2000;215-216:597-600
  54. 54. Kahler G, Torre ED. Application of simplified vector Preisach model to vector magnetizing process. Physica B. 2000;275:114-119
  55. 55. Torre ED, Pinzaglia E, Cardelli E. Vector modeling-Part I: Generalized hysteresis model. Physica B. 2006;372:111-114
  56. 56. Serpico C, d’Aquino M, Visone C, Davino D. A new class of Preisach-type isotropic vector model of hysteresis. Physica B. 2004;343:117-120
  57. 57. Fesenko EG, Danziger AYa, Kramarov OP, and Others. Polarization of Ceramics–Rostov on Don, Publishing House of Rostov State University, 1968. 136 p. (in Russia)
  58. 58. Skaliukh A. About hysteretic operators arising in the modelling of the polarization of ferroelectric materials. Proceedings of the International Conference on Mathematical Sciences, July 17-19, 2014, Chennai, India, Sathyabama University, Published By Elsevier, a division of Reed Elsevier India Private Limited, New Delhi; 2014. pp. 635-638
  59. 59. Shindo Y, Narita F, Horiguchi K, Magara Y, Yoshida M. Electric fracture and polarization switching properties of piezoelectric ceramic PZT studied by the modified small punch test. Acta Materialia. 2003;51:4773-4782
  60. 60. Tamm IE. Fundamentals of the Theory of Electricity: Textbook for high schools. 11th ed. Moscow: FIZMATLIT; 2003. 616 p. (in Russia)
  61. 61. Jiles DC, Atherton DL. Theory of magnetic hysteresis. Journal of Magnetism and Magnetic Materials. 1986;61:48-60
  62. 62. Sadowski N, Batistela NJ, Bastos JPA, Lajoie-Mazenc M. An inverse Jiles - Atherton model to take into account hysteresis in time-stepping finite-element calculations. IEEE Transactions on Magnetics. 2002;38(2):797-800
  63. 63. Jiles DC. Introduction to Magnetism and Magnetic Materials. New York: Chapman and Hall; 1991
  64. 64. Smith RC. A nonlinear optimal control method for magnetostrictive actuators. Journal of Intelligent Material Systems and Structures. 1998;9(6):468-486
  65. 65. Smith RC, Hom CL. A domain wall model for ferroelectric hysteresis. SPIE Conference on Mathematics and control in Smart Structures, SPIE Vol. 3667, Newport Beach, CA, March 1-4, 1999. pp. 150-161
  66. 66. Smith RC, Hom CL. Domain wall theory for ferroelectric hysteresis. Journal of Intelligent Material Systems and Structures. 1999;10(3):195-213
  67. 67. Smith RC, Ounaies Z. A domain wall model for hysteresis in piezoelectric materials. Journal of Intelligent Material Systems and Structures. 2000;11(1):62-79
  68. 68. Skaliukh A, Belokon A. Modeling strain and dielectric hysteretic type dependences in polycrystalline ferroelectrics by methods of two-level continuum. 1st International Conference on Rheology and Modeling of Materials (IC-RMM1). Journal of Physics: IOP Conference Series. 2015;602:012025
  69. 69. Parinov IA, editor. Piezoceramic Materials and Devices. Chapter 2. N.-Y.: NOVA Publishers; 2012. pp. 55-104
  70. 70. Benabou A, Clenet S, Piriou F. Comparison of Preisach and Jiles–Atherton models to take into account hysteresis phenomenon for finite element analysis. Journal of Magnetism and Magnetic Materials. 2003;261:139-160
  71. 71. Suzuki T, Matsumoto E. Comparison of Jiles–Atherton and Preisach models extended to stress dependence in magnetoelastic behaviors of a ferromagnetic material. Journal of Materials Processing Technology. 2005;161:141-145

Written By

Alexander Skaliukh

Submitted: 12 January 2018 Reviewed: 02 May 2018 Published: 03 October 2018