Modeling of a hysteretic deformation response in polycrystalline ferroelastics

You In the absence of an electric field a mathematical model describing the ferroelastic response of complete ferroelectrics ferroelastics on action of mechanical stresses is proposed. The modeling is based on the concept of a “ferroelastic” element, similar to the theory of plasticity where used the Saint-Venant element of “dry friction”. The constitutive relations for elastic and residual strains are constructed. The dependence of elastic compliance on the main values of the tensor of residual strains is established. For residual strains, the constitutive relations are obtained in differentials. The obtained constitutive equations can be used in finite element analysis of irreversible processes of deformation of polycrystalline ferroelastics. A number of numerical experiments were performed, which showed good agreement with the experimental data.


Introduction
Modeling of nonlinear deformation effects in ferroelastics is especially relevant for microminiature work elements and thin-film structures, where small efforts can lead to large mechanical stresses that can change the structure of the material due to changes in residual strain. In this paper, we will consider the ferroelastic properties of complete polycrystalline ferroelectrics -ferroelastics of the perovskite structure. Deformations in such materials substantially depend on the intensity of external stresses: under small influences, the deformations are small and the processes are reversible; at loads above threshold values, the process is irreversible and causes plastic deformations. Due to the similarity of the phenomena of magnetization, plasticity and polarization, mathematical models for magnets, plastic bodies and ferroelectrics have much in common. Nonlinear effects in polycrystalline ferroelectric materials are described in detail in [1]. Onedimensional models were numerically investigated in [2]. These include the Rayleigh model [3], models of the theory of plasticity, the Preisach model [4]. The Rayleigh model is valid for stress of medium intensity and describes axial deformation by a quadratic function tensile -compressive stress. In plasticity models, rheological formulas for elastic and plastic strains are conveniently described using differential inclusions [5]. In the Preisach model, the concept of "hysterone" is used [6][7], where the switching of each domain is described by a rectangular hysteresis loop with its values of the coercive ( c  ) and internal ( i  ) stress fields, as in the Mises farm model. Three-dimensional models are constructed mainly on the basis of a generalization of one-dimensional models. They include the construction of constitutive relations for both elastic and residual strains. Elastic components are state parameters; therefore, for them, the constitutive relations are linear algebraic operators, with elastic coefficients depending on the residual strain. Residual parameters are process parameters; therefore, for them, the constitutive relations are described by hysteresis type differential operators. The second part is the most difficult in modeling irreversible processes. At this stage, it is necessary to take into account the structure of the material, formulate the criteria for switching domains, carry out averaging, take into account the influence of neighboring domains on their switching, etc. A more detailed review of three-dimensional models can be found in [8].
In this paper, we consider the mathematical model of polycrystalline ferroelastics using the example of model a locked wall. This choice of model is due to the fact that it is in good agreement with the results of the experiment and being implanted into the finite element complex allows you to store the least amount of information.

The subject of the study
The subject of the study is complete ferroelectrics -ferroelastics of the perovskite type, from which polycrystalline ferroelastic materials or ceramics are made. By their nature, such materials contain a huge number of crystallites, in which there is also a large number of domains, the position of which in the initial state is arbitrary. For illustration the Figure 1 shows the structure of ceramics at different levels: sample (), particle ( * ), crystallite ( g ), unit cell  cr ).

Fig. 1.
Material structure:  -body volume;  * -representative volume (particle of level U1);  g is crystallite;  cr is a particle of level U 2 .
Ferroelectric ceramics [2] consists of a huge number of crystallites with sizes of the order of 10 -6  10 -5 m. On the linear crystallite size, an average of one to hundreds of domains is located, each of which is capable of containing hundreds or more atomic cells. Therefore, in both fine-grained and coarse-grained ceramics, the number of atomic cells on the linear size of crystallite is about 10 4 . Upon transition from the high-temperature to the lowtemperature phase, the unit cell deforms, passing from the cubic to the tetragonal phase, resulting in spontaneous strain fig. 2. At a constant temperature in the material, it is possible to change the orientation of the axes of spontaneous strain by applying external mechanical stress, for example, as shown in Fig. 3. To assess how mechanical stresses affect the structure of the material, we use the following technique. Consider a particle  * , in which unit cells have arbitrary crystallographic axes. We take some reduction point in space and associate with each axial deformation a unit length segment, the middle of which is at the selected point, as shown in Fig. 4. Under the action of tensile stresses, the deformations of all cells fall into the region of two cones, and under the action of compressive stresses, these deformations are distributed in a torus with a sectorial section. The irreversibility of the process lies in the fact that after stress relieving, all spontaneous deformations will remain in the switched position. In this regard, the structure of the material changes; we have a phase transition in the representative volume by "solid" -"solid" type . In other words, mechanical stresses can switch the orientation states of the crystal, change the anisotropy, or drive the position of domain walls. This property for transparent ceramics is used in optics [9,10]. In a representative volume containing N unit cells, we define the residual strain tensor by the averaging operation: Let an external mechanical stress  (determining parameter) be applied to a representative volume of a polycrystalline ferroelectric material, which we assume is the same in the entire representative volume. Stress of low intensity only deforms the walls of unit cells, which leads to the appearance in the representative volume of a reversible component of the strain tensor  e . Mechanical stress above a threshold value cause rotations of the axes of the tensor of spontaneous deformation, which leads to the appearance of irreversible deformation  0 . The total strain tensor  (the desired parameter) is subdivided into the reversible  e and irreversible  0 parts: 3 We will use the elements of thermodynamics of irreversible processes [2]. In the absence of residual deformations, the environment is an isotropic body, and in their presence -anisotropic. Therefore, the elastic characteristics change with a change in the residual strain  0 .

Constitutive relations for reversible deformations
where -mass density, mass density of internal energy, stress tensor, strain tensor, heat flux vector, power of internal heat sources, mass density of entropy, absolute temperature, Hamilton operator, respectively. The dot above denotes the substantial derivative with respect to time, which coincides with the partial derivative due to the validity of the hypothesis of geometrically linear mechanics. Internal energy is a function of external parameters, and entropy, i.e.
In it, for the case of the process of cold deformation, the thermal terms ( T s, ) were neglected. The cell switching dynamics is less than 10 -11 seconds. The process of switching a domain has a time scale of the order of 10 -8  10 -5 seconds. Therefore, when the stress varies within 10 -1  10 -2 seconds, we have a quasistatic process, and small deviations from thermodynamic equilibrium. It follows that the dissipation rate function i S can be considered linear with respect to the velocities of its parameters, i.e. , ε σ , as well as on the velocities of these parameters 0 , ε σ   , and they should ensure the implementation of dissipation rate inequality (2). For ceramic materials, the reversible parts of the strain e ε are quite small (approximately an order of magnitude smaller than the residual), which is why it is natural to construct the constitutive relations in the form of linear relations. To this end, we specify the Gibbs function by selecting it in a quadratic form: . Using the independence of the generalized velocities 0 ε σ  , , and having carried out standard mathematical operations in (2), on the basis of (4) we obtain a pair of relations. The first of them relates the strain tensor to the tensor of mechanical stresses. And the second equation relates the tensor of mechanical stresses and the derivatives of the compliance modules with respect to the residual strain. Avoiding excessive generalization, we neglect the viscoelastic properties, for which we put . Then the first ratio can be represented as: The second ratio can be written as follows: Considering e ε ε ε   0 , we can say that (5) are the constitutive relations for reversible parameters. It remains only to find out the nature of the dependence of physical characteristics on residual strains. To this end, we consider the generalized dissipative forces S 0 ε χ included in (6). For irreversible processes, we are not entitled to assume that they are analytical functions of their arguments. We assume that they consist of analytical and non-analytical parts, with the obligatory fulfillment of inequality (2). We denote the analytic parts by , and the non-analytic ones Obviously, the analytical parts cannot depend on the residual parameters, so they will be functions of only external parameters: ) ( It should also be noted that in the absence of external loads, these functions are equal to zero. All this allows us to decompose them into a Taylor series in the vicinity of an unloaded state. Performing such an expansion, up to quadratic components, we obtain: σ K σ σ A χ : : : Decomposition tensors 6 4 T ; T independent of the residual strain and stresses. The introduced analytical parts of dissipative forces do not change the dissipation inequality in any way (2). Substituting these expansions in (6) and taking into account the independence of the determining parameters, we obtain a system of equalities that connect these parameters with the non analytic parts of dissipative forces ) ( 0 na S  χ and with derivatives of physical characteristics, i.e. we have: Relations (7) describe the change in irreversible parameters. However, the construction of the non analytic parts of the generalized dissipative forces is a very difficult task, which, apparently, is not solved even in the simplest cases. Instead, other methods and approaches are used to find irreversible parameters [8], with the condition of satisfying inequality (2). However, for our purposes, relations (8) are important, which allow us to formulate a linear dependence of the tensor ) ( 0 ε S on the residual parameters up to constant tensor 0 S : Obviously, 0 S is the elastic compliance tensors for ceramics of the initial state, when 0 0  ε . The elastic compliance tensor 0 S in this case must be determined before the appearance of an irreversible deformation process. It remains only to specify the components of the tensor K , which is no longer possible to determine using the methods of thermodynamics; here, as in determining the physical characteristics of the material, experimental methods are necessary. Note that this tensor, being connected with physical modules, must possess the same symmetry properties as the compliance tensor S . In addition, it must be symmetric in the last pair of indices, which follows from (8). This circumstance allows us to use the Voigt matrix representation instead of the tensor representation.
To determine the K matrix coefficients, we also use the well-known values of elastic compliance sa t Ŝ at maximum plastic tensile strains along the Oz axis, generating transversal isotropy in the coordinate system with the Oz axis directed along the main axis of the strain tensor:           III  0  II  0  I  0  III  0  3  II  0  2  I  0  1 , where 0  is some factor whose dimension coincides with the dimension of mechanical stresses. Then . From here we determine  m , after which we find: In these relations 0 , since the incompressibility conditions are satisfied.
Because of this, the indefinite factor leaves 0  , and the expression takes the form: Note that only the largest principal value of the residual strain tensor is included in (12), and all matrix components on the right-hand side of (10), (11) are completely known.

Derivation of the ultimate dependence
By an ideal case, we mean a polycrystalline ferroelastic material in which there are no mechanisms for blocking the motion of domain walls. All unit cells in the field of intense stresses come into motion, which obeys certain statistical laws. The main axes of spontaneous strains of each cell in an undeformed state are oriented in a chaotic manner. Under the influence of mechanical stresses, spontaneous strains of each cell change. Therefore, the strain of the particle of continuum 1 U is found by averaging all spontaneous strains of the particles of the continuum level 2 U . Let  d is the volume of the particle of continuum 1 U , and  d let be the volume of the particle by continuum ). By averaging we mean simple summation, which, due to the huge number of cells, can be replaced by integration, i.e. where ξ1, ξ2 and ξ3 are local coordinates in the volume ω; f (ξ1, ξ2, ξ3) is the density function of the distribution of cells in the volume ω; F(ξ1, ξ2, ξ3) is the density function of the distribution of spontaneous strains in the volume ω. The main task of modeling residual strain is to construct either the function f (ξ1, ξ2, ξ3) or the function F(ξ1, ξ2, ξ3). Below, we construct the function F(ξ1, ξ2, ξ3). Let a particle of the level continuum 1 U be characterized by coordinates (x1, x2, x3). It is convenient to characterize the particle of the level continuum 2 U , i.e., the position the domain in this particle, with spherical coordinates (    , , ). It can be seen that domains with different coordinate  , but with the same coordinates  and  , make the same contribution to the total strain of the particle of the continuum level 1 U . This allows averaging not over all three coordinates (    , , ), but only over the last two. Then it suffices to introduce the unit sphere centered at the point (x1, x2, x3). As a result, the direction of the axis of spontaneous stretching deformation of the domain will be characterized by two angles  and  , and, and the averaging operation will be reduced to integration over this sphere.
Domain switching does not affect the stress field. This invariability means that its average value is equal to the stress field itself: σ σ   . In addition, spontaneous deformations do not change their values when the direction of the main axis changes to the opposite; therefore, all averaging operations are sufficient to carry out along the upper semi sphere. We derive the limiting dependence following the approach of [11].
The number of domains located in a particle of the continuum of the level 1 U is equal to the number of segments symbolizing a given deformation and reaching the surface of a unit sphere. In the initial state, the number of segments reaching the surface of a sphere element is proportional to the area of the element: , where c is a certain coefficient. According to Weiss theory, the process of switching domains in a particle of the continuum of the level 1 U is affected not by the field of true 7 stresses σ , but by the "effective" field 0 ε σ σ    ef . According to Boltzmann's theorem in a conservative field ef σ , the distribution of domain axes in a particle differs from their distribution in the absence of this field by the value . Thus, to count all switched domains in a particle, it is sufficient to carry out integration over the semi sphere: (13) To determine the strain of a particle of the continuum of the level 1 U , it is necessary to average the spontaneous strains of all the domains that make up the particle, and exclude N from the resulting ratio, we obtain: The resulting expression determines the most "favorable" case of switching domains in the field and is called the ultimate or a hysteresis one.

Derivation of the ultimate dependence
In real processes, the rotation of the axes of spontaneous strains begins only when the level of mechanical stresses is sufficient to break mechanisms locking domain walls. We'll estimate the energy required for this breakdown. Consider a certain level of effective stresses ef σ for which the cell has not switched.  a  a  a   a  a  a  a  a  a  a  a  a  a  a  a Assuming that such relations are valid for any domain, we carry out averaging at the volume of a particle of the continuum of the level Instead of expression (15), we take average values for the breakdown energy, where the value is determined from the previous relation. For a more general analysis, the average distribution density of the locked walls is introduced and the energy spent on breaking the locked walls in an arbitrary  is determined: where cr s

Estimation the work of mechanical stresses
Work in the process of development of deformation ε d under the action of stresses σ in an arbitrary volume  can be represented by the following dependence: Thus, expression (18) can be represented as a sum of five terms:  is a full differential of the functions of complete strain; 2 A  -the full differential of elastic strains with a changing tensor of elastic modules; 4 A  represents the full differential of plastic strains. If we consider the cyclic process, when the full and plastic strains return to their initial values, then . 0 ; 0 ; 0  represents an additional work of elastic strains due to changes of the elastic modules in an irreversible process. It can be seen that, within the framework of one circular (or cyclic) process, for each state with elastic strains e ε , there always exists a state in which the elastic strains have the opposite sign (e ε ). We call such states opposite. These conditions of symmetry allow us to assume that the increments of the tangent elastic modules in opposite states will have the same numerical values, but opposite signs. Given that the integrand is a quadratic function of elastic strains, it can be argued that increments of the additional work on elastic strains in opposite states will also have the same numerical values and opposite signs. The result for a circular process The full work in a circular process So, the fifth integral is responsible for energy losses in a cyclic process, i.e., they are determined by the expression (18)

Derivation of the energy ratio
If it were not necessary to spend work on destruction the mechanisms of locking the domain walls, then the energy loss in the system would be much less and would be determined by the ideal (or ultimate) case in which the rotation of an individual unit cell by mechanical stress would not depend on the influence and location of the remaining cells. To calculate the losses in the ideal case, it is enough in the previous expression to replace the residual strain tensor with the ultimate strain tensor. Denote such losses: At the same time, in a real process, energy losses in the system are much larger than in the limiting case, since part of the work is spent on the destruction of locking mechanisms. All this allows us to formulate the following condition of the energy balance: real energy losses (18)  . The parameter a appears when determining the ultimate strain (14). Considering the room temperature, the value of the Boltzmann constant, spontaneous strain, and the size of the unit cell volume, we can assume that this parameter varies within 7 10  a . The parameter k is directly proportional to the average energy density needed to rotate the main axes of spontaneous strains through an angle of 90 0 , the average density of the locked walls, and inversely proportional to value of spontaneous strain. Given that the average energy density for breaking the locked mechanisms is approximately equal to the product of the coercive stress and spontaneous strain, and 1 0 ). The influence of the parameter k most affects the shape, the area, and the amplitude of the hysteresis curves. This is due to the fact that it includes the average energy needed to break down the mechanisms locking the walls of domains. The lower this energy, the lower the intensity of the mechanical stress required to switch domains. Therefore, at lower values of the parameter, the loop is narrower. On the other hand, if less energy is spent on break down the mechanisms locking the walls of domains, the more domains will switch, and the greater the amplitude of the loop at the same intensity of mechanical stresses, as shown in the figures.

Conclusion
In this work, the constitutive relations are constructed for polycrystalline ferroelastics at arbitrary loads in the three-dimensional case. They include relationships for elastic and residual strain components. The constitutive relations for the elastic components are linear algebraic relations connecting stresses and strains using the Hooke law, however, the compliance coefficients depend on the residual strains. The nature of these dependences was investigated and it was found that they linearly depend on the principal values of the residual strain tensor. The constitutive relations for residual strains were obtained using elements of the mechanics of a two-level medium and energy estimates of the breakdown of mechanisms that prevent the switching of spontaneous deformation in unit cells, the energy costs for the work of mechanical stresses in real and ideal (ultimate) process of deformation . Equations in the differentials with respect to residual deformations are obtained from the energy balance. The boundaries of the variation of the parameters included in the equations are established, the influence of each of the parameters on the type, area and shape of the strain hysteresis curves is numerically studied. It is shown that with a certain choice of model parameters, it is possible to obtain curves that coincide with the experimental data not only qualitatively, but also quantitatively. Thus, the constructed model can be implanted into a finite element complex to calculate the stress-strain state of polycrystalline ferroelastics. It is noteworthy that the model requires storing a minimum of information that relates only to residual deformations. This work was supported by the Russian Foundation for Basic Research (grant 17-08-00860-a).