Variational inequalities for ferroelectric constitutive modeling

Abstract

This paper is concerned with modeling the polarization process in ferroelectric media. We develop a thermodynamically consistent model, based on phenomenological descriptions of free energy as well as switching and saturation conditions in form of inequalities. Thermodynamically consistent models naturally lead to variational formulations. We propose to use the concept of variational inequalities. We aim at combining the different phenomenological conditions into one variational inequality. In our formulation we use one Lagrange multiplier for each condition (the onset of domain switching and saturation), each satisfying Karush-Kuhn-Tucker conditions. An update for reversible and remanent quantities is then computed within one, in general nonlinear, iteration.

1. Introduction

In the current work we aim at describing the process of polarziation in ferroelectric media. Polarization has been described as a dissipative process in a thermodynamically consistent framework based on the Helmholtz free energy in a series of papers by Bassiouny et al. (1988a, 1988b) and Bassiouny and Maugin (1989a, 1989b). The notions introduced there are close to the theory of elasto-plasticity, including internal variables, yield conditions, and hardening moduli. The evolution of the internal variables is based on inequality constraints, such as the so-called switching condition. This condition defines the onset of remanent polarization in ferroelectric media.

McMeeking and Landis (2002) as well as Landis (2002) developed their theory of multi-axial polarization based on those early works. The internal variables are the polarization vector in McMeeking and Landis (2002), and an independent polarization strain is added in Landis (2002). They specified non-linear free energy functions that account also for the saturation phenomenon. Once reaching saturation, the polarization cannot grow any further, but may change direction or be reduced by (electric or mechanic) depolarization.

Kamlah and Tsakmakis (1999) chose a different set of conditions to describe the polarization process, see also the overview given in Kamlah (2001). They assume the reversible part of the free energy to be quadratic, and add separate switching and saturation conditions for remanent polarization and remanent polarization strain.

All of these frameworks have been implemented in finite element codes by various authors. Kamlah and Böhle (2001) provided the increments for the evolution polarization and polarization strain for the material described in Kamlah and Tsakmakis (1999). In Klinkel (2006) a return mapping algorithm is presented realizing Landis-type material properties. Semenov et al. (2010) provide such an algorithm for a finite element formulation including the vector (not scalar) potential of the dielectric displacement vector. Elhadrouz et al. (2005) implemented a material model including different switching functions for polarization and polarization strain, similar as proposed by Kamlah and Tsakmakis (1999). Zouari et al. (2011) proposed a quadratic element realizing a constitutive law, similar to the one proposed by Kamlah (2001).

A hybrid phenomenological model for ferroelectric ceramics is presented by Stark et al. (2016a, 2016b). Comments on the stability of ferroelectric constitutive models were given by Stark et al. (2016c) and Bottero and Idiart (2018). A general mathematical framework allowing for proofs of existence and uniqueness of solutions was recently presented by Pechstein et al. (2020). Humer et al. (2020) presented a formulation considering large deformation and hysteresis effects using threefold multiplicative decomposition of the deformation gradient tensor into a reversible mechanical, a reversible electrical, and an irreversible part.

We will follow Kamlah’s approach, but formulate the theory in the framework of variational inequalities. These inequalities arise naturally from energy-based constitutive models with dissipation, see Miehe et al. (2011) for the application to ferroelectric and also ferromagnetic materials. In our case, we add both switching and saturation condition. In the following derivations, we will see that, in the variational framework, these conditions are not independent as proposed in Kamlah (2001), but that the Lagrangian multiplier enforcing saturation becomes an additional term in the switching criterion. This way, both criteria can be satisfied at once.

The framework of variational inequalities allows the direct implementation of different inequality constraints in a thermodynamically consistent way. In contrast to return mapping algorithms the equations may be solved within one in general nonlinear iteration. An alternative approach to model the evaluation of internal variables is the construction of nonlinear dissipation functions as proposed by Landis (2002). An advantage of the latter method is, that it can be easily implemented as no inequalities have to be taken into account. A drawback of these dissipation functions is, that they are in general not differentiable.

Concerning the finite element implementation, we use mixed finite elements for the mechanical unknowns, adding independent stresses. We use non-standard TDNNS elements, where tangential displacements and the normal component of the normal stress are chosen as degrees of freedom (Pechstein and Schöberl, 2011, 2018). These elements were adapted to linear piezoelectricity in Pechstein et al. (2018) and Meindlhumer and Pechstein (2018). The TDNNS elements have been shown to be free from locking effects and therefore they are highly suitable for the efficient discretization of thin structures, which will be shown in one particular example.

Our contribution is organized as following: In the first section we provide a phenomenological material law in the spirit of Kamlah (2001). In the next section we show how this material law can be embedded into incremental variational formulations similar to Miehe et al. (2011). As switching and saturation are modeled via inequality constraints, a variational inequality is obtained. We show how these variational inequalities can be implemented in a finite element formulation. In the last section of this contribution we show several numerical results. We consider benchmark problems from the literature, including electrical polarization, mechanical depolarization, non-proportional loading and depolarization under bending stress. Additionally we show the polarization of a bimorph structure. Last, a short review is provided.

2. A phenomenological material law

In this section we briefly introduce the phenomenological material model in the spirit of Kamlah (2001), before deriving an energy-based approach in the next section. We consider a ferroelectric body in the domain $\mathrm{\Omega }\subset {\mathbb{R}}^3$ with boundary $\partial \mathrm{\Omega }$ . We are concerned with isothermal, rate independent and quasi-static deformation and polarization processes. We are interested in finding the displacement field $\mathbf{u}$ , the electric potential $\varphi$ and the remanent polarization ${\mathbf{P}}_I$ .

The electric field is introduced as the negative gradient of the electric potential $\mathbf{E}=-\nabla \varphi$ , whereas we use the linear strain tensor $\epsilon =\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T)$ . The mechanical balance equation and Gauss’ law are

$$-\mathrm{div}\mathit{\sigma }=\mathbf{f},$$

(1)

$$-\mathrm{div}\mathbf{D}=0,$$

(2)

for the stress tensor $\mathit{\sigma }$ and the dielectric displacement $\mathbf{D}$ . The mechanical boundary conditions are

$$\mathbf{u}=0\mathrm{on}{\mathrm{\Gamma }}_1\mathrm{and}\mathit{\sigma }·\mathbf{n}={\mathit{\sigma }}_0\mathrm{on}{\mathrm{\Gamma }}_2=\partial \mathrm{\Omega }\backslash {\mathrm{\Gamma }}_1,$$

(3)

and the electrical boundary conditions

$$\varphi ={\varphi }_0\mathrm{on}{\mathrm{\Gamma }}_3\mathrm{and}\mathbf{D}·\mathbf{n}=0\mathrm{on}{\mathrm{\Gamma }}_4=\partial \mathrm{\Omega }\backslash {\mathrm{\Gamma }}_3,$$

(4)

with $\mathbf{n}$ denoting the normal vector on the corresponding boundary. The constitutive model relates stress $\mathit{\sigma }$ and dielectric displacement $\mathbf{D}$ to strain $\epsilon$ and electric field $\mathbf{E}$ as well as to the remanent polarization ${\mathbf{P}}_I$ . A basic assumption is the additive decomposition of dielectric displacement and elastic strain into a reversible and a remanent or irreversible part. For the reversible parts, constitutive equations analogous to Voigt’s linear theory are assumed

$${\epsilon }_R=\epsilon -{\epsilon }_I={\mathbf{S}}^E·\mathit{\sigma }+{\mathbf{d}}^T·\mathbf{E},$$

(5)

$${\mathbf{D}}_R=\mathbf{D}-{\mathbf{P}}_I=\mathbf{d}·\mathit{\sigma }+{ϵ}^{\mathit{\sigma }}·\mathbf{E}.$$

(6)

The mechanical compliance tensor ${\mathbf{S}}^E$ as well as the dielectric tensor ${ϵ}^{\mathit{\sigma }}$ are assumed to be isotropic and constant. The components of the piezoelectric tensor $\mathbf{d}$ depend on ${\mathbf{P}}_I$ via

$$\begin{array}{l}{d}_{kji}=\frac{\left|{\mathbf{P}}_I\right|}{P_{sat}}[d_p{e}_i^P{e}_j^P{e}_k^P\\ +d_n\left({\delta }_{ij}-e_i^P{e}_j^P\right)e_k^P\\ +d_t\frac{1}{2}\left(\left({\delta }_{ki}-e_k^P{e}_i^P\right)e_j^P+\left({\delta }_{kj}-e_k^P{e}_j^P\right)e_i^P\right)]\end{array}$$

(7)

with unit vector in direction of polarization

$${\mathbf{e}}^P=\frac{{\mathbf{P}}_I}{|{\mathbf{P}}_I|},$$

(8)

and the Kronecker delta

$${\delta }_{\mathrm{ij}}=\{\begin{array}{cc}1,\hfill & \mathrm{if}i=j,\hfill \\ \multicolumn{1}{c}{0,}& \mathrm{if}i\ne j.\hfill \end{array}$$

(9)

The unpolarized material does not show any piezoelectric coupling effects. The dielectric tensor grows (linearly) with the norm of ${\mathbf{P}}_I$ . The constants $d_p$ , $d_n$ , and $d_t$ correspond to the standard parameters $d_{33}$ , $d_{31}$ , and $d_{15}$ , respectively, for the fully poled state with polarization in direction of ${\mathbf{e}}_3$ . As proposed by McMeeking and Landis (2002), we assume the remanent strain to be volume-preserving and depend only on the polarization via

$${\epsilon }_I=\frac{3}{2}\frac{S_{\mathrm{sat}}}{P_{\mathrm{sat}}^2}({\mathbf{P}}_I\otimes {\mathbf{P}}_I-\frac{1}{3}\mathbf{I}|{\mathbf{P}}_I{|}^2),$$

(10)

where $S_{\mathrm{sat}}$ characterizes the maximum possible remanent strain. The polarization strain ${\epsilon }_I$ is a deviatoric uni-axial strain and does not cause any further remanent strain. The evolution of the remanent polarization ${\mathbf{P}}_I$ is determined by switching and saturation conditions. Once the electric field is increased above the coercive field strength $E_C$ , the material will be polarized irreversibly. The switching condition is a condition on the electric driving force $\hat{\mathbf{E}}$ which will be specified in the sequel of this contribution. In the simplest case, the absolute value of $\hat{\mathbf{E}}$ may not exceed $E_C$ ,

$$f_P(\hat{\mathbf{E}})=|\hat{\mathbf{E}}|-E_C⩽0.$$

(11)

On the other hand, the saturation condition states that the remanent polarization is limited and must not exceed the saturation polarization $P_{\mathrm{sat}}$ , which resembles the fully polarized state. The condition reads

$$f_S\left({\mathbf{P}}_I\right)=|{\mathbf{P}}_I|-P_{\mathrm{sat}}⩽0.$$

(12)

3. Incremental optimization formulations

In this section we show how to embed the phenomenological material law from the previous section into an incremental optimization formulation for (coupled) irreversible processes in the spirit of Miehe et al. (2011). First, we briefly summarize the incremental formulation for purely reversible processes. This formulation will then be extended to nonlinear irreversible processes.

3.1. Reversible processes

The stored energy (or free energy) ${\mathrm{\Psi }}^R$ in the domain $\mathrm{\Omega }$ is given by

$${\mathrm{\Psi }}^R={\int }_{\mathrm{\Omega }}{\psi }^R\mathrm{d}\mathrm{\Omega },$$

(13)

with the (volumetric) energy density function ${\psi }^R$ . For coupled piezoelectric problems the energy density function ${\psi }^R$ reads

$${\psi }^R=\frac{1}{2}\epsilon :\mathit{\sigma }({\epsilon }_R,{\mathbf{D}}_R)+\frac{1}{2}\mathbf{D}·\mathbf{E}({\epsilon }_R,{\mathbf{D}}_R).$$

(14)

Now, let us consider a finite time interval $\tau :=t-t_n$ , ${\epsilon }_n$ and ${\mathbf{D}}_n$ denote strain and dielectric displacement at time $t_n$ , respectively. Further we assume body forces, boundary forces, and electric potential to be constant in the considered time interval. The work of external loads $W_{\mathrm{ext}}^{\tau }$ , splits up into two parts

$$W_{\mathrm{ext}}^{\tau }=W_{\mathrm{ext}}^{\tau ,V}+W_{\mathrm{ext}}^{\tau ,B},$$

(15)

first $W_{\mathrm{ext}}^{\tau ,V}$ , the external work according to body forces,

$$W_{\mathrm{ext}}^{\tau ,V}={\int }_{t_n}^t{\int }_{\mathrm{\Omega }}\mathbf{f}·\stackrel{·}{\mathbf{u}}\mathrm{d}\mathrm{\Omega }\mathrm{dt},$$

(16)

second $W_{\mathrm{ext}}^{\tau ,B}$ , external work according to boundary forces and applied electric potential on boundaries

$$\begin{array}{c}\multicolumn{1}{c}{W_{\mathrm{ext}}^{\tau ,B}={\int }_{t_n}^t{\int }_{{\mathrm{\Gamma }}_2}{\mathit{\sigma }}_0·\stackrel{·}{\mathbf{u}}d\mathrm{\Gamma }\mathrm{dt}}\\ +{\int }_{t_n}^t{\int }_{{\mathrm{\Gamma }}_4}{\varphi }_0\stackrel{·}{\mathbf{D}}·\mathbf{n}d\mathrm{\Gamma }\mathrm{dt}.\hfill \end{array}$$

(17)

We introduce, according to Miehe et al. (2011), the potential (for reversible processes) ${\mathrm{\Pi }}_R^{\tau }$ with its algorithmic representation

$${\mathrm{\Pi }}_R^{\tau }(\epsilon ,\mathbf{D}):={\mathrm{\Psi }}^R(\epsilon ,\mathbf{D})-{\mathrm{\Psi }}^R({\epsilon }_n,{\mathbf{D}}_n)-W_{\mathrm{ext}}^{\tau }.$$

(18)

The potential ${\mathrm{\Pi }}_R^{\tau }$ contains the difference between stored energy and work done by external loads $W_{\mathrm{ext}}^{\tau }$ in the considered time interval. As the path independence of work is an essential property of linear material behavior, $W_{\mathrm{ext}}^{\tau }$ can be determined as the difference of external works between the states at time $t_n$ and $t$ . For given external loads the associated strain and dielectric displacement is given by the constitutive minimization principle

$$\begin{array}{c}(\epsilon ,\mathbf{D})=\mathrm{Arg}\left\{\underset{\epsilon }{min}\underset{\mathbf{D}}{min}{\mathrm{\Pi }}_R^{\tau }\right\}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathrm{div}\mathbf{D}=0,\hfill \\ \multicolumn{1}{c}{\epsilon =\epsilon \left(\mathbf{u}\right).}\end{array}\hfill \end{array}$$

(19)

The solution of equation (19) is unique if the potential ${\mathrm{\Pi }}_R^{\tau }$ is convex, which is the case for stable piezoelectric materials. In the sequel of this contribution we will use, instead of strain $\epsilon$ and dielectric displacement $\mathbf{D}$ , mechanical stress $\mathit{\sigma }$ , and electric field $\mathbf{E}$ as free variables. The thermodynamic potential according to this specific choice of free variables is the free enthalpy $H^R$ given in analogy to equation (13), via the enthalpy density $h^R$ by

$$H^R={\int }_{\mathrm{\Omega }}{h}^R\mathrm{d}\mathrm{\Omega }.$$

(20)

The enthalpy density $h^R$ is related to the free energy via Legendre transformation,

$$h^R(\mathit{\sigma },\mathbf{E})=\underset{\epsilon ,\mathbf{D}}{min}({\psi }^R(\epsilon ,\mathbf{D})-\mathit{\sigma }:\epsilon -\mathbf{D}·\mathbf{E}).$$

(21)

For the case of linear piezoelasticity the enthalpy density reads

$$h^R(\mathit{\sigma },\mathbf{E})=-\frac{1}{2}\mathit{\sigma }:{\mathbf{S}}^E:\mathit{\sigma }-\mathit{\sigma }:\mathbf{d}·\mathbf{E}-\frac{1}{2}\mathbf{E}·{ϵ}^{\mathit{\sigma }}·\mathbf{E}.$$

(22)

In analogy to equation (18) we introduce for a (finite) time interval $\tau =t-t_n$ the potential ${\mathrm{\Pi }}_H^{\tau }$ with its algorithmic expression

$${\mathrm{\Pi }}_H^{\tau }=H^R(\mathit{\sigma },\mathbf{E})-H^R({\mathit{\sigma }}_n,{\mathbf{E}}_n),$$

(23)

where suffix $n$ denotes quantities at time $t_n$ . Analogously to equation (19) the associated stress and electric field is given by the constitutive maximization principle

$$\begin{array}{c}(\mathit{\sigma },\mathbf{E})=\mathrm{Arg}\left\{\underset{\mathbf{E}}{max}\underset{\mathit{\sigma }}{max}{\mathrm{\Pi }}_H^{\tau }\right\}\hfill \\ \mathrm{s}.\mathrm{t}.\{\begin{array}{c}\mathrm{div}\mathit{\sigma }=-\mathbf{f},\hfill \\ \multicolumn{1}{c}{\mathbf{E}=-\nabla \varphi .}\end{array}\hfill \end{array}$$

(24)

Note, that the potential ${\mathrm{\Pi }}_H^{\tau }$ in equation (23), contrary to ${\mathrm{\Pi }}_R^{\tau }$ in equation (18) does not contain the external work. External forces and applied potential are rather enforced via the side constraints and boundary conditions in equation (24). The first side constraint is the mechanical balance equation, the second ensures Faraday’s law by introducing the electric field as the negative gradient of the electric potential $\varphi$ .

The mechanical balance equation is typically treated by introducing the (mechanical) displacement $\mathbf{u}$ as corresponding Lagrangian multiplier. The enthalpy $H^R$ is extended and reads

$$\mathcal{L}={\int }_{\mathrm{\Omega }}{h}^R-(\mathrm{div}\mathit{\sigma }+\mathbf{f})·\mathbf{u}\mathrm{d}\mathrm{\Omega }.$$

(25)

The potential ${\mathrm{\Pi }}_{\mathcal{L}}^{\tau }$ corresponding to the extended enthalpy $\mathcal{L}$ is given by

$${\mathrm{\Pi }}_{\mathcal{L}}^{\tau }=\mathcal{L}(\mathit{\sigma },\mathbf{E},\mathbf{u})-\mathcal{L}({\mathit{\sigma }}_{\mathbf{n}},{\mathbf{E}}_n,{\mathbf{u}}_n).$$

(26)

The constitutive maximization problem (equation (24)) transforms into a saddle point problem. It reads

$$(\mathit{\sigma },\mathbf{E},\mathbf{u})=\mathrm{Arg}\left\{\underset{\mathbf{E}=-\nabla \varphi }{max}\underset{\mathit{\sigma }}{max}\underset{\mathbf{u}}{min}{\mathrm{\Pi }}_{\mathcal{L}}^{\tau }\right\}.$$

(27)

3.2. Dissipative processes

The material response in ferroelastic materials—for significantly high external loads—is characterized by non-reversible local response. This contribution focuses on polarization effects. The process of polarization is a dissipative process. Conducted work is not stored completely as free energy, but some of it is dissipated into heat.

We introduce the dissipated work $D^{\tau }$ for the time interval $\tau =t_n-t$ and the considered region $\mathrm{\Omega }$ . It is given by

$$D^{\tau }:={\int }_{t_n}^t{\int }_{\mathrm{\Omega }}\mathcal{D}\mathrm{d}\mathrm{\Omega }\mathrm{dt}⩾0,$$

(28)

with the local dissipation $\mathcal{D}$ . The second law of thermodynamics states that, for each time interval, $D^{\tau }⩾0$ has to be non-negative. This implies positive dissipation $\mathcal{D}⩾0$ for arbitrary processes. In this contribution we follow the approach of internal variables, and introduce irreversible polarization ${\mathbf{P}}_I$ . Consequently the stored energy density function for dissipative processes $\psi$ depends also on ${\mathbf{P}}_I$ ,

$$\psi =\psi (\epsilon ,\mathbf{D},{\mathbf{P}}_I).$$

(29)

Analogously to equation (13) the stored energy in the domain $\mathrm{\Omega }$ is given by

$$\mathrm{\Psi }={\int }_{\mathrm{\Omega }}\psi \mathrm{d}\mathrm{\Omega }.$$

(30)

Note, that in this work, remanent straining is always assumed to depend directly on the remanent polarization as proposed by for example, McMeeking and Landis (2002) or Linnemann et al. (2009). Therefore, the irreversible part of the strain ${\epsilon }_I$ is not to be considered as a free variable.

According to Miehe et al. (2011) we introduce the driving force $\hat{\mathbf{E}}$ . It is defined as the negative derivative of the stored energy density $\psi$ by the internal variables, in our case the irreversible polarization ${\mathbf{P}}_I$

$$\hat{\mathbf{E}}:=-\frac{\partial }{\partial {\mathbf{P}}_I}\psi .$$

(31)

The driving force $\hat{\mathbf{E}}$ is also referred to as internal constitutive force, as $\hat{\mathbf{E}}$ is the work conjugate of ${\mathbf{P}}_I$ . As shown in Miehe et al. (2011), the dissipation is related to the driving force and the evolution of the internal variables via

$$\mathcal{D}=\hat{\mathbf{E}}·{\stackrel{·}{\mathbf{P}}}_I⩾0.$$

(32)

The constitutive minimization principle (equation (19)) may be extended to the case of irreversible processes. We introduce, according to Miehe et al. (2011) for the case of dissipative (or irreversible) processes, the potential ${\mathrm{\Pi }}_D^{\tau }$ with its algorithmic expression

$$\begin{array}{c}{\mathrm{\Pi }}_D^{\tau }=\mathrm{\Psi }(\epsilon ,\mathbf{D},{\mathbf{P}}_I)-\mathrm{\Psi }({\epsilon }_n,{\mathbf{D}}_n,{\mathbf{P}}_{I,n})\hfill \\ +D^{\tau }-W_{\mathrm{ext}}^{\tau }.\hfill \end{array}$$

(33)

During the time interval $\tau$ the energy stored is $\mathrm{\Psi }(\epsilon ,\mathbf{D},{\mathbf{P}}_I)-\mathrm{\Psi }({\epsilon }_n,{\mathbf{D}}_n,{\mathbf{P}}_{I,n})$ and $D^{\tau }$ is dissipated.

Analogously to equation (19), for given external loads and states at time $t_n$ , the strain, dielectric displacement and irreversible polarization satisfy the constitutive minimization principle which reads for dissipative processes

$$\begin{array}{c}(\epsilon ,\mathbf{D},{\mathbf{P}}_I)=\mathrm{Arg}\left\{\underset{\epsilon }{min}\underset{\mathbf{D}}{min}\underset{{\mathbf{P}}_I}{min}{\mathrm{\Pi }}_D^{\tau }\right\}\hfill \\ s.t.\{\begin{array}{c}{f}_P(\hat{\mathbf{E}})⩽0,\hfill \\ \multicolumn{1}{c}{f_S({\mathbf{P}}_I)⩽0,}\\ \multicolumn{1}{c}{\mathrm{div}\mathbf{D}=0.}\end{array}\hfill \end{array}$$

(34)

Note that, contrary to equation (19), the minimization prinziple (equation (34)) for dissipative processes is constrained by inequality constraints. This has to be taken into account when solving for the free variables.

As in the reversible case, an equivalent enthalpy-based formulation can be found. We define the total enthalpy $H$ as in equations (20) and (21) but using the energy density $\psi$ ,

$$h(\mathit{\sigma },\mathbf{E},{\mathbf{P}}_{\mathbf{I}})=\underset{\epsilon ,\mathbf{D}}{min}(\psi (\epsilon ,\mathbf{D},{\mathbf{P}}_I)-\mathit{\sigma }:\mathbf{\epsilon }-\mathbf{D}·\mathbf{E}).$$

(35)

With the Lagrangian $\mathcal{L}$ and its potential ${\mathrm{\Pi }}_{\mathcal{L}}^{\tau }$ defined in analogy to equations (25) and (26) the following saddle point problem is obtained

$$\begin{array}{c}(\mathit{\sigma },\mathbf{E},\mathbf{u},{\mathbf{P}}_{\mathbf{I}})=\mathrm{Arg}\left\{\underset{\mathbf{E}=-\nabla \varphi }{max}\underset{\mathit{\sigma }}{max}\underset{\mathbf{u}}{min}\underset{{\mathbf{P}}_I}{min}{\mathrm{\Pi }}_{\mathcal{L}}^{\tau }\right\}\hfill \\ s.t.\{\begin{array}{c}{f}_P(\hat{\mathbf{E}})⩽0,\hfill \\ \multicolumn{1}{c}{f_S({\mathbf{P}}_I)⩽0.}\end{array}\hfill \end{array}$$

(36)

4. Incremental formulation of ferroelectricity

The focus of this section will be the realization of the energy density $\psi$ and the reformulation of the switching constraint $f_P(\hat{\mathbf{E}})⩽0$ . The energy density will be chosen such that it reflects the constitutive equations (5), (6) as well as the saturation condition (12).

We provide energy densities matching the constitutive equations from the previous section as a basis for all further deductions. According to the literature (see e.g. Kamlah, 2001; Landis, 2002; Tich, et al., 2010), we follow the approach of additive decomposition of free energy. It decomposes into two parts,

$$\psi ={\psi }^R+{\psi }^I,$$

(37)

with the reversible or stored part of the free energy ${\psi }^R$ and the additional contribution of the irreversible quantities to the free energy ${\psi }^I$ . The reversible part ${\psi }^R$ is given by

$$\begin{array}{c}{\psi }^R=\frac{1}{2}{\mathbf{D}}_R·{\left({ϵ}^{\epsilon }\right)}^{-1}·{\mathbf{D}}_R+{\mathbf{D}}_R·\mathbf{h}:{\epsilon }_R\hfill \\ +\frac{1}{2}{\epsilon }_R:{\mathbf{C}}^{\mathbf{D}}:{\epsilon }_R,\hfill \end{array}$$

(38)

with $\mathbf{h}$ , the tensor of piezoelectric constants, ${ϵ}^{\epsilon }$ the dielectric tensor at constant strain and ${\mathbf{C}}^{\mathbf{D}}$ the tensor of mechanical stiffness at constant dielectric displacement. These tensors are related to the material tensors in equation (5) in the standard way.

$${\mathbf{S}}^{\mathbf{D}}({\mathbf{P}}_I)={\mathbf{S}}^{\mathbf{E}}-\mathbf{d}({\mathbf{P}}_I{)}^T·{\beta }^{\mathit{\sigma }}·\mathbf{d}({\mathbf{P}}_I),$$

(39)

$${\beta }^{\epsilon }={\beta }^{\mathit{\sigma }}+{\beta }^{\mathit{\sigma }}·\mathbf{d}({\mathbf{P}}_I):{\mathbf{S}}^{\mathbf{D}}({\mathbf{P}}_I):\mathbf{d}({\mathbf{P}}_I{)}^T·{\beta }^{\mathit{\sigma }},$$

(40)

$$\mathbf{h}({\mathbf{P}}_I)={\beta }^{\epsilon }·\mathbf{d}({\mathbf{P}}_I):{\mathbf{C}}^{\mathbf{E}},$$

(41)

$${\mathbf{C}}^{\mathbf{D}}({\mathbf{P}}_I)={\mathbf{C}}^{\mathbf{E}}+{\mathbf{C}}^{\mathbf{E}}:\mathbf{d}({\mathbf{P}}_I)·{\beta }^{\epsilon }·\mathbf{d}({\mathbf{P}}_I{)}^T:{\mathbf{C}}^{\mathbf{E}},$$

(42)

with ${\beta }^{\mathit{\sigma }}={({ϵ}^{\mathit{\sigma }})}^{-1}$ and ${\beta }^{\epsilon }=({ϵ}^{\epsilon }{)}^{-1}$ and the tensor of mechanical compliance at constant dielectric displacement ${\mathbf{S}}^{\mathbf{D}}$ . This implies a non-trivial dependence on ${\mathbf{P}}_I$ for ${ϵ}^{\epsilon }$ , ${\mathbf{C}}^{\mathbf{D}}$ , and $\mathbf{h}$ . Note that, in the enthalpy-based formulations proposed in this work only much simpler structured material moduli ${\epsilon }^{\mathit{\sigma }}$ , ${\mathbf{S}}^{\mathbf{E}}$ , and $\mathbf{d}$ are needed. Of course non-trivial dependences on ${\mathbf{P}}_I$ may be avoided when assuming ${ϵ}^{\epsilon }$ and ${\mathbf{C}}^{\mathbf{D}}$ to be constant and isotropic.

The additional contribution of the free energy is given by

$${\psi }^I=\frac{1}{2}c{\mathbf{P}}_I·{\mathbf{P}}_I,$$

(43)

in the simplest case, where $c$ is a hardening parameter.

An additional term representing saturation is added in the following. The saturation condition $f_S⩽0$ (cmp. (12)) is reformulated as a variational inequality, involving the Lagrangian multiplier ${\lambda }_S$ ,

$${\lambda }_S{f}_S({\mathbf{P}}_I)⩽0\forall {\lambda }_S⩾0.$$

(44)

To ensure the saturation condition, we add the following supremum term to the energy function:

$$\psi ={\psi }^R+{\psi }^I+\underset{{\lambda }_S⩾0}{sup}{\lambda }_S{f}_S({\mathbf{P}}_I).$$

(45)

For states of irreversible polarization ${\mathbf{P}}_I$ not violating the saturation condition ( $f_S⩽0$ ), the supremum in equation (45) takes zero value, and $\psi$ is not altered in comparison to equation (37). Otherwise, if $|{\mathbf{P}}_I|$ exceeds the saturation polarization, the saturation condition is violated ( $f_S>0$ ), and therefore $\underset{{\lambda }_S⩾0}{sup}{\lambda }_S{f}_S$ tends to infinity, and with it the energy function $\psi$ . Note, that for all admissible values of ${\mathbf{P}}_I$ —not violating the saturation condition— $\psi$ takes finite values, while for inadmissible (improper) values of ${\mathbf{P}}_I$ , the free energy function $\psi$ tends to infinity. In the spirit of optimization problems, this leads to solutions with admissible states of irreversible polarization ${\mathbf{P}}_I$ .

For further deductions we treat the Lagrange multiplier ${\lambda }_S$ as a free variable. This is indicated by the index $S$ , as we use

$${\psi }^S={\psi }^R+{\psi }^I+{\lambda }_S{f}_S$$

(46)

In the final optimization problem ${\lambda }_S$ will be maximized.

4.1. Transformation to enthalpy

The enthalpy is a function of the free variables electric field $\mathbf{E}$ , mechanical stress $\mathit{\sigma }$ and irreversible polarization ${\mathbf{P}}_I$ and is related to the free energy via Legendre transformation (equation (35)). The enthalpy density $h^S$ including saturation condition (equation (46)) reads

$$\begin{array}{c}{h}^S=-\frac{1}{2}\mathbf{E}·{ϵ}^{\mathit{\sigma }}·\mathbf{E}-\frac{1}{2}\mathit{\sigma }:{\mathbf{S}}^{\mathbf{E}}:\mathit{\sigma }-\mathit{\sigma }:\mathbf{d}·\mathbf{E}\hfill \\ +\frac{1}{2}c{\mathbf{P}}_I·{\mathbf{P}}_I-{\mathbf{P}}_I·\mathbf{E}-{\epsilon }_I:\mathit{\sigma }\hfill \\ +{\lambda }_S{f}_S\left({\mathbf{P}}_I\right).\hfill \end{array}$$

(47)

In the enthalpy setting the driving force is defined as the negative derivative of the total enthalpy $h^S$ by the irreversible quantities, here by the irreversible polarization ${\mathbf{P}}_I$ . In contrast to the original model by Kamlah, the driving force $\hat{\mathbf{E}}$ contains saturation via the corresponding Lagrangian multiplier ${\lambda }_S$ . The driving force reads

$$\begin{array}{c}\hat{\mathbf{E}}=-\frac{\partial h^S}{\partial {\mathbf{P}}_I}=\mathbf{E}+\mathit{\sigma }:\left(\frac{\mathbf{d}}{\mathrm{d}{\mathbf{P}}_{\mathbf{I}}}\mathbf{d}\right)·\mathbf{E}\hfill \\ +\mathit{\sigma }:\left(\frac{\mathbf{d}}{\mathrm{d}{\mathbf{P}}_{\mathbf{I}}}{\mathbf{\epsilon }}_{\mathbf{I}}\right)-\mathbf{c}{\mathbf{P}}_{\mathbf{I}}-\frac{{\mathrm{df}}_{\mathbf{S}}\mathbf{(}{\mathbf{P}}_{\mathbf{I}}\mathbf{)}}{\mathrm{d}{\mathbf{P}}_{\mathbf{I}}}{\mathbf{\lambda }}_{\mathbf{S}}.\hfill \end{array}$$

(48)

Note that equation (48) implies that via equation (32) the saturation is part of the dissipation $\mathcal{D}$ .

4.2. Incremental optimization principle

Again we consider the (finite) time interval $\tau =t-t_n$ , with suffix $n$ denoting quantities at time $t_n$ . For the modified enthalpy $h^S$ we proceed to a Lagrangian including the equilibrium condition in the same way as equation (24)

$${\mathcal{L}}^S={\int }_{\mathrm{\Omega }}{h}^S-(\mathrm{div}\mathit{\sigma }+\mathbf{f})·\mathbf{u}\mathrm{d}\mathrm{\Omega }.$$

(49)

The algorithmic expression of the corresponding potential ${\mathrm{\Pi }}_S^{\tau }$ reads

$$\begin{array}{c}{\mathrm{\Pi }}_S^{\tau }={\mathcal{L}}^S(\mathit{\sigma },\mathbf{E},\mathbf{u},{\mathbf{P}}_{\mathbf{I}},{\mathbf{\lambda }}_{\mathbf{S}}\mathbf{)}\hfill \\ -{\mathcal{L}}^S({\mathit{\sigma }}_{\mathbf{n}},{\mathbf{E}}_n,{\mathbf{u}}_n,{\mathbf{P}}_{I,n},{\lambda }_{S,n})+D^{\tau }.\hfill \end{array}$$

(50)

The saddle point problem reads

$$\begin{array}{c}(\mathbf{E},\mathit{\sigma },\mathbf{u},{\mathbf{P}}_{\mathbf{I}},{\mathbf{\lambda }}_{\mathbf{S}}\mathbf{)}=\hfill \\ \mathrm{Arg}\left\{\underset{\mathbf{E}=-\nabla \mathbf{\varphi }}{max}\underset{\mathit{\sigma }}{max}\underset{\mathbf{u}}{min}\underset{{\mathbf{P}}_{\mathbf{I}}}{min}\underset{{\mathbf{\lambda }}_{\mathbf{S}}⩾0}{max}{\mathbf{\Pi }}_{\mathbf{S}}^{\mathbf{\tau }}\right\}\hfill \\ \mathrm{s}.\mathrm{t}.{\mathbf{f}}_{\mathbf{P}}\left(\hat{\mathbf{E}}\right)⩽0.\hfill \end{array}$$

(51)

In order to consider the switching condition, another Lagrangian multiplier ${\lambda }_P$ is introduced in analogy to equation (44) via

$${\lambda }_P{f}_P\left(\hat{\mathbf{E}}\right)⩽0\forall {\lambda }_P⩾0.$$

(52)

Consequently the Lagrangian ${\mathcal{L}}^S$ and the corresponding potential ${\mathrm{\Pi }}_S^{\tau }$ have to be updated by the terms in equation (52). The Lagrangian ${\mathcal{L}}^P$ , involving polarization, reads

$${\mathcal{L}}^P={\int }_{\mathrm{\Omega }}{h}^S+(\mathrm{div}\mathit{\sigma }+\mathbf{f})·\mathbf{u}+{\lambda }_P{f}_P\mathrm{d}\mathrm{\Omega }.$$

(53)

The algebraic form of the corresponding potential ${\mathrm{\Pi }}_P^{\tau }$ is given by

$$\begin{array}{c}{\mathrm{\Pi }}_P^{\tau }={\mathcal{L}}^P(\mathit{\sigma },\mathbf{E},\mathbf{u},{\mathbf{P}}_{\mathbf{I}},{\mathbf{\lambda }}_{\mathbf{S}},{\mathbf{\lambda }}_{\mathbf{P}}\mathbf{)}-\hfill \\ {\mathcal{L}}^P({\mathit{\sigma }}_{\mathbf{n}},{\mathbf{E}}_n,{\mathbf{u}}_n,{\mathbf{P}}_{I,n},{\lambda }_S,{\lambda }_{P,n})+D_{\tau },\hfill \end{array}$$

(54)

the corresponding saddle point problem, extended by the (independent) variable ${\lambda }_P$ , reads

$$\begin{array}{c}(\mathbf{u},\mathbf{E},\mathit{\sigma },{\mathbf{P}}_{\mathbf{I}},{\mathbf{\lambda }}_{\mathbf{P}},{\mathbf{\lambda }}_{\mathbf{S}})=\hfill \\ \mathrm{Arg}\left\{\underset{\mathbf{u}}{min}\underset{\mathbf{E}=-\nabla \varphi }{max}\underset{\mathit{\sigma }}{max}\underset{{\mathbf{P}}_I}{min}\underset{{\lambda }_P⩾0}{max}\underset{{\lambda }_S⩾0}{max}{\mathrm{\Pi }}_P^{\tau }\right\}.\hfill \end{array}$$

(55)

4.3. Variation of optimization problem

The saddle point problem (equation (55)) may be solved by the method of variation. As the Lagrangian multipliers ${\lambda }_S$ and ${\lambda }_P$ are constrained to be non-negative, the variation of the potential ${\mathrm{\Pi }}_P^{\tau }$ leads to a variational inequality. Variation with respect to $\mathbf{u}$ , $\mathbf{E}$ , $\mathit{\sigma }$ , and ${\mathbf{P}}_I$ can be done in the standard way. Arbitrary virtual values $\delta \mathbf{u}$ , $\delta \mathit{\sigma }$ , $\delta \mathbf{E}$ , and $\delta {\mathbf{P}}_I$ are admissible, only the essential (Dirichlet) boundary conditions on the corresponding boundaries ${\mathrm{\Gamma }}_1$ , ${\mathrm{\Gamma }}_2$ , ${\mathrm{\Gamma }}_3$ , and ${\mathrm{\Gamma }}_4$ have to be fulfilled, respectively. A variational equation is obtained via

$$\frac{\partial {\mathrm{\Pi }}_P^{\tau }}{\partial \mathbf{E}}·\delta \mathbf{E}+\frac{\partial {\mathrm{\Pi }}_P^{\tau }}{\partial \mathit{\sigma }}:\delta \mathit{\sigma }+\frac{\partial {\mathbf{\Pi }}_{\mathbf{P}}^{\mathbf{\tau }}}{\partial {\mathbf{P}}_{\mathbf{I}}}·\mathbf{\delta }{\mathbf{P}}_{\mathbf{I}}+\frac{\partial {\mathbf{\Pi }}_{\mathbf{P}}^{\mathbf{\tau }}}{\partial \mathbf{u}}·\mathbf{\delta }\mathbf{u}=0.$$

(56)

Note, that external (body) forces are already considered within the potential ${\mathrm{\Pi }}_P^{\tau }$ . The saturation and polarization condition, which are given as inequalities, cannot be treated in the standard way. As the set of admissible ${\lambda }_P$ and ${\lambda }_S$ is restricted by inequalities $\delta {\lambda }_P$ and $\delta {\lambda }_S$ are not arbitrary. These conditions directly lead to variational inequalities, we refer to the monograph (Han and Reddy, 1999) for an introduction into variational inequalities in the framework of elasto-plasticity. For both, ${\lambda }_P$ and ${\lambda }_S$ a test function ${\overline{\lambda }}_P$ and ${\overline{\lambda }}_S$ is introduced. The potential ${\mathrm{\Pi }}_P^{\tau }$ is maximized with respect to ${\lambda }_P$ if and only if for all other admissible choices ${\overline{\lambda }}_P⩾0$ there holds

$${\overline{\lambda }}_P{f}_P(\hat{\mathbf{E}})⩽{\lambda }_P{f}_P(\hat{\mathbf{E}})\forall {\overline{\lambda }}_P⩾0.$$

(57)

For the saturation and the corresponding Lagrangian multiplier ${\lambda }_S$ , the situation is a bit more involved, as the driving force $\hat{\mathbf{E}}$ depends on ${\lambda }_S$ . Writing

$$\stackrel{\overline{^}}{\mathbf{E}}=\hat{\mathbf{E}}({\overline{\lambda }}_S)$$

(58)

for some admissible test function ${\overline{\lambda }}_S⩾0$ , one finds that ${\mathrm{\Pi }}_P^{\tau }$ is maximized with respect to ${\overline{\lambda }}_S$ if and only if

$$\begin{array}{c}\stackrel{\overline{^}}{\mathbf{E}}·{\stackrel{·}{\mathbf{P}}}_I+{\overline{\lambda }}_S{f}_S\left({\mathbf{P}}_I\right)+{\lambda }_P{f}_P(\stackrel{\overline{^}}{\mathbf{E}})⩽\hfill \\ \multicolumn{1}{c}{\hat{\mathbf{E}}·{\stackrel{·}{\mathbf{P}}}_I+{\lambda }_S{f}_S\left({\mathbf{P}}_I\right)+{\lambda }_P{f}_P(\hat{\mathbf{E}})}\\ \multicolumn{1}{c}{\forall {\overline{\lambda }}_S⩾0.}\end{array}$$

(59)

Inserting the definition of $\hat{\mathbf{E}}$ (equation (48)), which implies that

$$\stackrel{\overline{^}}{\mathbf{E}}-\hat{\mathbf{E}}=-({\overline{\lambda }}_S-{\lambda }_S)\frac{{\mathrm{df}}_S}{\mathrm{d}{\mathbf{P}}_I},$$

(60)

the above inequality can further be reduced to

$$\begin{array}{c}({\overline{\lambda }}_S-{\lambda }_S)(-\frac{{\mathrm{df}}_S}{\mathrm{d}{\mathbf{P}}_I}·{\stackrel{·}{\mathbf{P}}}_I+f_S-{\lambda }_P\frac{{\mathrm{df}}_P}{\mathrm{d}\hat{\mathbf{E}}}\frac{{\mathrm{df}}_S}{\mathrm{d}{\mathbf{P}}_I})⩽0\hfill \\ \multicolumn{1}{c}{\forall {\overline{\lambda }}_S⩾0.}\end{array}$$

(61)

Summing up, the variational inequality to hold for all admissible $\delta \mathbf{E}$ , $\delta \mathit{\sigma }$ , $\delta {\mathbf{P}}_I$ , and $\delta \mathbf{u}$ as well as for all test functions ${\overline{\lambda }}_S⩾0$ and ${\overline{\lambda }}_P⩾0$ reads

$$\begin{array}{c}\frac{\partial {\mathrm{\Pi }}_P^{\tau }}{\partial \mathbf{E}}·\delta \mathbf{E}++\frac{\partial {\mathrm{\Pi }}_P^{\tau }}{\partial \mathit{\sigma }}:\delta \mathit{\sigma }+\hfill \\ \frac{\partial {\mathbf{\Pi }}_{\mathbf{P}}^{\mathbf{\tau }}}{\partial {\mathbf{P}}_{\mathbf{I}}}·\mathbf{\delta }{\mathbf{P}}_{\mathbf{I}}+\frac{\partial {\mathbf{\Pi }}_{\mathbf{P}}^{\mathbf{\tau }}}{\partial \mathbf{u}}·\mathbf{\delta }\mathbf{u}+\hfill \\ \mathbf{(}{\overline{\mathbf{\lambda }}}_{\mathbf{S}}-{\mathbf{\lambda }}_{\mathbf{S}}\mathbf{)}({\mathbf{f}}_{\mathbf{S}}-\frac{{\mathrm{df}}_{\mathbf{S}}}{\mathrm{d}{\mathbf{P}}_{\mathbf{I}}}·\mathbf{(}{\mathbf{\lambda }}_{\mathbf{P}}\frac{{\mathrm{df}}_{\mathbf{P}}}{\mathrm{d}\hat{\mathbf{E}}}+{\stackrel{·}{\mathbf{P}}}_{\mathbf{I}}\mathbf{)})+\hfill \\ \mathbf{(}{\overline{\mathbf{\lambda }}}_{\mathbf{P}}-{\mathbf{\lambda }}_{\mathbf{P}}\mathbf{)}{\mathbf{f}}_{\mathbf{P}}⩽0.\hfill \end{array}$$

(62)

4.4. Interpretation of the Lagrangian multipliers

In this section we explain the meaning of the Lagrangian multipliers ${\lambda }_P$ and ${\lambda }_S$ . For reasons of simplicity we restrict ourselves to the uncoupled, purely electrical problem. Note, that the conclusions of this section may be extended directly to the fully coupled problem, as the Lagrangian multipliers ${\lambda }_P$ and ${\lambda }_S$ , as well as the corresponding inequalities read the same as for the mechanical problem with stress states $\mathit{\sigma }=0$ .

For this electric case, as all mechanical quantities vanish, the driving force reduces to the simple form

$$\hat{\mathbf{E}}=\mathbf{E}-c{\mathbf{P}}_I-{\lambda }_S\frac{{\mathbf{P}}_I}{|{\mathbf{P}}_I|}.$$

(63)

The variational inequality (equation (62)) reads

$$\begin{array}{c}\delta \left[-\frac{1}{2}\mathbf{E}\cdot {ϵ}^{\mathit{\sigma }}\cdot \mathbf{E}-\mathbf{E}\cdot {\mathbf{P}}_I+\frac{c}{2}{\mathbf{P}}_I\cdot {\mathbf{P}}_I+\widehat{\mathbf{E}}\cdot {\stackrel{·}{\mathbf{P}}}_I\right]+\\ {\lambda }_P\left(\frac{d{f}_P}{d\widehat{\mathbf{E}}}\frac{\partial \widehat{\mathbf{E}}}{\partial \mathbf{E}}\delta \mathbf{E}+\frac{\partial f_P}{\partial {\mathbf{P}}_I}\delta {\mathbf{P}}_I\right)+\\ ({\overline{\lambda }}_S-{\lambda }_S)\left(f_S-\frac{d{f}_S}{d{\mathbf{P}}_I}\cdot ({\lambda }_P\frac{d{f}_P}{d\widehat{\mathbf{E}}}+{\stackrel{·}{\mathbf{P}}}_I)\right)+\\ ({\overline{\lambda }}_P-{\lambda }_P)f_P⩽0.\end{array}$$

(64)

We consider an update step with the (updated) solutions $\mathbf{E}$ , ${\mathbf{P}}_I$ , ${\lambda }_P$ , and ${\lambda }_S$ . We introduce $\mathrm{\Delta }{\mathbf{P}}_I$ as the update of the irreversible polarization between two calculation steps. As the considered process is rate independent and quasi-static the dissipation is given via $\hat{\mathbf{E}}·\mathrm{\Delta }{\mathbf{P}}_I$ . Splitting up (equation (64)) into individual equations and inequalities for the single variational quantities we get

$$\delta \mathbf{E}\left(-{ϵ}^{\mathit{\sigma }}\cdot \mathbf{E}-{\mathbf{P}}_I+\frac{\partial \widehat{\mathbf{E}}}{\partial \mathbf{E}}\mathrm{\Delta }{\mathbf{P}}_I+{\lambda }_P\frac{d{f}_P}{d\widehat{\mathbf{E}}}\frac{\partial \widehat{\mathbf{E}}}{\partial \mathbf{E}}\right)=0$$

(65)

$$\begin{array}{c}\delta {\mathbf{P}}_I(-\mathbf{E}+c{\mathbf{P}}_I+\hat{\mathbf{E}}+\frac{\partial \hat{\mathbf{E}}}{\partial {\mathbf{P}}_I}\mathrm{\Delta }{\mathbf{P}}_I+\hfill \\ {\lambda }_P\frac{{\mathrm{df}}_P}{\mathrm{d}\hat{\mathbf{E}}}\frac{\partial \hat{\mathbf{E}}}{\partial {\mathbf{P}}_I}+{\lambda }_S\frac{{\mathrm{df}}_S}{\mathrm{d}{\mathbf{P}}_I})=0\hfill \end{array}$$

(66)

$$({\overline{\lambda }}_P-{\lambda }_P)f_P⩽0$$

(67)

$$({\overline{\lambda }}_S-{\lambda }_S)(f_S-\frac{{\mathrm{df}}_S}{\mathrm{d}{\mathbf{P}}_I}·({\lambda }_P\frac{{\mathrm{df}}_P}{\mathrm{d}\hat{\mathbf{E}}}+\mathrm{\Delta }{\mathbf{P}}_I))⩽0.$$

(68)

To get an interpretation of ${\lambda }_P$ , consider a state, at which the polarization process has started, such that $f_P=0$ , $\hat{\mathbf{E}}\ne 0$ , and ${\lambda }_P>0$ non zero. The irreversible polarization is further considered to be (and to stay) below the saturation level. In this case, the derivative of the switching condition is given by

$$\frac{{\mathrm{df}}_P}{\mathrm{d}\hat{\mathbf{E}}}=\hat{\mathbf{E}}/|\hat{\mathbf{E}}|.$$

(69)

From the variation of the irreversible polarization $\delta {\mathbf{P}}_I$ in equation (66), taking into account the definition (equation (63)) of $\hat{\mathbf{E}}$ and (equation (69)) we get a relation between ${\lambda }_P$ and $\mathrm{\Delta }{\mathbf{P}}_I$ via

$$\frac{\partial \hat{\mathbf{E}}}{\partial {\mathbf{P}}_I}(\mathrm{\Delta }{\mathbf{P}}_I+{\lambda }_P\frac{\hat{\mathbf{E}}}{|\hat{\mathbf{E}}|})=0.$$

(70)

From equation (70) one deduces that the polarization update $\mathrm{\Delta }{\mathbf{P}}_I$ is in direction of $\hat{\mathbf{E}}$ . As $\frac{\hat{\mathbf{E}}}{|\hat{\mathbf{E}}|}$ is a unit vector the norm of $\mathrm{\Delta }{\mathbf{P}}_I$ equals ${\lambda }_P$ . On the other hand, if the coercive field is not reached and $f_P<0$ , equation (67) ensures that ${\lambda }_P=0$ , and further $\mathrm{\Delta }{\mathbf{P}}_I=0$ .

Before we give an interpretation for ${\lambda }_S$ , we first show that equation (68) really ensures the saturation condition. In case the coercive electric field is reached and $f_P=0$ , equation (70) holds, and thereby the last term in equation (68) reduces to zero, leaving KKT conditions for the saturation condition only. On the other hand, if the coercive field is not reached and $f_P<0$ , we find ${\lambda }_P=0$ and $\mathrm{\Delta }{\mathbf{P}}_I=0$ , which again reduces the last term in equation (68) to zero. Thus the saturation condition has to be satisfied in any case.

For the interpretation of ${\lambda }_S$ we consider a state where $f_S=0$ , that is, saturation polarization is reached. The driving electric field $\hat{\mathbf{E}}$ depending on ${\lambda }_S$ still has to satisfy the switching condition,

$$f_P=||\mathbf{E}-c{\mathbf{P}}_I-{\lambda }_S\frac{{\mathbf{P}}_I}{|{\mathbf{P}}_I|}||-E_C⩽0.$$

(71)

The Lagrangian multiplier ${\lambda }_S$ grows with the electric field $\mathbf{E}$ as soon as further growth of the remanent polarization is prohibited by the saturation condition. A graphical interpretation of ${\lambda }_P$ and ${\lambda }_S$ can be found in Figures 1 and 2, respectively.

Figure 1.

Graphical interpretation of Lagrange multipliers for switching condition.

Figure 2.

Graphical interpretation of Lagrange multipliers for saturation condition.

5. Finite element method

In this section, the proposed finite element discretization is described. For the electric potential continuous, nodal finite elements are used, while the polarization vector is approximated constant on each element. Also, the Lagrangian multipliers are realized taking one value per element. As already done in the derivations of the previous sections, displacement and stress are considered independent unknowns. To do so, a mixed finite element scheme is used. The specific approach taken in this work uses tangential displacements and normal components of the normal stress vector $(\mathit{\sigma }·\mathbf{n}\mathbf{)}·\mathbf{n}$ as degrees of freedom, which motivate the abbreviation “TDNNS.” This method was originally developed for elastic solids by Pechstein and Schöberl (2011, 2012) and later extended for linear piezoelectric materials (Meindlhumer and Pechstein, 2018; Pechstein et al., 2018) and for and for geometrically nonlinear electro-mechanically coupled problems (Pechstein, 2019). As the TDNNS method does not suffer from locking effects for elements of arbitrary aspect ratios, it is highly suitable for the discretization of thin structures.

This work is based on Pechstein et al. (2018) and Meindlhumer and Pechstein (2018) for linear piezoelectric materials. Note that, for the mixed TDNNS method, the $\mathbf{d}$ -tensor formulation can be used directly, and there is no need to transfer the electric permittivity at constant stress ${ϵ}^{\mathit{\sigma }}$ to that at constant strain ${ϵ}^{\epsilon }$ algebraically, as is necessary in standard formulations. As the underlying mathematics of the TDNNS method are rather involved, we skip it here, and only mention that work pairs such as ${\int }_{\mathrm{\Omega }}\mathit{\sigma }:\mathbf{\epsilon }$ or ${\int }_{\mathrm{\Omega }}\mathrm{div}\mathit{\sigma }·\mathbf{u}$ need to be considered in distributional sense and involve additional integration on element boundaries.

To solve the variational inequality, an active-set strategy is proposed. This means, in each iterative step, two active sets are identified beforehand: the switching active set $A_P$ where the switching condition shall be enforced as $f_P=0$ , and the saturation active set $A_S$ where the polarization condition is enforced as $f_S=0$ . On the other hand, for all elements not in the active set, the Lagrangian multipliers ${\lambda }_P$ and ${\lambda }_S$ are set to zero, respectively. For these fixed active sets, the nonlinear equations are solved by Newton’s method, with ${\lambda }_S$ and ${\lambda }_P$ free in their respective active sets. On convergence, the Kuhn-Tucker compatibility conditions are checked. All elements where either saturation or switching condition are violated are added to the respective active set. On the other hand, all elements where a Lagrangian multiplier is found negative are removed from the respective active set. An according algorithm can be found in Algorithm 1.

Algorithm 1: Active set strategy for solving the variational inequality.

Data: values for ${\mathbf{u}}_0,{\mathit{\sigma }}_0,{\varphi }_0,{\mathbf{P}}_{I0},{\lambda }_{S0},{\lambda }_{P0}$ and active sets $A_{S0},A_{P0}$ from the last converged time step Result: values for $\mathbf{u},\mathit{\sigma },\mathbf{\varphi },{\mathbf{P}}_{\mathbf{I}},{\mathbf{\lambda }}_{\mathbf{S}},{\mathbf{\lambda }}_{\mathbf{P}}$ , active sets $A_S,A_P$ initialize unknows with initial data, active sets $A_S=A_{S0},A_P=A_{P0}$ ; whilethe active sets change do  solve non-linear problem with ${\lambda }_S=0$ in $\mathrm{\Omega }\backslash A_S$ ,   ${\lambda }_P=0$ in $\mathrm{\Omega }\backslash A_P$ and ${\lambda }_S$ free in $A_S$ , ${\lambda }_P$ free in $A_P$ ;  using starting values from last iteration;  forall elements $T$ do   if $T\in A_P,T\notin A_S$ and ${\lambda }_P<{\delta }_P$ then     remove $T$ from $A_P$   else if $T\notin A_P$ and $f_P>{\delta }_P$ then     add $T$ to $A_P$ ;   end   if $T\notin A_S,T\in A_P$ and $f_S>{\delta }_S$ then     add $T$ to $A_S$   else if $T\in A_S$ and ${\lambda }_S<{\delta }_S$ then     remove $T$ from $A_S$ ;   end  end end

Note that, in order to avoid numerical instabilities Kuhn-Tucker compatibility conditions in Algorithm 1 are not checked for zero. The switching and saturation and saturation condition as well as the corresponding Lagrangian multipliers are compared to a sufficiently small value ${\delta }_P$ and ${\delta }_S$ , respectively. These small offsets shall prevent recurring switching of elements between active and inactive state. In our computations the values have been chosen ${\delta }_P=E_C·{10}^{-3}$ and ${\delta }_S=P_{\mathrm{sat}}·{10}^{-3}$ . These values lead to a maximum of three loops per load step for the three dimensional examples presented in the next chapter.

5.1. Implementation of dissipation

In order to obtain an incremental formulation the dissipation has to be reformulated. The rate of irreversible polarization in the time interval $\tau =t-t_n$ is given by

$${\stackrel{·}{\mathbf{P}}}_{I,n}=({\mathbf{P}}_I-{\mathbf{P}}_{I,n})/\tau =\mathrm{\Delta }{\mathbf{P}}_I/\tau .$$

(72)

Note that, in equation (72) a constant rate of polarization is assumed for one time interval. Taking into account (equation (32)) the time integral in equation (28) is replaced by multiplication by $\tau$ . The dissipated energy $D^{\tau }$ in the time interval $\tau$ is given by

$$D^{\tau }={\int }_{\mathrm{\Omega }}\hat{\mathbf{E}}·({\mathbf{P}}_I-{\mathbf{P}}_{I,n})\mathrm{d}\mathrm{\Omega }={\int }_{\mathrm{\Omega }}\hat{\mathbf{E}}·\mathrm{\Delta }{\mathbf{P}}_I\mathrm{d}\mathrm{\Omega }.$$

(73)

Note, that in equation (73) the driving force $\hat{\mathbf{E}}$ is taken at time $t$ . This implies the backward Euler method for the driving force.

6. Numerical results

In this section we show numerical results for several examples. First we show the standard hysteresis curves for the fully coupled homogeneous (uni-axial) problem. We show that the formulation allows mechanical depolarization. Our second example is the repolarization of ferroelastic material, that is initially polarized at a certain angle compared to the direction of applied electric field. Results are compared to the experiments carried out by Huber and Fleck (2001). Next we show the mechanical depolarization of a fully polarized ferroelectric beam. This example was introduced as a ferroelastic benchmark in Zouari et al. (2011). Finally we show the electric polarization of a bimorph structure. We use the same material data for all examples (except non-proportional loading) taken from Zouari et al. (2011) listed in Table 1.

Table 1.

Material parameters.

Young’s modulus $Y$	${10}^4\mathrm{MPa}$
Poisson’s ratio $\nu$	$0.3$
Coercive field $E_C$	$1\mathrm{MV}/m$
Saturation strain $S_{\mathrm{sat}}$	$0.002$
Saturation polarization $P_{\mathrm{sat}}$	$0.3C/m^2$
Hardening parameter $c$	$2·{10}^6\mathrm{Vm}/C$
Piezoelectric coefficient $d_p$	$5.93·{10}^{-10}m/V$
Piezoelectric coefficient $d_n$	$-2.74·{10}^{-10}m/V$
Piezoelectric coefficient $d_t$	$7.41·{10}^{-10}m/V$
Permittivity $ϵ$	$1.5·{10}^{-8}C/\mathrm{Vm}$

6.1. One-dimensional loading

In our first example we show that our implementation is capable of electrical polarization as well as electrical and mechanical depolarization. We consider an initially unpolarized unit square ( $1\mathrm{mm}$ by $1\mathrm{mm}$ ) with electrodes at top and bottom. The material parameters are listed in Table 1. Note, that for the two dimensional calculations plain strain is assumed. First the electrical polarization is shown. An electric field with a strength, for which the material will be fully polarized, here of two times the coercive field strength $2E_C$ is applied. Then the electric field is than reduced to $-2E_C$ , and again increased up to $2E_C$ . The resulting electric hysteresis, as well as the corresponding (mechanical) butterfly hysteresis are shown in Figures 3 and 4, respectively. In order to show the convergence of the method, calculations are performed with different step sizes. Note that for the step size $2E_C$ , full polarization and depolarization is reached within one load step, respectively.

Figure 3.

Electrical hysteresis of ferroelectric material.

Figure 4.

Mechanical butterfly hysteresis of ferroelectric material.

Next we show the effect of mechanical depolarization. The specimen is electrically polarized (field strength $2E_C$ ), then the electric field is reduced to $0$ . A mechanical compressive stress (aligned in polarization direction) is applied to the polarized material and mechanical depolarization can be observed. The mechanical depolarization curve is shown in Figures 5 and 6.

Figure 5.

Electrical hysteresis—depolarization of ferroelectric material.

Figure 6.

Mechanical butterfly hysteresis—depolarization of ferroelectric material.

6.2. Non-proportional loading

In our second example we verify the capability of non-proportional loading. A rectangular specimen is initially fully poled in a certain direction. Electrodes are located at top and bottom, applying an electrical field in vertical direction. The direction of polarization, referring to the electric field, is denoted by the angle $\alpha$ . The left and right surface of the specimen are insulated and considered to be free of charges. All boundaries are free of mechanical stresses, only rigid body motion is restricted. A sketch of the example, including geometry information, is shown in Figure 7. The material data for for this example are taken form Kamlah (2001) (electrical parameters) and Semenov et al. (2010) (mechanical and coupling parameters). They are summarized in Table 2.

Figure 7.

Sketch of repoling of ferroelectric material.

Table 2.

Material parameters for non-proportional loading.

Young’s modulus $Y$	$6.1·{10}^4\mathrm{MPa}$
Poisson’s ratio $\nu$	$0.31$
Coercive field $E_C$	$0.8\mathrm{MV}/m$
Saturation strain $S_{\mathrm{sat}}$	$0.0046$
Saturation polarization $P_{\mathrm{sat}}$	$0.23C/m^2$
Hardening parameter $c$	$0.9·{10}^6\mathrm{Vm}/C$
Piezoelectric coefficient $d_p$	$5.93·{10}^{-10}m/V$
Piezoelectric coefficient $d_n$	$-2.74·{10}^{-10}m/V$
Piezoelectric coefficient $d_t$	$7.41·{10}^{-10}m/V$
Permittivity $ϵ$	$6.2·{10}^{-8}C/\mathrm{Vm}$

As a consequence of the boundary conditions, which imply $\mathbf{D}·\mathbf{n}=0$ at insulated boundaries, for arbitrary angle $\alpha$ , homogeneous polarization is an incompatible state. Here this issue is solved by increasing the polarization in small steps and then solving for zero voltage (both electrodes grounded). Values of dielectric displacement are taken in the middle of the specimen, where boundary conditions do not affect the results. A more involved discussion about the effect of boundary conditions can be found in Stark et al. (2016c).

An electric field of size $2E_C$ is applied though the electrodes. Due to the applied electric field the direction of polarization is changed. The electric field and the change of dielectric displacement are measured at the center of the specimen. As a result the change of dielectric displacement in vertical direction over the applied electric field is shown for five particular angles alpha, ( $0^{°}$ , ${45}^{°}$ , ${90}^{°}$ , ${135}^{°}$ , ${180}^{°}$ ) in Figure 8. The results are visually compared to those from the experiments in Huber and Fleck (2001). The experimental data is shown in black lines in the back in Figure 8, our results are colored. We find good correlation between the measurements and the calculations.

Figure 8.

Non-proportional loading—repoling of ferroelectric material.

6.3. Depolarization of a ferroelectric beam

Our third example is a benchmark for ferroelastic media taken from Zouari et al. (2011). It is a cantilevered beam, fully polarized in its longitudinal direction and loaded with a tip force at its end. A sketch of the problem can be found in Figure 9. The computation is performed for a 3D-setting with square cross section ( $2\mathrm{mm}$ depth of the specimen). Grounded electrodes ( $\varphi =0$ ) are located at the clamped and at the tip face (illustrated in red). The clamping is realized by restricting longitudinal displacement and rigid body motions at the clamped face (blue dashed line).

Figure 9.

Sketch of fully polarized cantilever beam with tip load.

The tip force is realized as stress boundary condition, applying $\mathit{\sigma }\mathbf{xy}=2{\mathrm{N}/\mathrm{mm}}^2$ at the tip face. Taking into account the dimensions, this load equals the tip force of $8N$ in Zouari et al. (2011). The applied load causes bending of the beam, which leads to a compressive stress at the top face and therefore mechanical depolarization can be observed. The resulting polarization is shown in Figure 10. The resulting stresses are shown in Figure 11. Note that the irreversible polarization takes one (constant) value in each element, while the stress is interpolated. Due to the unsymmetrical remanent straining an unsymmetrical distribution of the bending stresses is to be observed. In Figure 12 the bending stress ${\mathit{\sigma }}_{\mathrm{xx}}$ on top and bottom of the beam, as well as the irreversible polarization on the top are plotted over the length of the beam. Note that, in contrast to the reference, only ferroelectric but no ferroelastic effects are considered in the constitutive model. While in Zouari et al. (2011) a decrease of the irreversible polarization of $0.82P_{\mathrm{sat}}$ is reported at the clamped face, here the minimum irreversible polarization is $0.5485P_{\mathrm{sat}}$ . Furthermore, in the reference the the top surface of the beam is fully polarized at locations $x>6\mathrm{mm}$ , here the state of full polarization is present at positions $x>8.78\mathrm{mm}$ .

Figure 10.

Depolarization of the ferroelectric beam under tip load.

Figure 11.

Stress ${\mathit{\sigma }}_{\mathrm{xx}}$ of the ferroelectric beam under tip load.

Figure 12.

Stress at top and bottom, polarization at top for polarized cantilever beam with tip load.

6.4. Bimorph structure

In our last example we show the polarization of a bimorph structure. A sketch of the structure, including geometry parameters is shown in Figure 13. The bimorph consists of two layers of different thickness, both of the same ferroelectric material with material parameters according to Table 1. Only the upper layer is electrically active, as the lowest and the middle electrode are grounded. The specimen is cantilevered in the same way as in the previous example, such that elongation within the fixed plane is not restricted. Figure 14 shows the mesh used for the calculation. As the used TDNNS elements do not suffer from locking effects, a very anisotropic (spatial) discretization can be chosen. In thickness direction in total six mesh layers are used, with very thin layers near the electrodes. The upper layer is fully polarized by applying a sufficient high electric field. After polarization the poling field is removed, all electrodes are grounded. Due to remanent straining only in the upper layer, the structure is deformed (Figure 14). Note that the final configuration is not free of mechanical stress, which further leads to (rather high) electric fields. The resulting stress component ${\mathit{\sigma }}_{\mathrm{xx}}$ , as well as the electric field in thickness direction $E_z$ are shown in Figures 15 and 16, respectively. The structure is shown in a deformed configuration with displacement scaled by a factor of $20$ . In Figure 17 the deflection $u_z$ is evaluated along two horizontal lines, one of them the bimorph’s center line ( $y=0$ , $z=0$ ), one on the lateral surface ( $z=0$ , $y=\pm w/2$ ). The deflection at the end of the center line is $u_z(x=b,y=0,z=0)=0.5379\mathrm{mm}$ , the deflection of the corner points is $u_z(x=b,y=\pm w/2,z=0)=0.5882\mathrm{mm}$ .

Figure 13.

Sketch of bimorph structure.

Figure 14.

Mesh for calculation of bimorph structure.

Figure 15.

Stress component ${\mathit{\sigma }}_{\mathrm{xx}}$ after polarization bimorph structure.

Figure 16.

Electric field component $E_z$ after polarization bimorph structure.

Figure 17.

Deflection of the bimorph structure after polarization at the center line ( $y=0$ , $z=0$ ) and on the lateral surfaces ( $y=\pm w/2$ , $z=0$ ).

7. Conclusion

In this paper we have shown the implementation of (known) material laws in the framework of variational inequalities. This allows the direct implementation of inequalities in frameworks based on variational concepts as for example finite elements. As shown, inequality constraints can be included in the calculation process and are taken into account directly (no return mapping is used). A finite element discretization is chosen using normal-normal stress and tangential displacements for the mechanical degrees of freedom, and electric potential and remanent polarization for the electrical degrees of freedom as well as two (additional) Lagrangian multipliers for the inequality constraints. The choice of free variables allows for direct use of the piezoelectric tensor $\mathbf{d}$ .

The performance of the method was shown within several numerical examples, where results according to the literature were provided.