Macroscopic constitutive model for ergodic and non-ergodic lead-free relaxors

Abstract

A fully electromechanically coupled, three dimensional phenomenological constitutive model for relaxor ferroelectric materials was developed for the use in a finite-element-method (FEM) solution procedure. This macroscopic model was used to simulate the macroscopic electromechanical response of lead-free ergodic $0.94\mathrm{N}{\mathrm{a}}_{1/2}\mathrm{B}{\mathrm{i}}_{1/2}\mathrm{Ti}{\mathrm{O}}_3-0.06\mathrm{BaTi}{\mathrm{O}}_3$ and non-ergodic $0.90\mathrm{N}{\mathrm{a}}_{1/2}\mathrm{B}{\mathrm{i}}_{1/2}\mathrm{Ti}{\mathrm{O}}_3-0.06\mathrm{BaTi}{\mathrm{O}}_3-0.04{\mathrm{K}}_{0.5}\mathrm{N}{\mathrm{a}}_{0.5}\mathrm{Nb}{\mathrm{O}}_3$ relaxor materials. The presented constitutive model is capable of accounting for the observed pinched hysteretic response as well as non-deviatoric polarization induced strain and internal order transitions. Time integration of the history dependent internal variables is done with a predictor-corrector integration scheme. The adaptability of the constitutive model regarding the pinching of the hystereses is shown. Simulations are compared to experimental observations.

Introduction

In the field of electromechanical transducer applications, lead titanate zirconate (PZT) dominates today’s market due to its excellent, tunable electromechanical properties. Environmental and health concerns related to the lead content of PZT, require nontoxic alternatives with equal or better performance. Over the last two decades, there have been a number of lead-free systems developed that display promising electromechanical properties (Liu and Ren, 2009; Li et al., 2013; Zhang et al., 2007). Among these lead-free alternatives, binary (Sakata and Masuda, 1974; Takenaka et al., 1991) and ternary (Anton et al., 2011; Zhang et al., 2007) solid solutions based on $\mathrm{N}{\mathrm{a}}_{1/2}\mathrm{B}{\mathrm{i}}_{1/2}\mathrm{Ti}{\mathrm{O}}_3$ (NBT) have been widely investigated for use in actuation systems due to their large strain behavior (Jo et al., 2009; Panda, 2009; Zhang et al., 2007). The origins of this behavior are related to the chemically disordered relaxor state of NBT-based materials, where previous studies have demonstrated that an external electric field (Daniels et al., 2009; Ma and Tan, 2011; Tan et al., 2011) or mechanical stress (Garg et al., 2013; Peng et al., 2018; Schader et al., 2016; Webber et al., 2017) can induced a long-range ferroelectric order and a corresponding jump in strain and polarization. Importantly, this effect can be modulated either chemically, through the introduction of additional end members to destabilize the ferroelectric order and thereby reduce the critical temperature, or thermally (Foronda et al., 2014; Jo et al., 2009). This field-induced long-range ferroelectric order, however, leads to problems that need to be addressed before using these materials in applications, including a relatively large critical electric field required to acquire the large strains and a significant hysteresis that indicates large energy loss and potential overheating during usage.

Lee et al. (2011, 2012) demonstrated a solution for both of these challenges by using ceramic-ceramic 0-3 type composites (connectivity nomenclature by Newnham et al. (1978)) of a nonergodic relaxor or ferroelectric phase (seed) and an ergodic relaxor (matrix). Here, the primary difference between the components is in the stability of the ferroelectric phase during the application of an electric field, which significantly affects the macroscopic polarization-electric field response. In an ergodic relaxor, the ferroelectric phase is unstable and forms back to a disordered relaxor state after the application of an electric field, reducing the remanent polarization. In contrast, nonergodic relaxors form an irreversible ferroelectric structure, making the material essentially a ferroelectric after the application of an electric field above the transition field $E^t$ with a corresponding large remanent polarization. Differences in the polarization and strain response of the components cause polarization- and strain-coupling effects, which can enhance the unipolar strain of the composite system. As such, throughout the manuscript ergodic relaxor will be called relaxor (RE) and the nonergodic relaxor, ferroelectric (FE). Recently, different 2-2 composites have been presented in order to better understand the underlying mechanisms that lead to the improved strain response. For instance, the importance of strain coupling in composites consisting of $0.94\mathrm{N}{\mathrm{a}}_{1/2}\mathrm{B}{\mathrm{i}}_{1/2}\mathrm{Ti}{\mathrm{O}}_3-0.06\mathrm{BaTi}{\mathrm{O}}_3$ (NBT-6BT) in ferroelectric (FE) phase and $0.90\mathrm{N}{\mathrm{a}}_{1/2}\mathrm{B}{\mathrm{i}}_{1/2}\mathrm{Ti}{\mathrm{O}}_3-0.06\mathrm{BaTi}{\mathrm{O}}_3-0.04{\mathrm{K}}_{0.5}\mathrm{N}{\mathrm{a}}_{0.5}\mathrm{Nb}{\mathrm{O}}_3$ (NBT-6BT-4KNN) in relaxor (RE) phase was experimentally shown in Martin et al. (2021). Despite previous work, however, there remain a number of important and open questions regarding the influence of variations in electromechanical properties between end members on the internal field distribution in layered composites. Through the development of a suitable constitutive model, the influence of such factors as well as the optimization of the composite structure can be addressed in numerical simulations.

Finite-element based phase-field models for ferroelectric and relaxor ferroelectric multilayer ceramics have been previously reported, where the importance of polarization and strain coupling in NBT-based multilayer composites was investigated (Franzbach, 2014; Wang et al., 2019). Although phase-field based models are suitable for simulating domain dynamics, they are computationally expensive, especially when simulations on the component length scale are required. An alternative is the class of micromechanical models. To date, they have been applied to macroscopic thermo-electromechanical constitutive behavior of ferroelectrics, where the dominate hysteretic process is domain wall motion, as well as field-induced structural phase transitions in relaxor single crystal ferroelectrics (Webber et al., 2008). A review of micromechanically motivated models, which aim mostly at length scales comparable to the grain size, like thin films, is given in Huber (2005). More suitable for structural mechanics analyses is the class of macroscopic phenomenological constitutive models introduced first in Chen and Peercy (1979) in a one-dimensional formulation and afterwards subsequently expanded to multi axial formulations (Elhadrouz et al., 2005; Ghandi and Hagood, 1997; Kamlah and Wang, 2003; Klinkel, 2006; Landis, 2002; Lynch, 1998; Schwaab et al., 2012). To the best of our knowledge, no available model can represent well the (pinched) hysteresis phenomena of NBT-6BT and NBT-6BT-4KNN in the form of a fully coupled electro-mechanical 3D constitutive model. The development of such an advanced constitutive model is the scope of the present publication.

Experimental methodology

Polycrystalline NBT-6BT and NBT-6BT-4KNN samples were synthesized by a mixed-oxide route using the starting powders from Alfa Aesar, (Karlsruhe, Germany): $\mathrm{B}{\mathrm{i}}_2{\mathrm{O}}_3$ (99.975% purity), $\mathrm{Ti}{\mathrm{O}}_2$ (99.6% purity), $\mathrm{N}{\mathrm{a}}_2\mathrm{C}{\mathrm{O}}_3$ (99.95% purity), ${\mathrm{K}}_2\mathrm{C}{\mathrm{O}}_3$ (99.95% purity), $\mathrm{N}{\mathrm{b}}_2{\mathrm{O}}_5$ (99.95% purity) and $\mathrm{BaC}{\mathrm{O}}_3$ (99.95% purity). Stoichiometrically weighed raw materials were ball milled for $24\mathrm{h}$ , before a two-step calcination at $900°\mathrm{C}$ for $3\mathrm{h}$ . Following cooling to room temperature, powder compacts were is ostatically pressed at $-180\mathrm{MPa}$ and sintered at $1150°\mathrm{C}$ for $3\mathrm{h}$ using sacrificial powder of same composition, with a $5\mathrm{K}/\min$ heating rate. The resulting samples had a dense and homogeneous microstructure, with a relative density $\ge 96\%$ and an average grain size of approximately $0.5\mu \mathrm{m}$ .

Prior to electrical measurements, samples were cut and ground to final dimensions of $4\mathrm{mm}\times 4\mathrm{mm}\times 4\mathrm{mm}$ , where Au electrodes were sputtered on opposing $4\mathrm{mm}\times 4\mathrm{mm}$ faces. The macroscopic ferroelectric behavior was characterized with a custom-built digital image correlation (DIC) system outfitted with an additional a linear variable differential transformer (LVDT) and a Sawyer-Tower circuit. A bipolar electric field of $\pm 4\mathrm{kV}/\mathrm{mm}$ was applied with a rate of $1.6\mathrm{kV}/(\mathrm{mm}\mathrm{s})$ or a frequency of $100\mathrm{mHz}$ at room temperature with a high voltage amplifier (Trek 20/20). To prevent arcing during the electrical loading, the samples have been immersed in silicon oil. The DIC system was equipped with a $2\times$ magnification lens (MVO-TML telecentric measuring lens, Edmund Optics Inc.), providing a resolution of $1.75\mu \mathrm{m}/\mathrm{pixel}$ . Artificial speckles were sprayed onto the observation surfaces by means of an air brush (AT-Airbrush Pistole Kit, AT-AK-02, Agora-Tec), creating patterns suitable for DIC tracking. During testing, the camera captured two images per second to record the displacement field during electrical loading. The images were analyzed with the commercially available DIC-program Veddac. Here, measurement points had a distance of $60\mathrm{pixel}$ or $105\mu \mathrm{m}$ to each other. In order to receive the displacement of the sample, a reference field of $120\times 120\mathrm{pixel}$ was chosen. Additionally, the hair wavelet function was chosen to increase accuracy and minimize the standard deviation in-between points. The total lateral and transverse strain was calculated by averaging the strain across the entire observation field. Noise level was determined to be $0.04\mathrm{pixel}$ or approximately 0.002 %, well below the measurement signal. Mechanical measurements were performed with cylindrical samples with dimensions of $6\mathrm{mm}$ in height and $5.8\mathrm{mm}$ in diameter. The stress-strain behavior of unpoled NBT-6BT and NBT-6BT-4KNN was measured during uniaxial compressive stress loading with an experimental setup previously described detailed in Webber et al. (2009). A preload of $-3.8\mathrm{MPa}$ was used to ensure mechanical contact, followed by an increase in the compressive stress to $-500\mathrm{MPa}$ and subsequent unloading. Both loading and unloading were done with a linear loading/unloading rate of $5\mathrm{MPa}/\mathrm{s}$ or a frequency of $1\mathrm{mHz}$ .

The macroscopic ferroelectric and ferroelastic response of NBT-6BT and NBT-6BT-4KNN is shown in Figure 1. These results correspond well to previous reports that show a reduction in the remanent polarization and strain with increase KNN content in NBT-6BT (Dittmer et al., 2012). This is understood to be due to the destabilization of the electrically induced long-range ferroelectric order. An important fingerprint of this response is the apparent pinching in the macroscopic polarization-electric field behavior. Similarly, stress-strain behavior also shows a clear reduction in the remanent strain for NBT-6BT-4KNN. In situ structural investigations have demonstrated an analogous mechanical behavior, where externally applied stress can induce ferroelectric ordering and the formation of domain structure (Martin et al., 2019).

Figure 1.

Measured $P_3-E_3$ , $S_3-E_3$ (blue) and $S_1-E_3$ (red) as well as $T_3-S_3$ hystereses of NBT-6BT (left) and NBT-6BT-4KNN (right).

Basics

In this publication, vectors, and tensors of second, third and fourth order are denoted with arrows, lower case bolt, lower case Fraktur and upper case blackboard bold symbols, respectively if not the index notation is used. The Voigt notation is used for the symmetric strain tensor $\mathit{S}$ , that is, $S_1=S_{11}$ , $S_2=S_{22}$ , $S_3=S_{33}$ , $S_4=2S_{23}$ , $S_5=2S_{13}$ and $S_6=2S_{12}$ , and for the stress tensor $\mathit{T}$ , that is, $T_1=T_{11}$ , $T_2=T_{22}$ , $T_3=T_{33}$ , $T_4=T_{23}$ , $T_5=T_{13}$ , and $T_6=T_{12}$ . For the third order piezoelectricity tensor $\mathfrak{d}$ , this yields $d_{33}=d_{333}$ , $d_{31}=d_{311}$ , $d_{15}=2d_{131}$ .

Quasi-electrostatic Maxwell equations

Due to quasi-electrostatic loading, magnetic effects are neglected in this theory. The fundamental Maxwell equations then reduce to the Gauss’s law in terms of the dielectric displacement $\overrightarrow{D}$ with the external free charge density $q^F$

$$\mathrm{div}(\overrightarrow{D})=q^F$$

(1)

and the Faraday’s law

$$\mathrm{rot}(\overrightarrow{E})=0$$

(2)

which states that the electric field $\overrightarrow{E}$ is free of rotations. This is equivalent to the statement that the electric field $\overrightarrow{E}$ can be derived in dependence of the electric potential $\phi$ to

$$\overrightarrow{E}=-\mathrm{grad}(\phi ).$$

(3)

Further details on the electrostatic assumption can be found in Maugin (1988). No electric conductivity is taken into account in the presented model.

Quasi-static mechanics equations

The mechanical part of this theory aims for ceramic materials under quasi-static loading conditions. Hence the small deformation theory is sufficient. The infinitesimal strain tensor

$$\mathit{S}=\frac{1}{2}(\mathrm{grad}\left(\overrightarrow{u}\right)+\mathrm{grad}{\left(\overrightarrow{u}\right)}^{\mathrm{T}})$$

(4)

is defined as the symmetric gradient of the displacement vector $\overrightarrow{u}$ . The superscript ${•}^{\mathrm{T}}$ denotes the transposition. From the balance of linear momentum, assuming quasi-static accelerations, the force equilibrium

$$\mathrm{div}\left(\mathit{T}\right)+{\overrightarrow{f}}^B=0$$

(5)

follows. Therein, $\mathit{T}$ and ${\overrightarrow{f}}^B$ are the Cauchy stress tensor and the body force. From the balance of angular momentum follows the symmetry of the stress tensor $\mathit{T}={\mathit{T}}^{\mathrm{T}}$ , which is already accounted for in the Voigt notation. Further details on continuum mechanics can be found for example in Holzapfel (2008).

Mixed finite-element formulation

For the implementation of the nonlinear constitutive model for non-conducting electromechanical continua, the mixed finite-element formulation

$${\int }_V\mathit{T}:\delta \mathit{S}\mathrm{d}V-{\int }_V{\overrightarrow{f}}^B·\delta \overrightarrow{u}\mathrm{d}V-{\int }_{S_f}{\overrightarrow{f}}^{S_f}·\delta {\overrightarrow{u}}^{S_f}\mathrm{d}{S}_f=0$$

(6)

$${\int }_V(\overrightarrow{E}+\mathrm{grad}\left(\phi \right))·\delta \overrightarrow{D}\mathrm{d}V=0$$

(7)

$$\begin{array}{c}{\int }_V\overrightarrow{D}·\mathrm{grad}\left(\delta \phi \right)\mathrm{d}V+{\int }_V{q}^F·\delta \phi \mathrm{d}V\hfill \\ +{\int }_{S_q}{q}^{S_q}·\delta {\phi }^{S_q}\mathrm{d}{S}_q=0\hfill \end{array}$$

(8)

according to Ghandi and Hagood (1997); Harper (1999); Schwaab et al. (2012) is used. This formulation is faster and more stable than the conventional element formulation for electromechanical problems, which only uses the displacement vector $\overrightarrow{u}$ and the electric potential $\phi$ as independent variable. The mixed formulation employed here uses the dielectric displacement $\overrightarrow{D}$ as an additional independent variable. $V$ indicates the volume, $S_f$ and $S_q$ the surfaces on which the surface traction ${\overrightarrow{f}}^{S_f}$ and the surface charge $q^{S_q}$ are acting. $\delta \overrightarrow{u}$ , $\delta \phi$ , $\delta \overrightarrow{D}$ are the related test functions. $\delta \mathit{S}$ is the symmetric gradient of $\delta \overrightarrow{u}$ , analog to equation (4). $\delta {\overrightarrow{u}}^{S_f}$ and $\delta {\phi }^{S_q}$ are the related test functions on the surfaces. The electric field $\overrightarrow{E}$ in equation (7) will be expressed by a constitutive equation. Consequently $\overrightarrow{E}$ is only equivalent to the definition in equation (3) in the weak average. To not confuse both definitions, the latter is redefined to

$${\overrightarrow{E}}_{\mathrm{grad}}=-\mathrm{grad}(\phi ).$$

(9)

Macroscopic constitutive model

Differences in the introduction mentioned phenomenological constitive models are the choice of the internal variables, which are needed to describe the piezoelectric hysteresis phenomena, how saturation of them is modeled, the time integration methods, and the required finite-element formulation. All classes of constitutive models for ferroelectric materials have their individual advantages and drawbacks. A further discussion on ferroelectric and ferroelastic constitutive models is given in the review article (Kamlah, 2001). In this work, we focus on a fully coupled, electromechanical multi axial phenomenological model that is suitable for finite-element implementation. For the representation of the dielectric, ferroelectric and ferroelastic hysteresis, it is important that the evolution of the polarization and the mechanical irreversible strain are described in the model. With such a model, full size specimens and later on multilayer composites can be simulated under various loading conditions. The above mentioned ferroelectric constitutive models are not able to capture well the behavior of the two materials investigated here accurately. Especially the order transition can’t be represented. Therefore, three extensions to the model from Schwaab et al. (2012) are made: First, the ability to account for pinched hysteresis phenomena. Second, account for an internal order transition phenomenon with primary and secondary load paths and according strain behavior. Third, account for non-deviatoric polarization induced strains. With these modifications, the extended constitutive model will be able to represent the behavior of NBT-6BT and NBT-6BT-4KNN although they have different characteristics.

Time continuous phenomenological constitutive model

The constitutive model to be developed here is an extension of the model in Schwaab et al. (2012). The ferroelectric material behavior is therein modeled with two switching and two saturation functions. In contrast to piezoelectricity, the ferroelectric material behavior is characterized by the load history dependence of internal variables. The main task of a ferroelectric model is to provide the evolution of these internal variables. The (total) polarization $\overrightarrow{P}$ and the (total) strain $\mathit{S}$ are decomposed additive in reversible (index r) and irreversible (index i) parts

$$\overrightarrow{P}={\overrightarrow{P}}^r+{\overrightarrow{P}}^i,$$

(10)

$$\mathit{S}={\mathit{S}}^r+{\mathit{S}}^i.$$

(11)

Once the irreversible quantities are known, the constitutive piezoelectricity equations

$$\overrightarrow{E}=-{\mathfrak{h}}^{\mathrm{T}}:(\mathit{S}-{\mathit{S}}^i)+{\beta }^S·(\overrightarrow{D}-{\overrightarrow{P}}^i)$$

(12)

$$\mathit{T}={\mathbb{C}}^D:(\mathit{S}-{\mathit{S}}^i)-\mathfrak{h}·(\overrightarrow{D}-{\overrightarrow{P}}^i)$$

(13)

are evaluated by reducing the overall strain and dielectric displacement by the irreversible strain and polarization parts. In this work, the piezoelectric equations are denoted in the so called h-form. Details and various other forms of piezoelectric constitutive equations can be found in the textbook Ikeda (1990). Equations (12) and (13) describe the material answer of electric field $\overrightarrow{E}$ and mechanical stress $\mathit{T}$ as functions of the reversible quantities $(\mathit{S}-{\mathit{S}}^i)$ and $(\overrightarrow{D}-{\overrightarrow{P}}^i)$ with the second order impermittivity tensor ${\beta }^S$ , the third order piezoelectric tensor $\mathfrak{h}$ , and the fourth order tensor of elasticity under constant dielectric displacement ${\mathbb{C}}^D$ .

The irreversible strain is composed of two components in the form

$${\mathit{S}}^i={\mathit{S}}^{\mathrm{ip}}+{\mathit{S}}^{\mathrm{im}}.$$

(14)

Here, ${\mathit{S}}^{\mathrm{ip}}$ is induced by irreversible polarization and can be described as a function of ${\overrightarrow{P}}^i$ , while ${\mathit{S}}^{\mathrm{im}}$ represents the irreversible strains, that is, the ferroelastic strains induced by mechanical stresses. Because of this, ${\overrightarrow{P}}^i$ and ${\mathit{S}}^{\mathrm{im}}$ can be taken equivalently as the internal variables of the model which are needed to describe the hysteresis behavior. To decrease the complexity of the model, simplifications and assumptions are made.

Ferroelectric ceramics show a history dependent anisotropy, in the sense that in the pristine unpoled state, they are macroscopically isotropic, while in the electrically poled state, they are of transversal anisotropy. As a matter of fact, this will be represented in our constitutive model in the sense, that, through its dependence on irreversible polarization, the piezoelectricity tensor (dd) represents a history dependent transversal anisotropy. In principle, the history dependent anisotropy of the tensor of elasticity ${\mathbb{C}}^E$ and the tensor of dielectric permittivity ${\mathbf{\kappa }}^T$ could be included as well by representing them as functions of appropriate internal variables such as irreversible strain. This would make the constitutive model and its numerical solution much more complex. Since elasticity and dielectric behavior are only modified by poling and do not completely disappear in the unpoled state, we choose, for simplicity, tensors ${\mathbb{C}}^E$ and ${\mathbf{\kappa }}^T$ to be isotropic constants. Additional, no frequency dependence is taken into account. The fourth order tensor of elasticity at constant electric field ${\mathbb{C}}^E$ is therefore expressed as ${\mathbb{C}}^E=\lambda \mathit{I}\otimes \mathit{I}+2\mu {\mathbb{I}}^{\mathrm{sym}}$ with the Lamé parameters $\lambda =\frac{Y\nu }{(1+\nu )(1-2\nu )}$ and $\mu =\frac{Y}{2(1+\nu )}$ which are functions of the Young’s modulus $Y$ and the Poisson’s ratio $\nu$ . ${\mathbb{I}}^{\mathrm{sym}}=\frac{1}{2}({\delta }_{\mathrm{ik}}{\delta }_{\mathrm{jl}}+{\delta }_{\mathrm{il}}{\delta }_{\mathrm{jk}})$ is the symmetric part of the fourth order unit tensor. ⊗ denotes the dyadic product, · the point product and : the double contraction. The second order tensor of dielectric permittivity at constant mechanical stress ${\mathbf{\kappa }}^T={\kappa }^T\mathit{I}$ is expressed by the scalar quantity ${\kappa }^T$ as well as the second order unit tensor $\mathit{I}={\delta }_{\mathrm{ij}}$ .

The components of the transposed third order piezoelectric tensor ${\mathfrak{d}}^{\mathrm{T}}$ , which describes a linear dependence between electric field and strain, are defined as a function of the irreversible polarization according to Kamlah (2001) to

$$\begin{array}{c}\multicolumn{1}{c}{d_{\mathrm{kij}}=\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}[d_{\parallel }\left(e_i^{P^i}{e}_j^{P^i}{e}_k^{P^i}\right)+d_{\perp }({\delta }_{\mathrm{ij}}-e_i^{P^i}{e}_j^{P^i})e_k^{P^i}}\\ +d_{=}\frac{1}{2}(({\delta }_{\mathrm{ki}}-e_k^{P^i}{e}_i^{P^i})e_j^{P^i}+({\delta }_{\mathrm{kj}}-e_k^{P^i}{e}_j^{P^i})e_i^{P^i})]\hfill \end{array}$$

(15)

with ${\overrightarrow{e}}^{P^i}=\frac{{\overrightarrow{P}}^i}{\Vert {\overrightarrow{P}}^i\Vert }$ and the three piezoelectric coefficients $d_{\parallel }$ , $d_{\perp }$ and $d_{=}$ which correspond in a state fully poled in 3-direction to $d_{33}$ , $d_{31}$ , and $d_{15}$ , respectively. The symmetry $d_{\mathrm{kij}}=d_{\mathrm{kji}}$ holds. In short, the material tensors are computed by the following equations

$$\begin{array}{c}{ϵ}^T={ϵ}_0\mathit{I}+{\mathbf{\kappa }}^T{ϵ}^S={ϵ}^T-{\mathfrak{d}}^T:{\mathbb{C}}^E:\mathfrak{d}\hfill \\ {\beta }^S=({ϵ}^S{)}^{-1}{\mathfrak{e}}^T={\mathfrak{d}}^T:{\mathbb{C}}^E\hfill \\ {\mathfrak{h}}^T={\beta }^S·{\mathfrak{e}}^T{\mathbb{C}}^D={\mathbb{C}}^E+\mathfrak{h}·{ϵ}^S·{\mathfrak{h}}^T,\hfill \end{array}$$

(16)

where the superscripts ${•}^T$ , ${•}^S$ , ${•}^E$ , and ${•}^D$ indicate constant stress, strain, electric field and dielectric displacement while ${ϵ}_0$ is the vacuum permittivity. Again, ${•}^{\mathrm{T}}$ means the transposition. The impermittivity ${\beta }^S$ is defined as the inverse of the permittivity tensor ${ϵ}^S$ and exists only if $det({ϵ}^S)\ne 0$ holds. For completeness, $\mathfrak{e}$ is also a third order piezoelectric tensor. The evolution of the internal variables will be presented now for the electrical and afterwards for the mechanical behavior.

Electrical behavior

As already mentioned, one aim of this constitutive model is to account for pinched hystereses. To do that, an approach from Linnemann (2008) for ferrimagnetic materials, where also comparable hystereses occur, is adapted. In the following figures, the load direction coincides with the 3-direction of the spatial coordinate system. Looking at Figure 2, one can see that the $D-E$ hysteresis is pinched so strongly, that two local hystereses occur. They look like normal but smaller dielectric hystereses. As a reference, the unpinched hysteresis is depicted by the black dotted line. The global behavior will now be expressed by the local hystereses via a coordinate transformation. This approach works for the electrical $D-E$ , $P^i-E$ , and $P-E$ hystereses as well as for the mechanical $T-S^{\mathrm{im}}$ and $T-S$ hystereses. Besides pinching, an internal order transition is taken into account for the electrical hystereses: The hystereses show a primary path, depicted by the dashed line as well as a secondary path, depicted by the solid line.

Figure 2.

Pinched $D_3-E_3$ hysteresis.

Looking at Figure 3, the upper right part of the $P^i-E$ hysteresis is depicted by the global and the newly introduced local coordinate system. To differentiate between global and local quantities, the latter are indexed with the additional uppercase subscript ${•}_L$ . The tilde $\widetilde{•}$ indicates path dependent quantities. Compared to the model in Schwaab et al. (2012), the newly introduced material parameters are the local electric coercive field $E_L^c$ , which influences the width of the local hystereses according to the physically reasonable range $\frac{E_L^c}{E^c}\in [0,1]$ as well as the order transition field $E^t$ and the order transition polarization $P^t$ . $E^c$ is the (global) electric coercive field and the nonnegative constant $c^e$ describes the slope of the $P^i-E$ hysteresis during poling. Based on these quantities, all other quantities for the electrical behavior of the model described in local coordinates will be derived from the hysteresis geometry.

Figure 3.

Upper right part of the pinched $P_3^i-E_3$ hysteresis.

The origin of the local coordinate system ${\overrightarrow{P}}_L^i-{\overrightarrow{E}}_L$ can either be in the first or third quadrant or as a degenerated case in the origin of the global coordinate system. This is expressed by the dielectric pinching function

$$\begin{array}{c}{\overrightarrow{g}}_d({\overrightarrow{P}}^i,\mathrm{grad}\left(\phi \right))=\hfill \\ \{\begin{array}{c}\frac{{\overrightarrow{P}}^i}{\Vert {\overrightarrow{P}}^i\Vert }\mathrm{if}\Vert {\overrightarrow{P}}^i\Vert >0\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{grad}\left(\phi \right)}{\Vert \mathrm{grad}\left(\phi \right)\Vert }\mathrm{else}\mathrm{if}\Vert {\overrightarrow{P}}^i\Vert =0\wedge \Vert \mathrm{grad}\left(\phi \right)\Vert >0}\\ \multicolumn{1}{c}{\overrightarrow{0}\mathrm{else},}\end{array}\hfill \end{array}$$

(17)

which implies three cases. In the first case, when $\Vert {\overrightarrow{P}}^i\Vert >0$ holds, the pinching function is equal to the direction of the irreversible polarization. In the second case, the negative direction of the gradient of the electric potential is used to avoid division by zero if $\Vert {\overrightarrow{P}}^i\Vert =0$ holds. If also $\Vert \mathrm{grad}\left(\phi \right)\Vert =0$ holds, the pinching function is set to zero. In the first two cases, $\Vert {\overrightarrow{g}}_d\Vert =1$ .

To distinguish between the primary and the secondary load path, the additional internal variable

$$p(t)=\underset{t}{max}\Vert {\overrightarrow{P}}^i(t)\Vert$$

(18)

with $p\in [0,P^{\mathrm{sat}}]$ is introduced. This variable remembers the maximum magnitude of irreversible polarization in the loading history. The primary load path passes the order transition point which is defined by the transition polarization $P^t$ and the electric transition field $E^t$ . This point is marked in Figures 3 and 4 by the circles. The local coercive field

$${\tilde{E}}_L^c=\{\begin{array}{c}{E}_L^c-(E^c-(E^t-c^e·P^t))\mathrm{if}p<P^{\mathrm{sat}}\hfill \\ \multicolumn{1}{c}{E_L^c\mathrm{else}\mathrm{if}p=P^{\mathrm{sat}},}\end{array}$$

(19)

Figure 4.

Upper right part of the pinched $D_3-E_3$ hysteresis.

depends on the internal variable $p$ which ensures that once the irreversible saturation polarization has been reached, the hysteresis path switches from the dashed to the solid line in Figure 3. To determine the position of the local coordinate system, the irreversible polarization offset

$$P_0^i=\frac{1}{2}{P}^{\mathrm{sat}},$$

(20)

and the electric field offset

$$E_0=E^c-E_L^c+c^e·P_0^i,$$

(21)

are needed. The transformation between global and local coordinates follows then via the offset and the dielectric pinching function for the irreversible polarization as

$${\overrightarrow{P}}_L^i={\overrightarrow{P}}^i-P_0^i·{\overrightarrow{g}}_d.$$

(22)

In Figure 4, the upper right part of the corresponding $D-E$ hystereses is depicted. For the transformation of the dielectric displacement between local and global coordinates, the offset

$$D_0=({ϵ}_0+{\kappa }^T)·E_0+P_0^i$$

(23)

is needed. The coordinate transformation for the dielectric displacement follows analog to equation (22) as

$${\overrightarrow{D}}_L=\overrightarrow{D}-D_0·{\overrightarrow{g}}_d.$$

(24)

With the local critical dielectric displacement

$${\tilde{D}}_L^c=({ϵ}_0+{\kappa }^T)·{\tilde{E}}_L^c,$$

(25)

the local dielectric switching function is defined to

$$f_L^d=\Vert {\overrightarrow{D}}_L-c^d·{\overrightarrow{P}}_L^i\Vert -{\tilde{D}}_L^c,$$

(26)

where the slopes $c^d=({ϵ}_0+{\kappa }^T)c^e+1$ and $c^{\mathrm{de}}=({ϵ}_0+{\kappa }^T)+\frac{1}{c^e}$ in Figure 4 are the kinematic hardening parameters. With the local saturation polarization

$$P_L^{\mathrm{sat}}=P^{\mathrm{sat}}-P_0^i·\Vert {\overrightarrow{g}}_d\Vert$$

(27)

which has to match the global saturation polarization $P^{\mathrm{sat}}$ for $\Vert {\overrightarrow{g}}_d\Vert =0$ , the local dielectric saturation function

$$h_L^d=\Vert {\overrightarrow{P}}_L^i\Vert -P_L^{\mathrm{sat}}$$

(28)

is defined.

The evolution equation of the internal variable ${\overrightarrow{P}}_L^i$ is defined as

$${\stackrel{·}{\overrightarrow{P}}}_L^i=\{\begin{array}{c}0\mathrm{if}{f}_L^d<0\wedge h_L^d<0\hfill \\ \multicolumn{1}{c}{{\lambda }_{f_L^d}\left(\frac{\partial f_L^d}{\partial {\overrightarrow{D}}_L}\right)\mathrm{else}\mathrm{if}{f}_L^d=0\wedge h_L^d<0}\\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{{\lambda }_{h_L^d}(-\frac{\partial h_L^d}{\partial {\overrightarrow{P}}_L^i})\mathrm{else}\mathrm{if}{h}_L^d=0}\end{array}$$

(29)

and includes three cases. In the first case, no evolution is possible. In the second case, the evolution of ${\overrightarrow{P}}_L^i$ takes place. The third case ensures, that $\Vert {\overrightarrow{P}}_L^i\Vert$ doesn’t exceed the saturation value $P_L^{\mathrm{sat}}$ . Equation (29) is a piece wisely defined, rate independent homogeneous ordinary differential equation (ODE) of order one. The Lagrange multiplicators ${\lambda }_{f_L^d}$ and ${\lambda }_{h_L^d}$ can be calculated by the consistency conditions ${\stackrel{·}{f}}_L^d({\lambda }_{f_L^d})=0$ and ${\stackrel{·}{h}}_L^d({\lambda }_{h_L^d})=0$ and have to fulfill the Kuhn-Tucker complementary conditions

$${\lambda }_{f_L^d}\ge 0,f_L^d\le 0,{\lambda }_{f_L^d}·f_L^d=0,$$

(30)

$${\lambda }_{h_L^d}\ge 0,h_L^d\le 0,{\lambda }_{h_L^d}·h_L^d=0.$$

(31)

The first condition is only valid if $h_L^d\le 0$ holds. Vice versa, if $h_L^d=0$ holds, $f_L^d>0$ is accepted.

In the evaluation of experiments, the order transition is associated with the inflection point of the primary load path in the $P-E$ and $S-E$ hystereses. If no order transition is observed in experiments, the transition polarization $P^t$ is set to the saturation polarization $P^{\mathrm{sat}}$ in the constitutive model. The order transition is modeled with the help of the smooth function

$$\beta (p)=\delta (\tanh \left(\gamma (\frac{p}{P^t}-1)\right)+1),$$

(32)

ranging between zero (primary state) and one (secondary state) in dependence on the internal variable $p$ . In Figure 5, this function is depicted for two cases wherein $\delta$ is defined in dependence of the transition polarization $P^t$ as

$$\delta =\{\begin{array}{c}\frac{1}{2}\mathrm{if}{P}^t<P^{\mathrm{sat}}\hfill \\ \multicolumn{1}{c}{1\mathrm{else}\mathrm{if}{P}^t=P^{\mathrm{sat}}.}\end{array}$$

(33)

Figure 5.

Evolution of $\beta (p)$ if $P^t<P^{\mathrm{sat}}$ (blue) and if $P^t=P^{\mathrm{sat}}$ (red).

The parameter $\gamma$ in equation (32) influences the sharpness of the order transition and has to be chosen sufficient large that the order transition is completed before $p=P^{\mathrm{sat}}$ , that is, the transition from the primary to the secondary load path, is reached. For materials with $P^t=P^{\mathrm{sat}}$ , both transitions occur simultaneously.

The evolution of irreversible polarization ${\overrightarrow{P}}^i$ causes the three dimensional and non-deviatoric irreversible polarization strain tensor

$$\begin{array}{c}{\mathit{S}}^{\mathrm{ip}}=\beta (p)·((\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}·(S_{\parallel }^{\mathrm{sat}}-S_{\parallel }^{\min })+S_{\parallel }^{\min })·({\overrightarrow{e}}^{P^i}\otimes {\overrightarrow{e}}^{P^i})\hfill \\ +(\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}·(S_{\perp }^{\mathrm{sat}}-S_{\perp }^{\min })+S_{\perp }^{\min })·(\mathit{I}-{\overrightarrow{e}}^{P^i}\otimes {\overrightarrow{e}}^{P^i}))\hfill \end{array}$$

(34)

which is strongly influenced by the order transition. This function depends linearly on $\beta (p)$ , on the irreversible polarization ${\overrightarrow{P}}^i$ and on the, compared to the formulations in Kamlah (2001); Schwaab et al. (2012), newly introduced material parameters. These are the saturation and minimum strain $S_{\parallel }^{\mathrm{sat}}\ge 0$ and $S_{\parallel }^{\min }\ge 0$ parallel, as well as $S_{\perp }^{\mathrm{sat}}\le 0$ and $S_{\perp }^{\min }\le 0$ perpendicular to the load direction. Figure 6 shows the $S-E$ hystereses for parallel (blue) and perpendicular (red) strain components. There, the total strains are the sums of corresponding elastic and irreversible polarization strains. As usual, the piezoelectric parameters $d_{33}$ and $d_{31}$ describe for a fully poled material the linear strain in, respectively perpendicular, to the poling direction induced by a parallel electric field. Before the order transition, $p<<P^{\mathrm{sat}}$ and thus $\beta (p)\approx 0$ , no irreversible polarization strain occurs at all. The hysteresis for the perpendicular strain component hast the same shape like the parallel strain component – but different sign and magnitude. The primary load paths, indicated by the dashed lines in Figure 6, rise respectively decline steeply when the electric field reaches $E^t$ . Due to the pinched behavior of the dielectric hysteresis, the $S-E$ hystereses are also pinched. In the secondary load path, the strain component parallel to the poling direction can’t fall below the strain level caused by the internal order transition, denoted here as $S_{\parallel }^{\min }$ , and the perpendicular strain component can’t be higher than the strain level denoted as $S_{\perp }^{\min }$ .

Figure 6.

$S_3-E_3$ (blue) and $S_1-E_3$ (red) hystereses with primary (dashed) and secondary (solid) load paths.

Materials without internal order transition but non-deviatoric polarization strain behavior can also be modeled by equation (34) by setting $\beta (p)=1$ and the material parameters to $S_{\parallel }^{\min }=S_{\perp }^{\min }=0$ . The relation then simplifies to

$${\mathit{S}}^{\mathrm{ip}}=\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}(S_{\parallel }^{\mathrm{sat}}{\overrightarrow{e}}^{P^i}\otimes {\overrightarrow{e}}^{P^i}+S_{\perp }^{\mathrm{sat}}(\mathit{I}-{\overrightarrow{e}}^{P^i}\otimes {\overrightarrow{e}}^{P^i})).$$

(35)

Setting furthermore $S^{\mathrm{sat}}=S_{\parallel }^{\mathrm{sat}}=-2S_{\perp }^{\mathrm{sat}}$ , equation (34), becomes the deviatoric irreversible polarization strain tensor

$${\mathit{S}}^{\mathrm{ip}}=\frac{3}{2}{S}^{\mathrm{sat}}\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}({\overrightarrow{e}}^{P^i}\otimes {\overrightarrow{e}}^{P^i}-\frac{1}{3}\mathit{I})$$

(36)

according to Kamlah (2001); Schwaab et al. (2012).

For a discussion of the formulation of the irreversible polarization strain tensor consider the following example: Assuming a fully poled material in for example, spatial 3-direction, ${\overrightarrow{P}}^i=[00P^{\mathrm{sat}}{]}^{\mathrm{T}}$ with $\frac{\Vert {\overrightarrow{P}}^i\Vert }{P^{\mathrm{sat}}}=1$ and $\beta (p)=1$ , then from equation (34) the maximum irreversible stain tensor

$${\mathit{S}}_{\max }^i=\left[\begin{array}{ccc}{S}_{\perp }^{\mathrm{sat}}\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& S_{\perp }^{\mathrm{sat}}\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & S_{\parallel }^{\mathrm{sat}}\hfill \end{array}\right]$$

(37)

follows with the norm $\Vert {\mathit{S}}_{\max }^i\Vert =\sqrt{{(S_{\parallel }^{\mathrm{sat}})}^2+2(S_{\perp }^{\mathrm{sat}}{)}^2}$ . The deviator of the maximum irreversible stain tensor is then

$$\begin{array}{c}{\mathit{S}}_{\max }^{i,\mathrm{Dev}}=\hfill \\ \left[\begin{array}{ccc}\frac{1}{3}(-S_{\parallel }^{\mathrm{sat}}+S_{\perp }^{\mathrm{sat}})\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{1}{3}(-S_{\parallel }^{\mathrm{sat}}+S_{\perp }^{\mathrm{sat}})\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \frac{2}{3}(S_{\parallel }^{\mathrm{sat}}-S_{\perp }^{\mathrm{sat}})\hfill \end{array}\right],\hfill \end{array}$$

(38)

with the norm

$$\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert =\sqrt{\frac{2}{3}}(S_{\parallel }^{\mathrm{sat}}-S_{\perp }^{\mathrm{sat}}).$$

(39)

Setting again $S^{\mathrm{sat}}=S_{\parallel }^{\mathrm{sat}}=-2S_{\perp }^{\mathrm{sat}}$ in equations (38) and (39), then $\Vert {\mathit{S}}_{\max }^i\Vert =\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert =\sqrt{\frac{3}{2}}{S}^{\mathrm{sat}}$ is obtained as it is known from Kamlah (2001); Schwaab et al. (2012). The norm of a tensor can be express via its first invariant and hence $\Vert {\mathit{S}}_{\max }^i\Vert$ as well as $\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert$ are independent of the polarization direction. The quantity $\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert$ is needed in the now following description of the mechanical behavior of the constitutive model.

Mechanical behavior

For the mechanical part of the constitutive model also the pinching effect is taken into account for the $T-S^{\mathrm{im}}$ and the $T-S$ hystereses analog to the pinched dielectric hystereses. In mechanically loaded experiments, only the stress-strain relation under compression has been explored due to the brittleness of the considered bulk ceramic materials. Hence, we focus on compressive states. However, the model represents a symmetric stress-strain hysteresis. The pinched ferroelastic hysteresis depicted in Figure 7 shows the mechanical strain and stress parallel to the load direction. This hysteresis is pinched so strongly, that local hystereses occur. As a reference, the unpinched hysteresis is also depicted, marked by the black dotted line. As before, the load direction coincides in the following figures with the 3-direction of the spatial coordinate system.

Figure 7.

Pinched $T_3-S_3$ hysteresis.

Analog to the dielectric hysteresis, the global hysteresis can be expressed with local hystereses and coordinate transformations using the mechanical pinching function. This function is defined as

$$\begin{array}{c}{\mathit{g}}_m({\mathit{S}}^{\mathrm{im}},{\mathit{T}}^{\mathrm{Dev}})=\hfill \\ \{\begin{array}{c}\sqrt{\frac{3}{2}}\frac{{\mathit{S}}^{\mathrm{im}}}{\Vert {\mathit{S}}^{\mathrm{im}}\Vert }\mathrm{if}\Vert {\mathit{S}}^{\mathrm{im}}\Vert >0\hfill \\ \multicolumn{1}{c}{\sqrt{\frac{3}{2}}\frac{{\mathit{T}}^{\mathrm{Dev}}}{\Vert {\mathit{T}}^{\mathrm{Dev}}\Vert }\mathrm{else}\mathrm{if}\Vert {\mathit{S}}^{\mathrm{im}}\Vert =0\wedge \Vert {\mathit{T}}^{\mathrm{Dev}}\Vert >0}\\ \multicolumn{1}{c}{0\mathrm{else}}\end{array}\hfill \end{array}$$

(40)

and depends on the internal variable of irreversible mechanical strain ${\mathit{S}}^{\mathrm{im}}$ and on the deviatoric stress tensor ${\mathit{T}}^{\mathrm{Dev}}=\mathit{T}-\frac{1}{3}\mathrm{tr}(\mathit{T})·\mathit{I}$ . Three cases are implied: In the first case, when $\Vert {\mathit{S}}^{\mathrm{im}}\Vert >0$ , ${\mathit{g}}_m$ is equal to the normalized irreversible mechanical strain tensor. In the second case, the normalized stress deviator is used to avoid division by zero if $\Vert {\mathit{S}}^{\mathrm{im}}\Vert =0$ . In the third case, when also $\Vert {\mathit{T}}^{\mathrm{Dev}}\Vert =0$ , the pinching function is set to $0$ . The factor $\sqrt{\frac{3}{2}}$ in the first and second cases forces the entries of ${\mathit{g}}_m$ to be $\pm 1$ , $\mp \frac{1}{2}$ or $0$ . Relying on deviatoric quantities in modeling the mechanical behavior is motivated by pure switching processes being of shape changing nature only.

The origin of the local coordinate system $T_L-S_L^{\mathrm{im}}$ can lie in the first (tensile) and third (compressive) quadrant or in the origin of the global coordinate system $T-S^{\mathrm{im}}$ , depending on the equation (40). Let’s consider now the compressive part of the global $T_3-S_3^{\mathrm{im}}$ hysteresis, depicted in Figure 8, which shows the irreversible mechanical strain and stress parallel to the load direction. There, the local $T_L-S_L^{\mathrm{im}}$ coordinate system is placed in the middle of the local hysteresis. For completeness, the corresponding $T_3-S_3$ hysteresis is depicted in Figure 9, containing additionally the reversible strain. The height of the local $(T_L-S_L^{\mathrm{im}})$ hysteresis in Figure 8 and the local $(T_L-S_L)$ hysteresis in Figure 9, respectively, are determined from the newly introduced local mechanical coercive stress $T_L^c$ . All other quantities for the mechanical behavior of the model described in local coordinates will be derived from the hysteresis geometry.

Figure 8.

Compressive part of the pinched $T_3-S_3^{\mathrm{im}}$ hysteresis.

Figure 9.

Compressive part of the pinched $T_3-S_3$ hysteresis.

The irreversible strain offset

$$S_0^i=\frac{1}{2}(\sqrt{\frac{2}{3}}\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert -\sqrt{\frac{2}{3}}\Vert {\mathit{S}}^{\mathrm{ip},\mathrm{Dev}}\Vert )$$

(41)

is a function of the maximum possible irreversible deviatoric strain state according to equation (39), lowered by the amount of irreversible strain induced by the polarization.

As reported among others in Kamlah (2001); Schwaab et al. (2012), the switching of domains under mechanical load will either be stabilized or destabilized by a superposed electric field, depending on the strength of $\overrightarrow{E}$ and the relative orientation between $\overrightarrow{E}$ and ${\overrightarrow{P}}^i$ . This behavior is taken into account by introducing the electric field dependent mechanical coercive stress

$${\hat{T}}^c(\overrightarrow{E},{\overrightarrow{P}}^i)=(T^c-T_{\min }^c)〈1+\frac{h}{T^c}\frac{\overrightarrow{E}·{\overrightarrow{e}}^{P^i}}{E^c}〉+T_{\min }^c.$$

(42)

This is the stress level where domains start non- ${180}^{°}$ switching. A superposed electric field decreases ( $\overrightarrow{E}$ and ${\overrightarrow{P}}^i$ in opposite directions) or increases ( $\overrightarrow{E}$ and ${\overrightarrow{P}}^i$ in same directions) this stress level. This function depends on the mechanical coercive stress $T^c$ , on the minimum coercive stress $T_{\min }^c>0$ to ensure always ${\hat{T}}^c\ge T_{\min }^c>0$ to prevent numerical issues, and on a coupling parameter $h\ge 0$ which describes the strength of the influence of the electric field on mechanically induced switching. The denominators of the fractions in the bracket are introduced to get the proper units. $〈x〉$ is the bracket function

$$〈x〉=\{\begin{array}{c}x\mathrm{if}x\ge 0\hfill \\ \multicolumn{1}{c}{0\mathrm{else}\mathrm{if}x<0.}\end{array}$$

(43)

Without superposed electric field, ${\hat{T}}^c=T^c$ . If the directions of the electric field and the polarization are the same, then ${\hat{T}}^c>T^c$ . If the directions are opposite to each other, than $T^c>{\hat{T}}^c\ge T_{\min }^c>0$ . It is worth mentioning that equation (42) is a generalization of the corresponding function in Schwaab et al. (2012) and by setting $T_{\min }^c=0$ , these functions are identical. To incorporate the electric field dependent global mechanical coercive stress from equation (42) into the local hystereses, the electric field dependent local mechanical coercive stress

$${\hat{T}}_L^c={\hat{T}}^c-T^c+T_L^c$$

(44)

is determined. With the stress offset

$$T_0=T^c-T_L^c+\frac{3}{2}{c}^m·S_0^i,$$

(45)

depending on the irreversible strain offset $S_0^i$ defined in equation (41), the position of the local coordinate system $(T_L-S_L^{\mathrm{im}})$ is determined. The slope of the linear kinematic hardening is described by $c^m>0$ in the $(T_L-S_L^{\mathrm{im}})$ hysteresis and

$$c^t={(\frac{1}{Y}+\frac{1}{\frac{3}{2}{c}^m})}^{-1}$$

(46)

in the $(T_L-S_L)$ hysteresis, respectively. With the help of ${\mathit{g}}_m$ from equation (40) and the strain and stress offsets equations (41) and (45), the local stress tensor

$${\mathit{T}}_L=\mathit{T}-T_0·\frac{2}{3}·{\mathit{g}}_m$$

(47)

as well as the local irreversible mechanical strain tensor

$${\mathit{S}}_L^{\mathrm{im}}={\mathit{S}}^{\mathrm{im}}-S_0^i·{\mathit{g}}_m$$

(48)

are determined.

Now, the local mechanical switching function

$$f_L^m=\Vert {({\mathit{T}}_L-c^m{\mathit{S}}_L^{\mathrm{im}})}^{\mathrm{Dev}}\Vert -\sqrt{\frac{2}{3}}{\hat{T}}_L^c,$$

(49)

is assumed analogously to the classical yield function for associate von-Mises plasticity with kinematic hardening. The last missing quantity for the local mechanical saturation function is the local irreversible saturation strain

$$S_L^{\mathrm{sat}}=\{\begin{array}{c}\sqrt{\frac{2}{3}}\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert -\sqrt{\frac{2}{3}}\Vert {\mathit{S}}^{\mathrm{ip},\mathrm{Dev}}\Vert -S_0^i\mathrm{if}\Vert {\mathit{g}}_m\Vert >0\hfill \\ \multicolumn{1}{c}{\sqrt{\frac{2}{3}}\Vert {\mathit{S}}_{\max }^{i,\mathrm{Dev}}\Vert -\sqrt{\frac{2}{3}}\Vert {\mathit{S}}^{\mathrm{ip},\mathrm{Dev}}\Vert \mathrm{else}\mathrm{if}\Vert {\mathit{g}}_m\Vert =0}\end{array}$$

(50)

which contains two cases in dependence on $\Vert {\mathit{g}}_m\Vert$ . The second case is only needed if $\Vert {\mathit{g}}_m\Vert =0$ to match the global saturation strain. The local mechanical saturation function is then defined as

$$h_L^m=\Vert {\mathit{S}}_L^{\mathrm{im}}\Vert -\sqrt{\frac{3}{2}}{S}_L^{\mathrm{sat}}.$$

(51)

The evolution equation for the internal variable ${\mathit{S}}_L^{\mathrm{im}}$ is now given by

$${\stackrel{·}{\mathit{S}}}_L^{\mathrm{im}}=\{\begin{array}{c}0\mathrm{if}{f}_L^m<0\wedge h_L^m<0\hfill \\ \multicolumn{1}{c}{{\lambda }_{f_L^m}\left(\frac{\partial f_L^m}{\partial {\mathit{T}}_L}\right)\mathrm{else}\mathrm{if}{f}_L^m=0\wedge h_L^m<0}\\ \multicolumn{1}{c}{{\lambda }_{h_L^m}(-\frac{\partial h_L^m}{\partial {\mathit{S}}_L^{\mathrm{im}}})\mathrm{else}\mathrm{if}{h}_L^m=0}\end{array}$$

(52)

and includes three cases. In the first case, the switching as well as the saturation function are lower than zero and no evolution is possible. In the second case, when the switching function is equal to zero and the saturation function is less than zero the evolution of ${\mathit{S}}_L^{\mathrm{im}}$ takes place. The third case ensures that $\Vert {\mathit{S}}_L^{\mathrm{im}}\Vert$ doesn’t exceed $S_L^{\mathrm{sat}}$ even if the switching function may be higher then zero. Equation (52) is a rate independent homogeneous ODE of order one. The Lagrange multiplicators ${\lambda }_{f_L^m}$ and ${\lambda }_{h_L^m}$ can be calculated by the consistency conditions ${\stackrel{·}{f}}_L^m({\lambda }_{f_L^m})=0$ and ${\stackrel{·}{h}}_L^m({\lambda }_{h_L^m})=0$ and have to fulfill the Kuhn-Tucker complementary conditions

$${\lambda }_{f_L^m}\ge 0,f_L^m\le 0,{\lambda }_{f_L^m}·f_L^m=0,$$

(53)

$${\lambda }_{h_L^m}\ge 0,h_L^m\le 0,{\lambda }_{h_L^m}·h_L^m=0.$$

(54)

The first condition is only valid if $h_L^m\le 0$ holds. Vice versa, if $h_L^m=0$ holds, $f_L^m>0$ is accepted. Even if the evolution of the internal variables is calculated by help of the local hystereses, only global variables will enter the balance equations.

Time discrete integration algorithm

As proposed in Schwaab et al. (2012), the two nonlinear ODE’s for ${\overrightarrow{P}}_L^i$ and ${\mathit{S}}_L^{\mathrm{im}}$ , equations (29) and (52), are integrated computationally via a predictor-corrector procedure called return mapping algorithm according to Simo and Hughes (2000) in the framework of the iterative finite-element-method solving procedure. This approach was introduced to model plastic behavior of metals using a so called yield function with a corresponding yield surface. If the predictor lies outside of the yield surface, it will be projected radially back to it via the corrector. This approach will now be adapted for the switching and saturation surfaces. Before the calculation starts, the internal variables ${}^{0}p=0$ , ${{}^0\overrightarrow{P}}^i=\overrightarrow{0}$ , and ${{}^0\mathit{S}}^{\mathrm{im}}=0$ are initialized to zero and ${{}^0\hat{T}}^c=T^c$ on each integration point. Superscripts ${}^n•$ and ${}^{n+1}•$ denote the known quantities at the end of the former and the current time step, respectively. The current primary variables ${}^{n+1}\overrightarrow{u}$ , ${}^{n+1}\phi$ and ${}^{n+1}\overrightarrow{D}$ result from the finite-element iteration. At the beginning of every iteration, ${\lambda }_{f_L^d}$ , ${\lambda }_{h_L^d}$ , ${\lambda }_{f_L^m}$ , and ${\lambda }_{h_L^m}$ are set to zero and only if the according switching and saturation functions are not fulfilled, these quantities differ from zero. The local irreversible polarization and the local irreversible mechanical strain at the end of the current time step,

$${{}^{n+1}\overrightarrow{P}}_L^i{=}^n{\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,h}$$

(55)

$${{}^{n+1}\mathit{S}}_L^{\mathrm{im}}{=}^n{\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},h}$$

(56)

are computed from the quantities of the former time step, and the corresponding correctors $\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}$ , $\mathrm{\Delta }{\overrightarrow{P}}_L^{i,h}$ , $\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}$ , and $\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},h}$ .

The hierarchically ordered evaluation of the switching and saturation criteria starts with the dielectric switching function, equation (26). With the help of the discrete form

$$\begin{array}{c}{{}^{n+1}\overrightarrow{g}}^d({{}^n\overrightarrow{P}}^i,{}^{n+1}\phi )=\hfill \\ \{\begin{array}{c}\frac{{{}^n\overrightarrow{P}}^i}{\Vert {{}^n\overrightarrow{P}}^i\Vert }\mathrm{if}\Vert {{}^n\overrightarrow{P}}^i\Vert >0\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{grad}\left({}^{n+1}\phi \right)}{\Vert \mathrm{grad}\left({}^{n+1}\phi \right)\Vert }\mathrm{else}\mathrm{if}\Vert {{}^n\overrightarrow{P}}^i\Vert =0\wedge \Vert \mathrm{grad}\left({}^{n+1}\phi \right)\Vert >0}\\ \multicolumn{1}{c}{\overrightarrow{0}\mathrm{else}}\end{array}\hfill \end{array}$$

(57)

of the dielectric pinching function in equation (17), the local quantities ${{}^n\overrightarrow{P}}_L^i$ and ${{}^{n+1}\overrightarrow{D}}_L$ are determined from equations (22) and (24).

If the dielectric switching function exceeds zero, $f_L^d\left({{}^{n+1}\overrightarrow{D}}_L{,}^n{\overrightarrow{P}}_L^i\right)>0$ , the first corrector is computed to

$$\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}={\lambda }_{f_L^d}\frac{\partial f_L^d}{\partial {\overrightarrow{D}}_L}$$

(58)

with the magnitude

$${\lambda }_{f_L^d}=\frac{1}{c^d}(\Vert {{}^{n+1}\overrightarrow{D}}_L-c^d{·}^n{\overrightarrow{P}}_L^i\Vert -{\tilde{D}}_L^c)$$

(59)

and the direction

$$\frac{\partial f_L^d}{\partial {\overrightarrow{D}}_L}=\frac{{{}^{n+1}\overrightarrow{D}}_L-c^d{·}^n{\overrightarrow{P}}_L^i}{\Vert {{}^{n+1}\overrightarrow{D}}_L-c^d{·}^n{\overrightarrow{P}}_L^i\Vert }$$

(60)

to satisfy $f_L^d({{}^{n+1}\overrightarrow{D}}_L{,}^n{\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f})=0$ . As a major simplification, $\frac{\partial f_L^d}{\partial {\overrightarrow{D}}_L}$ is considered to be constant during the time step.

If the dielectric saturation function, equation (28), exceeds zero, $h^d({{}^n\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f})>0$ , the second corrector is computed with the help of the already known corrector $\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}$ to

$$\mathrm{\Delta }{\overrightarrow{P}}_L^{i,h}={\lambda }_{h_L^d}(-\frac{\partial h_L^d}{\partial {\overrightarrow{P}}_L^i}),$$

(61)

with the magnitude

$${\lambda }_{h_L^d}=\Vert {{}^n\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}\Vert -P_L^{\mathrm{sat}}$$

(62)

and the direction

$$\frac{\partial h_L^d}{\partial {\overrightarrow{P}}_L^i}=\frac{{{}^n\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}}{\Vert {{}^n\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}\Vert }$$

(63)

to satisfy $h^d({{}^n\overrightarrow{P}}_L^i+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,f}+\mathrm{\Delta }{\overrightarrow{P}}_L^{i,h})=0$ . The saturation criterion is of higher priority than the switching criterion. Consequently, if the saturation criterion is fulfilled, the switching criterion is allowed to exceed $f_L^d=0$ . Knowing the two dielectric correctors, ${{}^{n+1}\overrightarrow{P}}_L^i$ is evaluated according to equation (55) and the current global irreversible polarization is obtained as

$${{}^{n+1}\overrightarrow{P}}^i=P_0^i{·}^{n+1}{\overrightarrow{g}}_d{+}^{n+1}{\overrightarrow{P}}_L^i.$$

(64)

Now, the direction vector of the irreversible polarization can be calculated:

$$\begin{array}{c}{{}^{n+1}\overrightarrow{e}}^{P^i}\left({{}^{n+1}\overrightarrow{P}}^i{,}^{n+1}\phi \right)=\hfill \\ \{\begin{array}{c}\frac{{{}^{n+1}\overrightarrow{P}}^i}{\Vert {{}^{n+1}\overrightarrow{P}}^i\Vert }\mathrm{if}\Vert {{}^{n+1}\overrightarrow{P}}^i\Vert >0\hfill \\ \multicolumn{1}{c}{-\frac{\mathrm{grad}\left({}^{n+1}\phi \right)}{\Vert \mathrm{grad}\left({}^{n+1}\phi \right)\Vert }\mathrm{else}\mathrm{if}\Vert {{}^{n+1}\overrightarrow{P}}^i\Vert =0\wedge \Vert \mathrm{grad}\left({}^{n+1}\phi \right)\Vert >0.}\\ \multicolumn{1}{c}{\overrightarrow{0}\mathrm{else}}\end{array}\hfill \end{array}$$

(65)

Note the different time indices of ${\overrightarrow{P}}^i$ in equation (65) in contrast to equation (57). The internal variable

$${}^{n+1}p=max({}^{n}p,\Vert {{}^{n+1}\overrightarrow{P}}^i\Vert )$$

(66)

then gets updated and with

$${}^{n+1}\beta \left({}^{n+1}p\right)=\delta (\tanh \left(\gamma (\frac{{}^{n+1}p}{P^t}-1)\right)+1)$$

(67)

the current irreversible polarization strain ${{}^{n+1}\mathit{S}}^{ip}{(}^{n+1}\beta {,}^{n+1}{\overrightarrow{P}}^i)$ is determined from equation (34).

The dielectric internal variables are now determined and the mechanical ones follows now in the same way to determine the irreversible mechanical strain ${{}^{n+1}\mathit{S}}^{\mathrm{im}}$ . For the return mapping procedure, the global $\mathrm{trial}$ stress

$$\begin{array}{c}\multicolumn{1}{c}{{}^{\mathrm{trial}}\mathit{T}{=}^{n+1}{\mathbb{C}}^D:{(}^{n+1}\mathit{S}{-}^{n+1}{\mathit{S}}^{\mathrm{ip}}{-}^n{\mathit{S}}^{\mathrm{im}})}\\ {-}^{n+1}\mathfrak{h}·\left({}^{n+1}\overrightarrow{D}{-}^{n+1}{\overrightarrow{P}}^i\right),\hfill \end{array}$$

(68)

is evaluated according to the constitutive equation (13) with the so far known quantities. Only the internal variable ${{}^n\mathit{S}}^{\mathrm{im}}$ comes from the former time step. With the help of the discrete form

$$\begin{array}{c}{{}^{n+1}\mathit{g}}_m({{}^n\mathit{S}}^{\mathrm{im}},{{}^n\mathit{T}}^{\mathrm{Dev}})=\hfill \\ \{\begin{array}{c}\sqrt{\frac{3}{2}}\frac{{{}^n\mathit{S}}^{\mathrm{im}}}{\Vert {{}^n\mathit{S}}^{\mathrm{im}}\Vert }\mathrm{if}\Vert {{}^n\mathit{S}}^{\mathrm{im}}\Vert >0\hfill \\ \multicolumn{1}{c}{\sqrt{\frac{3}{2}}\frac{{{}^n\mathit{T}}^{\mathrm{Dev}}}{\Vert {{}^n\mathit{T}}^{\mathrm{Dev}}\Vert }\mathrm{else}\mathrm{if}\Vert {{}^n\mathit{S}}^{\mathrm{im}}\Vert =0\wedge \Vert {{}^n\mathit{T}}^{\mathrm{Dev}}\Vert >0}\\ \multicolumn{1}{c}{0\mathrm{else}}\end{array}\hfill \end{array}$$

(69)

of the mechanical pinching function, depending on irreversible mechanical strain and deviatoric stress tensors from the former time step, the local quantities ${{}^n\mathit{S}}_L^{\mathrm{im}}$ and ${{}^{\mathrm{trial}}\mathit{T}}_L$ are determined by equations (47) and (48) and used in the mechanical switching function equation (49). If this function exceeds zero, $f_L^m\left({{}^{\mathrm{trial}}\mathit{T}}_L{,}^n{\mathit{S}}_L^{\mathrm{im}}{,}^n{\hat{T}}_L^c\right)>0$ , the third corrector is computed as

$$\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}={\lambda }_{f_L^m}\frac{\partial f_L^m}{\partial {\mathit{T}}_L}$$

(70)

with the magnitude

$${\lambda }_{f_L^m}=\left(\frac{\mathit{A}:\mathit{B}-\sqrt{{(\mathit{A}:\mathit{B})}^2-(\mathit{B}:\mathit{B})·(\mathit{A}:\mathit{A}-{\left(\sqrt{\frac{2}{3}}{}^n{\hat{T}}_L^c\right)}^2)}}{\mathit{B}:\mathit{B}}\right)$$

(71)

wherein the abbreviations

$$\mathit{A}={({{}^{\mathrm{trial}}\mathit{T}}_L-c^m{·}^n{\mathit{S}}_L^{\mathrm{im}})}^{\mathrm{Dev}}$$

(72)

and

$$\mathit{B}={(}^{n+1}{\mathbb{C}}^D:\mathit{N}+c^m·\mathit{N}{)}^{\mathrm{Dev}}$$

(73)

and the direction

$$\frac{\partial f_L^m}{\partial {\mathit{T}}_L}=\frac{{({{}^{\mathrm{trial}}\mathit{T}}_L-c^m{·}^n{\mathit{S}}_L^{\mathrm{im}})}^{\mathrm{Dev}}}{\Vert {({{}^{\mathrm{trial}}\mathit{T}}_L-c^m{·}^n{\mathit{S}}_L^{\mathrm{im}})}^{\mathrm{Dev}}\Vert }=\mathit{N}$$

(74)

are used to satisfy $f_L^m({{}^{\mathrm{trial}}\mathit{T}}_L{,}^n{\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}{,}^n{\hat{T}}^c)=0$ . As a simplification, $\frac{\partial f_L^m}{\partial {\mathit{T}}_L}$ is considered to be constant during the time step.

If the mechanical saturation function in equation (51) exceeds zero, $h_L^m({{}^n\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}{,}^{n+1}{\mathit{S}}^{\mathrm{ip}})>0$ , the fourth corrector is computed as

$$\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},h}={\lambda }_{h_L^m}(-\frac{\partial h_L^m}{\partial {\mathit{S}}_L^{\mathrm{im}}})$$

(75)

with the magnitude

$${\lambda }_{h_L^m}=\Vert {{}^n\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}\Vert -\sqrt{\frac{3}{2}}{S}_L^{\mathrm{sat}}$$

(76)

and the direction

$$\frac{\partial h_L^m}{\partial {\mathit{S}}_L^{\mathrm{im}}}=\frac{{{}^n\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}}{\Vert {{}^n\mathit{S}}_L^{\mathrm{im}}+\mathrm{\Delta }{\mathit{S}}_L^{\mathrm{im},f}\Vert }$$

(77)

to satisfy $h_L^m{(}^n{\mathit{S}}_L^{im}+\mathrm{\Delta }{\mathit{S}}_L^{im,f}+\mathrm{\Delta }{\mathit{S}}_L^{im,h}{,}^{n+1}{\mathit{S}}^{ip})=0$ . The mechanical saturation criterion is of higher priority than the mechanical switching criterion. This means, if the saturation criterion is fulfilled, the switching criterion is allowed to exceed $f_L^m=0$ . Knowing the two mechanical correctors, ${{}^{n+1}\mathit{S}}_L^{\mathrm{im}}$ is evaluated according to equation (56) and the current global irreversible mechanical strain is computed as

$${{}^{n+1}\mathit{S}}^{\mathrm{im}}=S_0^i{·}^{n+1}{\mathit{g}}_m{+}^{n+1}{\mathit{S}}_L^{\mathrm{im}}.$$

(78)

The material tensors ${}^{n+1}\mathfrak{h}$ , ${{}^{n+1}\beta }^S$ , and ${{}^{n+1}\mathbb{C}}^D$ are also functions of ${{}^{n+1}\overrightarrow{P}}^i$ and therefore they need to be evaluated in every iteration. With the current strain tensor

$${}^{n+1}\mathit{S}=\frac{1}{2}(\mathrm{grad}\left({}^{n+1}\overrightarrow{u}\right)+\mathrm{grad}{\left({}^{n+1}\overrightarrow{u}\right)}^{\mathrm{T}})$$

(79)

the constitutive equations

$$\begin{array}{c}\multicolumn{1}{c}{{}^{n+1}\overrightarrow{E}={-}^{n+1}{\mathfrak{h}}^{\mathrm{T}}:\left({}^{n+1}\mathit{S}{-}^{n+1}{\mathit{S}}^{\mathrm{ip}}{-}^{n+1}{\mathit{S}}^{\mathrm{im}}\right)}\\ {+}^{n+1}{\beta }^S·\left({}^{n+1}\overrightarrow{D}{-}^{n+1}{\overrightarrow{P}}^i\right)\hfill \end{array}$$

(80)

$$\begin{array}{c}\multicolumn{1}{c}{{}^{n+1}\mathit{T}={{}^{n+1}\mathbb{C}}^D:\left({}^{n+1}\mathit{S}{-}^{n+1}{\mathit{S}}^{\mathrm{ip}}{-}^{n+1}{\mathit{S}}^{\mathrm{im}}\right)}\\ {-}^{n+1}\mathfrak{h}·\left({}^{n+1}\overrightarrow{D}{-}^{n+1}{\overrightarrow{P}}^i\right)\hfill \end{array}$$

(81)

according to equations (12) and (13) are evaluated.

Now the current generalized mechanical coercive stress

$${{}^{n+1}\hat{T}}^c=(T^c-T_{\min }^c)〈1+\frac{h}{T^c}\frac{{}^{n+1}\overrightarrow{E}·{{}^{n+1}\overrightarrow{e}}^{P^i}}{E^c}〉+T_{\min }^c$$

(82)

is updated and saved for the next time step. The flowchart of the time discrete integration algorithm, which needs to be evaluated at all integration points in all elements for all (pseudo) times steps, is depicted in Figure 10.

Figure 10.

Flowchart of time discrete integration algorithm.

The here presented time discrete material model is not a fully implicit integration algorithm. Therefore, small time steps are needed to get sufficient and reliable simulation results.

Examples

For the following examples, the three dimensional constitutive model is implemented in the commercial finite-element program COMSOL Multiphysics^™ using the equation based modeling module for partial differential equations. Simulations of simple cubes consisting of two quadratic tetrahedral elements are conducted. Minimal necessary mechanical boundary conditions are applied.

The first example, Figure 11, demonstrates quantitatively the influence of the local coercive field strength $E_L^c$ and the local coercive stress $T_L^c$ on the pinching effect of the corresponding hystereses. Fictitious material parameters are used and a cyclic electric field or a cyclic mechanical stress in spatial 3-direction, respectively, is applied. Let’s focus on the left colum in Figure 11, containing the $P-E$ hystereses for several $\frac{E_L^c}{E^c}$ ratios. For all ratios, the primary load paths are identical. The differences in the secondary load path will be discussed now. For $\frac{E_L^c}{E^c}=1$ the hysteresis is, except for the primary path, identical to the model in Schwaab et al. (2012) and shows no pinching effect. For $\frac{E_L^c}{E^c}=\frac{3}{4}$ , a slightly pinched hysteresis is visible. By reducing the ratio to $\frac{E_L^c}{E^c}=\frac{1}{2}$ , two local hystereses occur, which touch each other in the origin. Further decrease of the ratio, here shown for $\frac{E_L^c}{E^c}=\frac{1}{4}$ , causes two smaller and separated local hystereses. Let’s focus now on the center column in Figure 11, containing the corresponding $S-E$ ferroelectric hystereses for the strain component parallel to the load direction. For all $\frac{E_L^c}{E^c}$ ratios, the evolution of the internal order transition is identical. After the order transition, the minimum strain is $S_{\parallel }^{\min }$ which is assumed to be $0.5·S_{\parallel }^{\mathrm{sat}}$ here. For $\frac{E_L^c}{E^c}=1$ , a hysteresis with sharp tips as well as lines with linear slopes on the upper region is observed. These two lines intersect each other at zero electric field. Decreasing the ratio to $\frac{E_L^c}{E^c}=\frac{3}{4}$ , plateaus on minimum strain level form on the left and right hand side of the hysteresis. The slope of the upper lines shows a drop before zero electric field is reached. For $\frac{E_L^c}{E^c}=\frac{1}{2}$ , the plateaus unite and two connected local hystereses are observed. Further decrease of the ratio leads to more and more separated and also thinner local hystereses, as can be seen for $\frac{E_L^c}{E^c}=\frac{1}{4}$ .

Figure 11.

Influence of $E_L^c$ on $P_3-E_3$ (left) and $S_3-E_3$ (middle) hystereses, influence of $T_L^c$ on $T_3-S_3$ (right) hystereses. Fictitious material parameters.

Let’s focus now on the right column in Figure 11, containing the ferroelastic $T-S$ hystereses for several $\frac{T_L^c}{T^c}$ ratios. Depicted are the stress and strain components parallel to the load direction. For $\frac{T_L^c}{T^c}=1$ , no pinching is observed and the hysteresis is identical to the models in Kamlah (2001); Schwaab et al. (2012). Decreasing the ratio to $\frac{T_L^c}{T^c}=\frac{3}{4}$ , the hysteresis is slightly pinched. For $\frac{T_L^c}{T^c}=\frac{1}{2}$ , two local hystereses occur, which touch each other in the origin. Hystereses with a ratio $\frac{T_L^c}{T^c}\le 0.5$ show a pseudoelastic behavior known from shape-memory alloys Morin et al. (2011).

The second example, Figure 12, shows the constitutive model (black) adjusted to the experimental data (blue/red) of NBT-6BT and NBT-6BT-4KNN, respectively. The related material parameters are listed in Table 1. From the conducted experiments, the parameters $\nu$ , $d_{15}$ , $T_{\min }^c$ , and $h$ can’t be derived and therefore they are assumed. The simulated $P-E$ hysteresis of NBT-6BT shows an internal order transition with primary and secondary load path and matches the experiment well. Also for NBT-6BT-4KNN the simulated $P-E$ is in good agreement with the experiment. Here, the pinching effect is visible, but no internal order transition. The $S-E$ hystereses (parallel as well as perpendicular to the load direction) from simulation and experiment match well for both materials. Only minor derivations at the outer tips of the hysteresis of NBT-6BT and on the primary load path of NBT-6BT-4KNN are visible. The simulated $T-S$ hystereses match also well, except for the unloading path for NBT-6BT. This is caused by the linearity of the kinematic hardening of the switching and saturation functions. Concluding, all simulated hystereses are very close to the experiment and the constitutive model works well for the lead-free ferroelectric and the relaxor materials investigated here.

Figure 12.

Simulated (black) and measured $P_3-E_3$ , $S_3-E_3$ (blue) and $S_1-E_3$ (red) as well as $T_3-S_3$ hystereses of NBT-6BT (left) and NBT-6BT-4KNN (right).

Table 1.

Material parameter.

	NBT-6BT	NBT-6BT-4KNN	Units
$E^c$	$2.20·{10}^6$	$1.60·{10}^6$	$\mathrm{V}/\mathrm{m}$
$E_L^c$	$2.20·{10}^6$	$1.36·{10}^6$	$\mathrm{V}/\mathrm{m}$
$E^t$	$2.44·{10}^6$	$3.57·{10}^6$	$\mathrm{V}/\mathrm{m}$
$P^{\mathrm{sat}}$	$4.20·{10}^{-1}$	$2.90·{10}^{-1}$	$\mathrm{C}/{\mathrm{m}}^2$
$P^t$	$3.16·{10}^{-1}$	$2.90·{10}^{-1}$	$\mathrm{C}/{\mathrm{m}}^2$
$c^e$	$2.75·{10}^6$	$6.80·{10}^6$	${\mathrm{V}}^2/\mathrm{N}$
$d_{33}$	$0.26·{10}^{-9}$	$0.42·{10}^{-9}$	$\mathrm{m}/\mathrm{V}$
$d_{31}$	$-0.10·{10}^{-9}$	$-0.12·{10}^{-9}$	$\mathrm{m}/\mathrm{V}$
$d_{15}$	$0.20·{10}^{-9}$	$0.20·{10}^{-9}$	$\mathrm{m}/\mathrm{V}$
${\kappa }^T$	$14.00·{10}^{-9}$	$26.00·{10}^{-9}$	$\mathrm{N}/{\mathrm{V}}^2$
$S_{\parallel }^{\mathrm{sat}}$	$27.50·{10}^{-4}$	$5.00·{10}^{-4}$	−
$S_{\perp }^{\mathrm{sat}}$	$-14.10·{10}^{-4}$	$-0.55·{10}^{-4}$	−
$S_{\parallel }^{\min }$	$17.25·{10}^{-4}$	$1.37·{10}^{-4}$	−
$S_{\perp }^{\min }$	$-9.85·{10}^{-4}$	$-0.10·{10}^{-4}$	−
$T^c$	$0.12·{10}^9$	$0.22·{10}^9$	$\mathrm{N}/{\mathrm{m}}^2$
$T_L^c$	$0.12·{10}^9$	$0.22·{10}^9$	$\mathrm{N}/{\mathrm{m}}^2$
$T_{\min }^c$	$6.00·{10}^6$	$11.00·{10}^6$	$\mathrm{N}/{\mathrm{m}}^2$
$c^m$	$65.00·{10}^9$	$500.00·{10}^9$	$\mathrm{N}/{\mathrm{m}}^2$
$Y$	$117.00·{10}^9$	$104.00·{10}^9$	$\mathrm{N}/{\mathrm{m}}^2$
$h$	$0.06·{10}^9$	$0.11·{10}^9$	$\mathrm{N}/{\mathrm{m}}^2$
$\nu$	$0.33$	$0.33$	−
$\gamma$	$15.00$	$8.00$	−

Summary and outlook

The here presented three dimensional and fully electromechanically coupled constitutive model is able to capture well the behavior of the lead-free ferroelectric and relaxor materials NBT-6BT and NBT-6BT-4KNN, respectively. It takes into account the pinching effect, the internal order transition and the non-deviatoric polarization induced strain tensor observed in experiments. In principle, this model may also work for lead based ferroelectrics, and it covers materials with unpinched hysteresis as a special case. Time integration of the history dependent internal variables is done with a predictor-corrector integration scheme in the framework of a finite-element-method solution procedure. The model is implemented in the partial differential equation interface in COMSOL Multiphysics^TM.

In future work, now knowing the material parameters of ferroelectric NBT-6BT and relaxor-type NBT-6BT-4KNN, simulations of layered composites will follow. Ferroelectric/relaxor bilayers and trilayers will be considered to investigate polarization and strain coupling effects.