Variational inequalities for ferroelectric constitutive modeling

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.


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, Ghaleb and Maugin [Bassiouny et al., 1988a, Bassiouny et al., 1988b, Bassiouny and Maugin, 1989a, Bassiouny and Maugin, 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 and co-workers [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].
We will follow Kamlah's approach, but formulate the theory in the frame-work of variational inequalities. These inequalities arise naturally from energybased 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. 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 Schöberl, 2011, Pechstein andSchöberl, 2018]. These elements were adapted to linear piezoelectricity in [Pechstein et al., 2018, Meindlhumer andPechstein, 2018].
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 form the literature, including electrical polarization, mechanical depolarization, nonproportional loading and depolarization under bending stress. Additionally we show the polarization of a bimorph structure. Last, a short review is provided.

A phenomenological material law
In this section we briefly introduce the phenomenological material model in the spirit of [Kamlah, 2001], before derivation an energy-based approach in the next section. We consider a ferroelectric body in the domain Ω ∈ R 3 with the boundary ∂Ω. We are concerned with isothermic, rate independent and quasi-static deformation and polarization processes. We are interested in finding the displacement field u, the electric potential φ and the remanent polarization P I .
The electric field is introduced as the negative gradient of the electric potential E = −∇φ, whereas we use the linear strain tensor ε = 1 2 ∇u + ∇u T . The mechanical balance equation and Gauss' law are for the stress tensor σ and the dielectric displacement D. The mechanical boundary conditions are and the electrical boundary conditions with n denoting the normal vector on the corresponding boundary. The constitutive model relates stress σ and dielectric displacement D to strain ε and electric field E as well as to the remanent polarization 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 The mechanical compliance tensor S E as well as the dielectric tensor ǫ σ are assumed to be isotropic and constant. The components of the dielectric tensor depend on P I via with unit vector in direction of polarization and the Kronecker delta The unpolarized material does not show any piezoelectric coupling effects. The dielectric tensor grows (linearly) with the norm of 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 sate with polarization in direction of e 3 . As proposed by [McMeeking and Landis, 2002], we assume the remanent strain to be volumepreserving and depend only on the polarization via where S Sat characterizes the maximum possible remanent strain. The polarization strain ε I is a deviatoric uni-axial strain and does not cause any further remanent strain.
The evolution of the remanent polarization 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Ê which will be specified in the sequel of this contribution. In the simplest case, the absolute value ofÊ may not exceed E C , On the other hand, the saturation condition states that the remanent polarization is limited and must not exceed the saturation polarization P sat , which resembles the fully polarized state. The condition reads

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.

Reversible processes
The stored energy (or free energy) Ψ R in the domain Ω is given by with the (volumetric) energy density function ψ R . For coupled piezoelectric problems the energy density function ψ R reads Now, let us consider a finite time interval τ := t − t n , ε n and 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 τ ext , splits up into two parts first W V ext , the external work according to body forces, second W B ext , external work according to boundary forces and applied electric potential on boundaries We introduce, according to [Miehe et al., 2011], the potential (for reversible processes) Π τ R with its algorithmic representation The potential Π τ R contains the difference between stored energy and work done by external loads W τ ext in the considered time interval. As the path independence of work is an essential property of linear material behavior, W τ ext 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 The solution of (21) is unique if the potential Π τ R is convex, which is the case for stable piezoelectric materials. In the sequel of this contribution we will use, instead of strain ε and dielectric displacement D, mechanical stress σ and electric field 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 (15), via the enthalpy density h R by The enthalpy density h is related to the free energy via Legendre transformation, For the case of linear piezoelasticity the enthalpy density reads In analogy to (20) we introduce for a (finite) time interval τ = t − t n the potential Π H τ with its algorithmic expression where suffix n denotes quantities at time t n . Analogously to (21) the associated stress and electric field is given by the constitutive maximization Note, that the potential Π H τ in (25), contrary to Π R τ in (20) does not contain the external work. External forces and applied potential are rather enforced via the side constraints in (26). 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 φ.
The mechanical balance equation is typically treated by introducing the (mechanical) displacement u as corresponding Lagrangian multiplier. The enthalpy H is extended and reads The potential Π L τ corresponding to the extended enthalpy L is given by The constitutive maximization problem (26) transforms into a saddle point problem. It reads (29)

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 τ for the time interval τ = t n − t and the considered region Ω. It is given by with the local dissipation D . The second law of thermodynamics states that, for each time interval, D τ ≥ 0 has to be non-negative. This implies positive dissipation D ≥ 0 for arbitrary processes.
In this contribution we follow the approach of internal variables, and introduce irreversible polarization P I . Consequently the stored energy density function for dissipative processes ψ depends also on P I , Analogously to (15) the stored energy in the domain Ω is given by Note, that in this work, remanent straining is always assumed to depend directly on the remanent polarization as proposed by e.g. [McMeeking and Landis, 2002] or [Linnemann et al., 2009]. Therefore the irreversible part of the strain ε I is not to be considered as a free variable.
According to [Miehe et al., 2011] we introduce the driving forceÊ. It is defined as the negative derivative of the stored energy density ψ by the internal variables, in our case the irreversible polarization P Î The driving forceÊ is also referred to as as internal constitutive force, asÊ is the work conjugate of 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 The constitutive minimization principle (21) 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 Π τ D with its algorithmic expression During the time interval τ the energy stored is Ψ (ε, D, P I )−Ψ (ε n , D n , P I,n ) and D τ is dissipated.
Analogously to (21), 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 Note that, contrary to (21), the minimization prinziple (36) 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 (22), (23) but using the energy density ψ, With the Lagrangian L and its potential Π L τ defined in analogy to (27), (28) the following saddle point problem is obtained 4 Incremental formulation of ferroelectricity The focus of this section will be the realization of the energy density ψ and the reformulation of the switching constraint f P (Ê) ≤ 0. The energy density will be chosen such that it reflects the constitutive equations (5), (6) as well as the saturation condition (14). 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. [Landis, 2002, Kamlah, 2001, Tichỳ et al., 2010), we follow the approach of additive decomposition of free energy. It decomposes into two parts, with the reversible or stored part of the free energy ψ R and the additional contribution of the irreversible quantities to the free energy ψ I . The reversible part ψ R is given by As stated in (7) all material tensors may depend on the irreversible polarization P I . The additional contribution of the free energy is given by 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. (14)) is reformulated as a variational inequality, involving the Lagrangian multiplier λ S , To ensure the saturation condition, we add the following supremum term to the energy function: For states of irreversible polarization P I not violating the saturation condition (f S ≤ 0), the supremum in (43) takes zero value, and ψ is not altered in comparison to (39). Otherwise, if the |P I | exceeds the saturation polarization, the saturation condition is violated (f S > 0), and therefor sup λ S ≥0 λ S f S tends to infinity, and with it the energy function ψ. Note, that for all admissible values of P I -not violating the saturation conditionψ takes finite values, while for inadmissible (improper) values of P I , the free energy function ψ tends to infinity. In the spirit of optimization problems, this leads to solutions with admissible states of irreversible polarization P I .
For further deductions we treat the Lagrange multiplier λ S as an free variable. This is indicated by the index S, as we use In the final optimization problem λ S will be maximized.

Transformation to enthalpy
The enthalpy is a function of the free variables electric field E, mechanical stress σ and irreversible polarization P I and is related to the free energy via Legendre transformation (37). The enthalpy density h S including saturation condition (44) reads 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 P I . In contrast to the original model by Kamlah, the driving forceÊ contains saturation via the corresponding Lagrangian multiplier λ S . The driving force readŝ Note that (46) implies that via (34) the saturation is part of the dissipation D.

Incremental optimization principle
Again we consider the (finite) time interval τ = 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 (26) The algorithmic expression of the corresponding potential Π τ S reads The saddle point problem reads In order to consider the switching condition, another Lagrangian multiplier λ P is introduced in analogy to (42) via Consequently the Lagrangian L S and the corresponding potential Π S τ have to be updated by the terms in (50). The Lagrangian L P , involving polarization, reads The algebraic form of the corresponding potential Π P τ is given by the corresponding saddle point problem, extended by the (independent) variable λ P , reads

Variation of Optimization Problem
The saddle point problem (53) may be solved by the method of variation. As the Lagrangian multipliers λ S and λ P are constrained to be non-negative, the variation of the potential Π P τ leads to a variational inequality. Variation with respect to u, E, σ and P I can be done in the standard way. Arbitrary virtual values δu, δσ, δE and δP I are admissible, only the essential (Dirichlet) boundary conditions on the corresponding boundaries Γ 1 , Γ 2 , Γ 3 and Γ 4 have to be fulfilled, respectively. A variational equation is obtained via Note, that external (body) forces are already considered within the potential Π P τ . The saturation and polarization condition, which are given as inequalities, cannot be treated in the standard way. As the set of admissible λ P and λ S is restricted by inequalities δλ P and δλ 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, λ P and λ S a test functionλ P and λ S is introduced. The potential Π P τ is maximized with respect to λ P if and only if for all other admissible choicesλ P ≥ 0 there holds For the saturation and the corresponding Lagrangian multiplier λ S , the situation is a bit more involved, as the driving forceÊ depends on λ S . Writinḡ for some admissible test functionλ S ≥ 0, one finds that Π P τ is maximized with respect toλ S if and only if Inserting the definition ofÊ (46), which implies that the above inequality can further be reduced to Summing up, the variational inequality to hold for all admissible δE, δσ, δP I and δu as well as for all test functionsλ S ≥ 0 andλ P ≥ 0 reads (60)

Interpretation of the Lagrangian multipliers
In this section we explain the meaning of the Lagrangian multipliers λ P and λ 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 λ P and λ S , as well as the corresponding inequalities read the same as for the mechanical problem with stress states σ = 0.
For this electric case, as all mechanical quantities vanish, the driving force reduces to the simple formÊ The variational inequality (60) reads We consider an update step with the (updated) solutions E, P I , λ P and λ S . We introduce ∆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Ê · ∆P I . Splitting up (62) into individual equations and inequalities for the single variational quantities we get To get an interpretation of λ P , consider a state, at which the polarization process has started, such that f P = 0,Ê = 0 and λ 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 From the variation of the irreversible polarization δP I in (64), taking into account the definition (61) ofÊ and (67) we get a relation between λ P and ∆P I via From (68)  Before we give an interpretation for λ S , we first show that (66) really ensures the saturation condition. In case the coercive electric field is reached and f P = 0, (68) holds, and thereby the last term in (66) 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 λ P = 0 and ∆P I = 0, which again reduces the last term in (66) to zero. Thus the saturation condition has to be satisfied in any case.
For the interpretation of λ S we consider a state where f S = 0, i.e. saturation polarization is reached. The driving electric fieldÊ depending on λ S still has to satisfy the switching condition, The Lagrangian multiplier λ S grows with the electric field E as soon as further growth of the remanent polarization is prohibited by the saturation condition. A graphical interpretation of λ P and λ S can be found in Figure 1 and Figure

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 normal stresses as degrees of freedom, which motivate the abbreviation "TDNNS". This method was originally developed for elastic solids Schöberl, 2011, Pechstein andSchöberl, 2012] and later extended for linear piezoelectric materials [Pechstein et al., 2018, Meindlhumer andPechstein, 2018] and for geometrically nonlinear electro-mechanically coupled problems [Pechstein, 2019]. This work is based on [Pechstein et al., 2018, Meindlhumer andPechstein, 2018] for linear piezoelastics. Note that, for the mixed TDNNS method, the d-tensor formulation can be used directly, and there is no need to transfer the electric permittivity at constant stress ǫ σ to that at constant strain ǫ ε 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 Ω σ : ε or Ω div σ · 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 λ P and λ S are set to zero, respectively. For these fixed active sets, the nonlinear equations are solved by Newton's method, with λ S and λ 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 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.

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 τ = t − t n is given byṖ I,n = (P I − P I,n ) /τ = ∆P I /τ.
Note that, in (70) a constant rate of polarization is assumed for one time interval. Taking into account (34) the time integral in (30) is replaced by multiplication by τ . The dissipated energy D τ in the time interval τ is given by Note, that in (71) the driving forceÊ is taken at time t. This implies the backward Euler method for the driving force.
Data: values for u 0 , σ 0 , φ 0 , P I0 , λ S0 , λ P 0 and active sets A S0 , A P 0 from the last converged time step Result: values for u, σ, φ, P I , λ S , λ P , active sets A S , A P initialize unknows with initial data, active sets A S = A S0 , A P = A P 0 ; while the active sets change do solve non-linear problem with λ S = 0 in Ω\A S , λ P = 0 in Ω\A P and λ S free in A S , λ P free in A P ; using starting values from last iteration; for all elements T do if T ∈ A P , T / ∈ A S and λ P < 0 then remove T from A P ; else if T / ∈ A P and f P > 0 then add T to A P ; end if T / ∈ A S , T ∈ A P and f S > 0 then add T to A S ; else if T ∈ A S and λ S < 0 then remove T from A S ; end end end Algorithm 1: Active set strategy for solving the variational inequality.

NUMERICAL RESULTS
In this section we show numerical results for several examples. First we show the standard hysteresis curves for the fully coupled homogeneous (uniaxial) 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. 1 MV/m Saturation strain S sat 0.002 Saturation polarization P sat 0.3 C/m 2 Hardening parameter c 2 · 10 6 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/Vm

On-dimensional loading
In our first example we show, that our implementation is capable of both, electrical polarization as well as electrical and mechanical depolarization. We consider an initially unpolarized unit square (1 mm by 1 mm) with electrodes at top and bottom. The material parameters are listed in Table 1. 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 Figure 3 and Figure 4, respectively. 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 Figure 5 and Figure 6.

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 α. The left and right surface of the specimen are isolated 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. Hardening parameter c 0.9 · 10 6 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/Vm As a consequence of the boundary conditions, which imply D · n = 0 at isolated boundaries, for arbitrary angels α, homogeneous polarization, is an incompatible state. Here this issue is solved by applying 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., 2016].
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

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 mm depth of the specimen). Grounded electrodes (φ = 0) are located at the clamped and at the tip face (illustrated in red). The clamping is realized by restricting longitudinal displacement at the clamped face (blue dashed line), furthermore all rigid body motions are restricted via the left face.
The tip force is realized as stress boundary condition, applying σ xy = 2 N/mm 2 at the tip face. Taking into account the dimensions, this load equals the tip force of 8 N 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. It decreases at the clamped face to 55% of the saturation polarization. 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. The results fit to [Zouari et al., 2011]. Note that, in contrast to the reference, only ferroelectric but no ferroelastic effects are considered in the constitutive model.

Bimorph structure
In our last example we show the polarization a bimorph structure. A sketch of the structure, including geometry parameters is shown in Figure 12. 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 12 shows the mesh used for the calculation. As the used TDNNS elements do not suffer form locking effects, a very unisotropic (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. Note that the final configuration is not free of mechanical stress, which further leads to (rather high) electric field. The resulting stress component σ xx , as well as the electric field in thickness direction E z are shown in Figure 14 and Figure 15, respectively. The structure is shown in a deformed configuration with displacement scaled by a factor of 20.

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 e.g. 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 d.
The performance of the method was shown within several numerical examples, where results according to the literature were provided.