The polarization process of ferroelectric materials analyzed in the framework of variational inequalities

We are concerned with the mathematical modeling of the polarization process in ferroelectric media. We assume that this dissipative process is governed by two constitutive functions, which are the free energy function and the dissipation function. The dissipation function, which is closely connected to the dissipated energy, is usually non-differentiable. Thus, a minimization condition for the overall energy includes the subdifferential of the dissipation function. This condition can also be formulated by way of a variational inequality in the unknown fields strain, dielectric displacement, remanent polarization and remanent strain. We analyze the mathematical well-posedness of this problem. We provide an existence and uniqueness result for the time-discrete update equation. Under stronger assumptions, we can prove existence of a solution to the time-dependent variational inequality. To solve the discretized variational inequality, we use mixed finite elements, where mechanical displacement and dielectric displacement are unknowns, as well as polarization (and, if included in the model, remanent strain). It is then possible to satisfy Gauss' law of zero free charges exactly. We propose to regularize the dissipation function and solve for all unknowns at once in a single Newton iteration. We present numerical examples gained in the open source software package Netgen/NGSolve.


Introduction
A thermodynamical framework for the description of ferroelectric materials based on the Helmholtz free energy was originally provided in the series of papers [1,2,3,4] by Bassiouny, Ghaleb and Maugin. Their theory allows to describe multiaxial electromechanical loading procedures. The introduced notions are similar to elasto-plasticity, including internal variables, yield (or switching) criteria and hardening moduli. Explicit choices of energy and switching criteria were provided by Cocks and McMeeking [8] for the onedimensional case. The multi-dimensional case followed in the works of McMeeking and Landis [15] and Landis [14]. In the former reference, the remanent polarization vector is the only internal unknown, and the polarization stress is linked directly to the remanent polarization. Contrarily, in the latter work remanent polarization and strain are independent of each other, but are determined by a common switching condition. With these theories, not only hysteresis loops can be tracked, but also butterfly hystereses are predicted correctly. The models were validated against measurements provided by Huber and Fleck [10] for non-proportional loading procedures.
Another approach based on the thermodynamic framework due to the group around Maugin is that by Kamlah and Tsakmakis [12], see also [11]. They use a set of different switching and saturation conditions to determine the evolution of remanent polarization and polarization strain.
Miehe, Rosato and Kiefer [16] introduced an incremental variational principle for the even more general case of coupled electro-magneto-mechanics. They distinguish between energy-based and enthalpy-based models. For the latter, the independent unknowns are strain and electric field, whereas for the former, strain and dielectric displacement are independent. While most finite element formulations use enthalpy-based models discretizing the electric potential, we provide theory for an energy-based setting with an independent dielectric displacement.
Semenov et al. [17] present an energy-based formulation using a vector potential for the dielectic displacement. Then, the dielectric displacement vector satisfies Gauss' law of zero divergence automatically. Contrarily, we propose to use H(div) conforming finite elements as can be found in the context of mixed methods [5]. Gauss' law is then introduced as a constraint enforced exactly by a Lagrangian multiplier. This way accuracy of the electric unknowns is improved, as no derivatives have to be taken. In some sense, the approach by Klinkel [13] may be seen as diametrically opposed, where an irreversible electric field is introduced instead of the remanent polarization vector.
In the current paper, we aim at proving existence and uniqueness of a solution to the problem of finding an update solution in the polarization process of ferroelectric media. To this end, we describe the polarization of ferroelectric media as a dissipative process in the mathematical framework of variational inequalities. We start from an energy-based model of the problem, from which we derive time-dependent variational equations and inequalities. The independent unknowns are then strain (or displacement), dielectric displacement, remanent polarization, and, if included in the model, remanent polarization strain. We provide an abstract mathematical framework for this problem. Similar abstract problems have been analyzed in context of elasto-plasticity (see e.g. the monograph by Han and Reddy [9]) or contact mechanics (we refer to Sofonea and Matei [18]). We proceed in a similar manner as in the first reference [9] and see that, under the standard assumption of convexity of the free energy function, existence and uniqueness of a solution to the time-discrete update problem can be shown. Two standard material models are demonstrated to fit into this framework. These results fit well with the findings in [19,6], where stability issues are discussed. Under further assumptions on the energy, also existence of a time-dependent solution can be guaranteed. However, we see that these stronger assumptions are only satisfied for very simple models not including saturation. Whether it is possible to weaken these requirements will be subject of further research. This paper is organized as follows: in Section 2 the underlying energy-based consitutive models are described, and remanent quantities, dissipation function and dissipative driving forces are introduced on a physical level. These quantities are embedded into an abstract mathematical framework of variational inequalities in Section 3, where also all assumptions are stated. The time-discrete update equation is introduced in Section 4. Existence and uniqueness of a solution update are shown. In Section 5, this result is used to gain existence of a time-dependent solution to the original variational inequality under stronger assumptions. Two standard material models are analyzed with respect to the abstract theory in Section 6. Computational aspects of a finite element implementation are discussed in Section 7, and numerical results are provided in Section 8.

Energy-based constitutive modeling
In the following, we present a variational inequality that describes the problem of polarization of ferroelectric media. We postpone mathematical exactness concerning solution spaces, (weak) differentiability and other issues to Section 3, in order not to complicate matter too much. Let Ω ⊂ R d , d = 2, 3 denote the body of interest. Concerning the mechanical quantities, we use u for the displacement field and T for the total stress. We assume small deformations, thus we identify deformed and undeformed configuration, and use the linear strain S = 1 2 (∇u + ∇u T ). Additionally, we assume a quasistatic regime. The electric part of the problem is then characterized by the electric potential φ, while E = −∇φ is the electric field. The dielectric displacement vector shall be denoted as D.
Last, in the polarization problem we are interested in finding the remanent polarization P i and the remanent strain S i .
We assume that the material is characterized by a Helmholtz free energy function which consists of a reversible or stored part Ψ r , and a part Ψ i that is associated only with the internal variables of remanent polarization and strain. This part may tend to infinity as the polarization approaches saturation. For a similar characterization, we refer to [14]. The independent reversible unknowns are strain and dielectric displacement, where the latter additionally has to satisfy Gauss' law, div D = 0. ( Electric field and stress are dependent quantities, and defined as derivatives of the free energy with respect to dielectric displacement and strain, As strain and dielectric displacement are not constrained, or constrained to the linear subspace of divergence-free functions in the latter case, there holds the following variational equation, E · δD + T · δS = δW ext for all δD, div δD = 0 and δS = S(δu).
In (4), all virtual strains and divergence-free virtual dielectric displacements are considered, where δD and δu have to satisfy the respective boundary conditions. On the right hand side, all virtual work by external forces is summarized as δW ext . This formula is well-known as principle of virtual works. Polarization of ferroelectric materials is a dissipative process. The dissipative driving forces dual to the irreversible quantities arê The dissipation is then given by the inner product of these driving forces and the dissipative fluxesṖ i andṠ i , The dissipation function Φ relates the driving forcesÊ,T to the driving ratesṖ i ,Ṡ i . Typically, the dissipation function is not smooth, but weakly lower semicontinuous and allows for a subdifferential 1 . The driving forces are contained in the respective subdifferentials, i.e.Ê 1 The subdifferential of some function j with respect to u is denoted by ∂uj and represents a set, namely By definition of the subdifferential, eq. (8) is equivalent to the variational inequality To arrive at one single variational inequality, we add up (4) and (9). In (4), we use the admissible virtual dielectric displacement δD =D −Ḋ, where divD = divḊ = 0. The virtual strain is defined from the admissible displacement updateũ andu by δS = S(ũ) − S(u). With these choices, we deduce for allũ,D,P i andS i ,

A mathematical framework
We define the variational inequality that describes the polarization process in a mathematical framework. Therefore, we introduce compact notation, which is compatible with the literature on variational inequalities, especially with the monograph [9]. We use A priori, we assume the different quantities to live in the following Hilbert spaces, Compound spaces for the compact notation are then, Depending on the definition of the irreversible energy Ψ i , the free energy Ψ may tend to infinity as the material approaches polarization saturation. We introduce the effective domain of Ψ by (cf. [9, p 73]) We define the nonlinear operator A : X → H * by its action Note that, in physical meaning according to (3) and (5), the operator A maps reactions to forces, i.e.
The work of external forces δW ext shall be represented by the linear functional ℓ ∈ H * . Note that, for the present choice of independent unknowns u and D, the external work also contains boundary conditions on the electric potential. In the exemplary case of a body under mechanical volume load f , surface tractions t on the boundary part Γ trac , and a prescribed potential φ 0 on the electrodes Γ el , this functional is defined as Thus, the variational inequality (10) translates to the abstract variational inequality of the form: Find w : The variational inequality can be extended to hold for all z ∈ H, i.e. also for dielectric displacement updates with non-zero divergence. To this end, the dissipation function needs to be augmented by this restriction, see [9]. Indeed, the augmented dissipation maps all dielectric displacements with non-zero divergence to infinity. In the following, it shall be denoted by j : H → R + := R ∪ {+∞} and is defined as Note that, when restricted to H 0 , j(·) and Φ(·) are equivalent. Therefore, we will use j(·) in the following, and consider the variational inequality (equivalent to (22)),

Assumptions
We collect assumptions on the various functionals that we need to proof existence and uniqueness of the update in a time-discrete scheme. In Section 6, we show whether these assumptions are satisfied for some standard material models. The following definitions are taken from Brezis [7].
The functional φ is called lower semicontinuous if and only if for all u ∈ E, and all sequences u n → u in E, It is weakly lower semicontinuos if (26) holds for all weakly convergent sequences u n ⇀ u.
is continuous.
As a minimal assumption we demand the Helmholtz free energy to be lower semicontinuous and convex. Moreover, its Frechet (or full) derivative shall exist and define the operator A: Ψ : X → R is convex, lower semicontinuous and Frechet differentiable with A = DΨ (28) To get convergence estimates and stability bounds, we further need that A : X 0 → H * is strongly monotone, i.e. there exist m > 0 such that for all w 1 , w 2 ∈ X 0 , The dissipation function Φ : H → R as well as its augmented counterpart j : H → R + need not be differentiable, but non-negative, proper, and positively homogeneous, j(w) ≥ 0 for all w, z ∈ H, and it exists at least one w ∈ H with j(w) < ∞, and (30) Additionally, we assume j : H → R + is convex and lower semicontinuous, and (32) In Section 5, we aim at showing existence of a solution to the time-dependent variational inequality. To accomplish this task, we need an additional assumption on the energy, namely Lipschitz-continuity of A,

Time discrete update equation
We use a uniform partitioning of the time interval [0, T ] into N sub-intervals, For N fixed, we will define a sequence {w n } N n=0 ∈ [H 0 ] N +1 as consecutive solutions to (spatial but time-independent) variational inequalities. We use the backward difference ∆w n = w n − w n−1 and ℓ n = ℓ(t n ). We show that this sequence is defined uniquely, and that certain stability estimates are satisfied. We use the following existence result by Brezis: Theorem 1 (Corollaire 30 in [7]). Let E be a reflexible Banach space, and let X ⊂ E be closed and convex with 0 ∈ X. Let A : X → E * be weakly pseudo-monotone and φ : X →] − ∞, +∞] be convex lower semicontinuous with φ(0) < ∞. If then for ℓ ∈ E * there exists a solution u ∈ X to If A is additionally strongly monotone, the solution is unique.
The main result of this section is the following: Theorem 2. Let Ψ : X → R be convex, lower semicontinuous and Frechet differentiable as in (28), and let A : X → H * be its derivative. Let j : H → R + be non-negative, proper, convex, positively homogeneous and lower semicontinuous as in (30), (31) and (32). Let moreover A : X 0 → H * be strongly monotone on X 0 as in (29). Then, for N fixed and any given {ℓ n } N n=0 with ℓ n ∈ H * and ℓ 0 = 0, there exists a unique sequence {w n } N n=0 such that w n ∈ X 0 and ∆w n ∈ H 0 and The set of test functions can be equivalently restricted to z ∈ H 0 . With the constant of monotonicity m from (29), the solution satisfies the stability estimate Proof. We use Theorem 1 from convex analysis to show existence and uniqueness of the solutions. For the stability estimates, we progress along the lines of proof of [9] and see that some of their assumptions can be weakened. To show existence and uniqueness of the sequence {w n } N n=0 , we proceed inductively from w 0 = 0, assuming w n−1 to be known. We rewrite the variational inequality (38) in terms of the unknown ∆w n and w n−1 , We show that we can apply Theorem 1 to obtain existence of a solution ∆w n . As solution space we choose E = H 0 . We see that admissible updates are in the set X 0,wn−1 : Proposition 1] we know that convexity of Ψ implies monotonicity and hemicontinuity of A on thepossibly open -set X 0,wn−1 , which further implies pseudo-monotonicity of A. Obviously, strong monotonicity of A and positivity of j imply condition (36). But still we cannot use X 0,wn−1 directly for X in Theorem 1, as this set is not necessarily closed. Instead, we use the parameter-dependent closed sub-set with the choice of c > 0 still to be determined. We show that this set is closed: Assume a sequence z k ∈ X 0,wn−1 (c) that converges strongly to some z ∈ H. To show closedness of X 0,wn−1 (c) we need to prove that z ∈ X 0,wn−1 (c). Starting from the defining condition of (41) applied for z k , we apply lim inf k→∞ on both sides. Proposition 6 in [7] and continuity of ℓ n ensure that Thus we have shown z ∈ X 0,wn−1 (c) due to (41), and further that X 0,wn−1 (c) is closed for any fixed c > 0. The functional φ = j satisfies the conditions of Theorem 1. With all assumptions of Theorem 1 satisfied, there exists a solution ∆w c n to the parameter-dependent variational inequality As A is strongly monotone, the solution is unique . It remains to be shown that the variational inequality holds also for test functions z ∈ H 0 \X 0,wn−1 (c) provided c is larger than some fixed value not depending on the solution. To preserve uniqueness, we have to show that no w ∈ X 0 \X 0,wn (c) can be an additional solution.
For the first task, choose c n−1 = c(w n−1 ) depending on the previous iterate such that This constant exists since for w n−1 ∈ X 0 fixed, the derivative A(w n−1 ) is a continuous linear operator. Now set c = 2c n−1 , and let z ∈ X 0 \X 0,wn−1 (2c n−1 ), i.e.
From the hemicontinuity of A we deduce that there exists some t 0 > 0 such that for all t ∈ [0, t 0 ] i.e. tz ∈ X 0,wn−1 (2c n−1 ). But then the variational inequality (44) is satisfied for tz, and due to linearity and positive homogeneity of j, it is thus satisfied also for test function z ∈ H.
For the second task, we have to show that there cannot exist any ∆w ∈ X 0 \X 0,wn−1 (2c n−1 ) that is also solution to (40). But for such a function ∆w we know − A(w n + ∆w), ∆w + ℓ n , ∆w < −2c n−1 w < 0.
Together with the positivity of dissipation, one immediately obtains that the variational inequality (40) is not satisfied for z = 0, namely it holds Thus ∆w n is the only solution to (40). We proceed to the stability estimates (39). We additionally assume strong monotonicity of A as in (29). Note that in this case, condition (36) is trivially satisfied. In (40), we set z = 0 to obtain A(∆w n + w n−1 ), ∆w n + j(∆w n ) ≤ ℓ n−1 + ∆ℓ n , ∆w n .

Existence of a time-dependent solution
We proceed to finding a time-dependent solution from series of time-discrete solutions. A similar approach can be found in the framework of elasto-plasticity [9]. We note that this approach is intrinsically different from the procedure used by Sofonea and Matei [18] in contact mechanics. This approach [18] uses viscosity to ensure existence and uniqueness, which is not present in our problem.
We generate a time-dependent solution w N (t) from the series {w n } N n=0 from the previous section. We do so by linear interpolation in time, As, for any time t, w N (t) is a convex combination of w n ∈ X 0 , we observe that w N (t) ∈ X 0 . The following lemma proves that this interpolated solution satisfies a modified timedependent variational inequality for a certain kind of step functions z N .
Let A : X 0 → H * be additionally Lipschitz continuous as in (34). Then, there exists some constant c > 0 such that the following time-dependent inequality is satisfied, Proof. For each time step 1 ≤ n ≤ N , the variational inequality (38) is satisfied, if we choose test function z = ∆T /2(z n + z n+1 ) with z N +1 = 0, and use the positive homogeneity of j, Summing over n leads to In [9] the following estimates have been shown, Using the Lipschitz continuity of A, we see that the first part of the first sum in (60) is close to a corresponding integral. We estimate Using Cauchy-Schwarz inequality and estimate (54), we further derive (70) By similar conclusions we estimate the difference between the second part of the first sum in (60) and the corresponding integral, Putting these estimates together, and observing the positivity of j(z 1 ), we arrive at the desired inequality.
Next, we provide a candidate for the time-dependent solution of the original variational inequality. We observe that w N is bounded in the sense that, for any N , This follows directly from Lemma 1. Additionally, from the Lipschitz continuity of A we know that Now, consider a fixed step size N 0 and the according sequence of step sizes N l = 2 −l N 0 . To this sequence, there exists a weakly convergent subsequence, without loss of generality again denoted by N such that Applying lim sup and lim inf to the left and right hand side of the variational inequality from Lemma 2, we know that (with ∆T = T /N ), Obviously, the limit of the right hand side is zero, as ∆T = T /N → 0. In [9] it is shown that Due to the weak convergence of A(w N ) (77), we see Last, we observe from the chain rule of differentiation, the weak lower semicontinuity of the energy Ψ and the weak-star convergence of w N in L ∞ (0, T ; H) that lim sup Collecting these results, we find that w satisfies the following variational inequality Approximating z ∈ L(0, T ; H 0 ) by step functions and using a localization argument in time, the weak limit w can be shown to satisfy the time-dependent variational inequality (24). For details on this procedure, we refer to [9, p. 165]. Thereby, we arrive at the desired result: There exists a solution w ∈ H 1 (0, T ; H) solving the time-dependent variational inequality Proof. We have seen that w above satisfies the variational inequality, and that w ∈ L ∞ (0, T ; H) andẇ ∈ L 2 (0, T ; H). From the Sobolev embedding theorem we deduce that then w ∈ H 1 (0, T ; H).

Application to different ferroelectric material models
In the following, we motivate in how far the assumptions of the previous sections are reasonable, and if they hold for some standard material models. We do so first for a simple model without saturation, and for the ferroelectric material model for non-remanent straining as proposed by Landis [14]. We derive the form of the dissipation function Φ(Ṗ i ) for a given switching surface depending onÊ. We will see that the assumptions for existence of an update solution (Theorem 2) are met in both cases if material parameters are in common ranges, but that Lipschitz continuity lacks for the latter model, and thereby convergence in time is not guaranteed by our deductions. In all cases, we assume the reversible part of the energy as proposed by [14], In the theoretical considerations below, we assume the case of non-remanent straining, i.e. S i = 0 and c, h and β depend on the remanent polarization. Quite similarly, one might introduce a kinematic assumption for the remanent strain S i = S i (P i ), as done in [15,16]. We consider this case in our numerical examples. The different material models will be characterized via different irreversible energies. We assume that the material constants are such that the compound material tensor has positive eigenvalues bounded away from zero independently of P i . This issue is treated in detail by Stark et al. [19], they give suitable conditions. However, they claim that these conditions are not met by all data-sets provided for commercially available ferroelectric materials. Considering the special form of Ψ r , in case these conditions are met, the derivative A r := ∂Ψ r /∂w is not strictly monotone, but satisfies Before we turn to the different material models, we first characterize the dissipation function. Usually, it is defined by a threshold or switching surface, i.e. a condition on the forcesÊ andT. We assume a condition of the form Then the dissipation function Φ is given by In [9] it is shown that dissipation functions of the above format satisfy all assumptions such as convexity, lower semicontinuity and positive homogeneity. Indeed, in case of elasto-plasticity the dissipation function is of the same abstract form, where the coercive electric field E 0 resembles the yield stress.

A simple ferroelectric material model
A very simple material model in the spirit of Kamlah [11] assumes a quadratic dependence of the irreversible energy on the remanent polarization P i via Ψ i = Ω H 0 P i · P i dx. Saturation is not included in this model. In this case, the solution space is the whole space, X = H. As a quadratic form, obviously Ψ i is convex and Lipschitz continuous on H, and also in the notation of the previous sections.
We can now prove strict monotonicity of A = A r + A i . One immediately deduces by (89) and (92) that The last, essential ingredient to strict monotonicity is the fact that we restricted the full space of all dielectric displacements D to those which are divergence-free, div D = 0, and that this space D 0 is closed. In this case, the L 2 norm is equivalent to the full H(div) norm on D 0 , and (93) is sufficient for strict monotonicity of A on H 0 .

A saturating ferroelectric material model
For the second material model, the irreversible part of the energy is assumed as [14, Section 5] In [14] it has been shown that the free energy is differentiable, and all differentiations are provided analytically. The irreversible part A i := ∂Ψ i /∂w is strictly monotone in the remanent polarization, Thus, strict monotonicity of A = A r +A i follows the same way as in Section 6.1. However, A i is not Lipschitz continuous, as Although we can prove solvability of the time-discrete update equation, existence of a time-dependent solution is not guaranteed by our line of proof. However, in numerical examples, we did not meet any convergence problems.
7 Finite element implementation

Finite element spaces
We propose to use conforming finite element spaces for the discretization of the variational inequality. For a simplicial finite element mesh T = {T } and k ≥ 0, we use the nodal space of order k + 1 for the displacements, the Brezzi-Douglas-Marini space BDM k+1 for the divergence-conforming dielectric displacement (see e.g. [5]), and piecewise defined remanent polarizations of order k, Gauss' law of divergence free dielectric displacements is enforced by a Lagrangian multiplier in the sense of (23), which coincides with the electric potential φ and is discretized also by piecewise order k functions, Note that it is essential to choose dielectric displacement and its Lagrangian multiplier in a stable combination of spaces, such that not only but also the discrete inf-sup condition holds independently of the mesh size, This condition is satisfied for the pair BDM k+1 and piecewise order k functions. For a thorough theoretical background we refer the interested reader to the exhaustive monograph [5] on mixed problems. We mention that other stable choices of finite element pairs exist. In our numerical results, we used a subspace of divergence-free BDM k+1 elements. Then the electric potential is discretized using only one degree of freedom per element, regardless of the approximation order of D. In this case, it is impossible to evaluate the electric field as a derivative of ϕ. In any case, we recommend to use the constitutive law (3), as this leads to more accurate results.

Regularization of the dissipation function
For solving variational inequalities, various numerical algorithms are proposed in the literature. Well-known for dual variational inequalities is the return-mapping algorithm in different variants. There, after a predictor step, the generalized stress is projected back to the admissible set in the corrector step. Also for primal variational inequalities, as derived in this work, predictor/corrector iterations have been analyzed e.g. in the application of elasto-plasticity [9,Section 12.2]. In all these methods, after an "reversible", i.e. linear, predictor step the remanent quantities are altered in a consistent way in the corrector step.
In contrast, we propose to regularize the non-differentiable dissipation function, such that the problem can be solved "all at once" in a single Newton iteration. This regularization technique has been analyzed for convergence and accuracy in [9,Section 12.4]. Briefly, for a given regularization parameter ε, the non-differentiable dissipation j is replaced by a smooth function j ε , which differs from j only by ε (see Section 8 for a special choice). If the findings from [9] can be transferred to the ferroelectric polarization problem, one can expect that for the regularization parameter ε, • the solution w ε converges to w in H and However, we do not aim at proving these convergence estimates for the present problem.

Implementation in Netgen/NGSolve
We use the software package Netgen/NGSolve available at https://ngsolve.org. Netgen/NGSolve is an all-purpose finite element code, where high-order hierarchical finite elements for all element types (segments, triangles, quadrilaterals, tetraderda, hexahedra, prisms, . . . ) and many different spaces (continuous or discontinuous, curl or div conforming, . . . ) are implemented. Via a Python interface, variational equations or even energy formulations can be entered symbolically. The (symbolic) equations are differentiated automatically, and a Newton iteration can be realized in a straightforward manner, without need to implement tedious tangent stiffnesses etc. by hand for each formulation.
In the present manuscript, an energy formulation was used, where the free energy was entered analytically, and a regularized version of the dissipation was added,

Numerical results
We provide a patch test example, where we reproduce known hysteresis effects and mechanical depolarization on a ferroelectric cube. In the second example, a ferroelectric cantilever is polarized by an applied electric field, and partially depolarized in bending.
In both examples, we use the energies Ψ r and Ψ i as described in (88), (94) and (95). We assumed the permittivity at constant strain β = β S = ǫ −1 I to be independent of P i . The stiffness at constant electric field c E shall also be isotropic and independent of P i , and is characterized by Young's modulus E Y and Poisson ratio ν. The piezoelectric tensor d depends on P i in the standard way as given in [14, eq. (4.3)]. Then the coupling tensor h and the stiffness at constant dielectric displacement can be computed algebraically from β S , d and c E , compare [14]. We include remanent straining, where S i is assumed to depend directly on the polarization. We use the following formula provided in [15],

Patch test
Consider a cube of side length 2 mm whose normal displacement is fixed at the three coordinate planes. The cube is electroded on top and bottom, the other faces are electrically insulated. We use material constants derived from the dimensionless constants proposed by [15]. We set E 0 = 1000 V mm −1 , P 0 = 0.3 C m −2 , S 0 = 0.002, m = 2, The regularization parameter from Section 7.3 is set to ε = P 0 × 10 −6 .
We provide hysteresis curves for the standard load case of electric polarization and depolarization by an electric field of 1.5E 0 in Figure 1. Moreover, Figure 2 shows the effect of mechanical depolarization by a compressive load of 200 N mm −2 applied to the top surface of the cube.

Ferroelectric cantilever
The second example is that of a ferroelectric cantilever beam of length 2 mm and cross section 2 × 2 mm, which was proposed by [20]. The clamped end as well as the tip of the beam are electroded. In the first loading cycle, the beam is polarized applying an electric potential to the beam tip while keeping the other electrode grounded. The electric field is applied in 12 load steps amounting to three times the coercive electric field, then lowered back to ground. In the second step, a vertical tip force of 16 N is applied to the tip surface. Due to the compression in the upper part of the beam, the material depolarizes mechanically in this section.
The obtained values cannot be compared directly to the original work of [20], as in this reference a different material model based on Kamlah's work [11] is used. However, we chose material constants close to their values, using E 0 = 1000 V mm −1 , P 0 = 0.3 C m −2 ,  is set to ε = P 0 × 10 −6 . We used two different unstructured tetrahedral meshes -a coarse one consisting of 97 elements, and a fine one consisting of 6208 elements. For the reference solution, we chose order k = 1 as described in Section 7.1 on the fine mesh. This means second order displacement elements and first order polarization/dielectric displacements and leads to a total of 151786 degrees of freedom. The absolute value of the remanent polarization |P i | after bending, and the corresponding stress distribution T xx , are depicted in Figure 3. As observed in [20], the cantilever depolarizes in the region close to the clamped end where compressive stresses arise. This depolarization reduces the stress level there to a maximum of −128.22 N mm −2 , which compares well to the values listed in the original reference. However, due to the different description of the ferroelectric material, the region of depolarization is much larger. Indeed, the material already depolarizes when going back to zero voltage on the electrodes, which is due to the curved form of the hysteresis (see Figure 1). We compare these results to a second computation on the very coarse mesh using high order k = 2, i.e. third order displacement elements and second order polarizations/dielectric displacements. In this case, we end up with 6786 degrees of freedom in total, while maintaining good accuracy (see Figure 4). Note that, in both cases, neither polarization nor stresses have been post-processed in any way, but the finite element solution is depicted directly.

Conclusion
In this contribution, we have formulated the polarization problem in ferroelectric media as a variational inequality. In this framework, we were able to show existence and uniqueness of a solution to the time-discrete update problem under reasonable assumptions on the free energy. Under stronger assumptions, it is possible to prove existence of a solution to the time-dependent problem. We propose to choose finite elements such that these assumptions are satisfied also in the discretized setting. To solve the discrete problems, we regularize the non-differentiable dissipation function. Then is is possible to solve the optimization problem all at once by a single Newton iteration. All numerical results provided in this contributions have been generated in the open-source software package Figure 4: Mechanical depolarization of a polarized cantilever beam under a vertical tip forceabsolute value of remanent polarization |P i | (left) and stress distribution T xx . Finite elements as described in Section 7.1 for order k = 2 are used.
Netgen/NGSolve, which provides all the non-standard elements of arbitrary order as well as automatic differentiation of the energies.