High-order mixed finite elements for an energy-based model of the polarization process in ferroelectric materials

An energy-based model of the ferroelectric polarization process is presented in the current contribution. In an energy-based setting, dielectric displacement and strain (or displacement) are the primary independent unknowns. As an internal variable, the remanent polarization vector is chosen. The model is then governed by two constitutive functions: the free energy function and the dissipation function. Choices for both functions are given. As the dissipation function for rate-independent response is non-differentiable, it is proposed to regularize the problem. Then, a variational equation can be posed, which is subsequently discretized using conforming finite elements for each quantity. We point out which kind of continuity is needed for each field (displacement, dielectric displacement and remanent polarization) is necessary to obtain a conforming method, and provide corresponding finite elements. The elements are chosen such that Gauss' law of zero charges is satisfied exactly. The discretized variational equations are solved for all unknowns at once in a single Newton iteration. We present numerical examples gained in the open source software package Netgen/NGSolve.


Introduction
Piezoceramics are widely used as so-called smart materials for high-precision actuation and sensing. Due to the piezoelectric effect, mechanical strains give rise to dielectric displacements, whereas applied voltages introduce stresses. To quantify the piezoelectric effect, the polarization state of the smart material has to be known. In most applications, the polarization state is assumed as constant unidirectional. This initial polarization is obtained by applying a high electric field to the ferroelectric specimen. Polarization aligned with the electric field builds up, and a remanent polarization is maintained as this poling field is removed.
However, the remanent polarization state can change under sufficiently high electrical or mechanical loadings. Also, in applications such as some macro-1 fibre composites (MFCs), the poling electric field is not unidirectional, which will lead to non-standard polarization patterns. Therefore, the modeling of the polarization process in ferroelectric materials and the subsequent numerical simulation are an active field of research.
More precisely, macroscopic modeling of the ferroelectric behavior has been pursued for a long time, as a macroscopic mathematical model is prerequisite to any simulation using finite elements or other concepts. The model is required to catch the nonlinear material behavior especially under high electric fields and also depolarization under high pressures. Two inherently different approaches have been pursued successfully so far: microscopic-based and phenomenological models.
Microscopic-based models work on the level of microscopic cells or crystals, from which a macroscopic behaviour is derived by averaging over a large number of unit cells. This approach has been discussed e.g. by Huber et al. [20]. Finite element methods based on microscopic models date back to Chen and Lynch [10]. With increasing computing power, methods using one finite element per grain have become feasible, we cite non-exhaustively [1,29,23,33].
Another approach, which is adopted in this work, relies on a phenomenological description of the macroscopic material behavior. Internal variables then resemble macroscopic quantities, such as the remanent polarization or remanent polarization strain. We differ between thermodynamically consistent and non-consistent models.
Models based on a Preisach operator are not thermodynamically consistent. Parameters are chosen to emulate hysteretic curves. This approach has been chosen e.g. by [22,21,17]. Non-linear finite elements for the latter model have been proposed by [24]. A thermodynamically consistent model was first presented in the series of papers by the group around Maugin [2,3,4,5]. Kamlah [25] developed a theory for multi-axial loading based on several switching criteria defining the polarization and strain hystereses. A finite element implementation was proposed by Kamlah and Böhle [26]. More recently, Elhadrouz et al. [14] proposed a model including two different switching criteria for ferroelectric and ferroelastic switching.
Another family of thermodynamically consistent models has been developed, where ferroelectric and ferroelastic switching are combined in only one switching function. McMeeking and Landis [34] as well as Schröder and Romanowski [38] proposed formulations where the remanent polarization strain depends one-toone on the remanent polarization. Nevertheless, these models are capable of mechanic depolarization, see the numerical examples from [34]. Linnemann et al. [32] added magnetostriction in their work. Another approach with an independent polarization strain has been proposed by Landis [28]. Also Klinkel [27] uses two independent remanent quantities remanent electric field and remanent strain in his model. Laxman et al. [30] additionally modeled creep behavior in their approach.
Miehe et al. [35] presented a variational framework for electro-magnetomechanics. They introduced a concept of minimizing the work increment. They differ between approaches based on the free energy and electric enthalpy based approaches. Depending on this base, the free unknowns are either strain and dielectric displacement, or strain and electric field. While all of the models mentioned above are electric enthalpy based, only few energy-based models have been used so far. Sands and Guz [37] describe a one-dimensional model based on the free energy. Semenov et al. [39] provide a multi-directional energy based model using a vector potential for the dielectric displacement. They propose to discretize the vector potential in the finite element method, and provide consistent tangent moduli for the implementation. They observed very fast convergence of the Newton-Raphson iteration. Stark et al. [42,43] combine a phenomenological and micro-mechanically motivated model in their work.
The authors of the present work have analyzed such an energy-based formulation in the framework of variational inequalities in [36]. This concept has been used in the mathematical analysis of other nonlinear problems in mechanics, such as elasto-plasticity by Han and Reddy [16] or contact mechanics by Sofonea and Matei [41]. In the present work, a finite element method for the energy-based formulation is proposed. The utilized finite elements for displacement, dielectric displacement and remanent polarization are described in detail. These elements are chosen such that Gauss' law of zero charges inside the domain is satisfied exactly, what is usually not the case in electric enthalpy based methods. Numerical results for the problem of non-proportional loading and the homogenization of the polarization of a macro-fibre composite are presented.
The paper is organized as follows: First, the energy-based model in analogy to [35] is presented. Algorithmic representations for energy and dissipation are provided subsequently. A mathematical model of the time-discrete polarization problem is provided, which is non-differentiable due to the choice of dissipation function. Its regularization is addressed, and also the stability of the resulting variational equation. The next section is devoted to the finite element discretization, the choice of finite elements and the implementation in the open-source software package Netgen/NGSolve. The last section comprises numerical results.

Modeling dissipative material behavior
In the following, we assume to be in the standard setting of small-deformation mechanics in three space dimensions. The body of interest shall be denoted by Ω, which is in total or in part made of ferroelectric material. We use u : Ω × [0, T ] → R 3 to denote the displacement field, and use S = 1 2 (∇u + (∇u) T ) for the linearized strain tensor.
In an electrostatic regime on a simply connected domain, the electric field E can be shown to be a gradient field due to Faraday's law. We use standard nomenclature from the literature, introducing the electric potential ϕ with E = −∇ϕ. The dielectric displacement vector D satisfies Gauss' law of zero charge in the material, div D = 0 in Ω.
Last, we use the remanent polarization P i as an internal variable. The reversible material response is characterized by the internal energy Ψ or its density ψ, which shall depend on strain, dielectric displacement as well as remanent polarization vector, Total stress T and electric field E are conjugate to strain S and dielectric displacement D in the following sense, Conjugate to the remanent polarization is the (negative) electric driving forcê We denote by P ext = d dt W ext the power added from external sources such as applied forces or electric fields. As in mixed methods the external work often differs from standard approaches, we provide an exemplary expression below. For a given volume load f in Ω, a surface load t on Γ t and an applied electric potential V 0 on Γ el , the power from external sources is given by A reversible process is fully characterized by internal energy Ψ and external work W ext , and there are no internal variables such as the remanent polarization present. In this contribution however, we are concerned with irreversible, dissipative processes. As stated by Miehe et al. [35], the dissipative behavior can be described by the dissipation function Φ, which depends only on the rate of the internal variable. In our case, the dissipation function is thus a function of the remanent polarization rate, Depending on the dissipation function, rate-dependent or rate-independent material behavior is modeled. In the next section, we motivate our special choice of dissipation function, which leads to rate-independent behavior.
Using the notions introduced above, we cite the following incremental minimization problem for any time interval (t 0 , t 1 ), For a more detailed introduction, we refer to the extensive derivations originally provided by [35].

Choice of internal energy
Different choices of the internal energy density have been proposed in the literature. In many contributions, the internal energy is decomposed into two additive terms: the first of these terms coincides with the energy in linear problems for fixed remanent polarization. It is quadratic in the reversible parts of dielectric displacement and strain, but material tensors c D , h and β S may depend on the remanent polarization. The second term depends then only on the remanent polarization, and contains the electric hystersis shape parameters h 0 and m and the saturation polarization P 0 .
ψ =ψ r +ψ i , We use a one-to-one relationship between remanent polarization and polarization strain as suggested by [34], Note that this choice is volume preserving.

Choice of dissipation function
Apart from the free energy densityψ, the dissipation functionφ has to be specified to obtain a valid model through the minimization problem (9). The dissipation function is closely related to the dissipation D. The dissipation describes the amount of energy dissipated into heat and is always non-negative, as stated by the Clausius-Duhem inequality The choice of dissipation functionφ is motivated by the principle of maximum dissipation. It is closely related to the switching condition, which states in its simplest form that the remanent polarization rate is non-trivial only if the driving electric fieldÊ exceeds the coercive field E 0 , The dissipation function is chosen such that dissipation is maximized over all possible values forÊ,φ This supremum is of course taken forÊ = E 0 e P , with e P = P i /|P i | the unit vector in direction of remanent polarization, and forÊ = 0 in case |Ṗ i | = 0. Thus, the dissipation function can be simplified to read Note that this dissipation function is positively homogeneous (compare (20)). Such a dissipation function leads to rate-independent behavior, see [35].

Mathematical modeling
To get a well-defined problem in a mathematical sense, we introduce Hilbert spaces for the various quantities. Then, the time-discrete minimization problem (9) can be solved minimizing in the respective closed spaces. We use the standard H 1 Sobolev space of weakly differentiable functions with derivative in L 2 for the displacements. (Homogeneous) displacement boundary conditions on the boundary part Γ f ix are included, such that we are searching for The choice of the dielectric displacement space is less standard. We resort to the H(div) space of L 2 vector fields with weak divergence. Gauss' law (1) is then included explicitely. Additionally, insulation boundary conditions of vanishing surface charges D · n = 0 are posed explicitely. This means that these surface charge boundary conditions are essential boundary conditions, while we have already seen in the power statement (6) that voltage boundary conditions on electrodes are natural. Summing up, we search for 6 The remanent polarization does not sport any smoothness, nor does it allow for boundary conditions. Therefore we choose the vector-valued Lebesgue space We will now formulate the incremental variational problem for some small but finite time increment ∆t, starting from time t 0 to t 1 = t 0 + ∆t. To this end, we assume u 0 = u(t 0 ), D 0 = D(t 0 ) and P i 0 = P i (0) known. We are interested in finding the increments ∆u, ∆D and ∆P i such that u(t 1 ) = u 0 +∆u, D(t 1 ) = D 0 +∆D and P i (t 1 ) = P i 0 +∆P i . Miehe et al. [35] stated that these increments minimize the potential increment, In the special case of a rate-independent response, where the dissipation function is positively homogeneous, the time step size cancels out in (19). The parameter t acts as pseudo-time then, the minimization problem reads

Variational inequality
As the dissipation function is not differentiable atṖ i = 0, this minimization problem cannot be transformed into a variational equation. However, it can be treated in the framework of variational inequalities. Inspired by the mathematical analysis of elasto-plasticity (see e.g. [16]) or contact problems (see e.g. [41]), a thorough analysis of the arising time-dependent variational inequality has been carried out in [36]. In this contribution, we state the variational inequality which is equivalent to the optimization problem (21). We use the definitions for total stress T and electric field E (4) and for the electric driving forceÊ (5) evaluated at t 1 , Then, ∆u ∈ V, ∆D ∈ D 0 and ∆P i ∈ P satisfy the following variational inequality

Regularization
In the framework of elasto-plasticity, Han and Reddy [16] propose to regularize the dissipation function. The dissipation function is changed slightly, such that it becomes differentiable everywhere. When using the regularized dissipation function, the variational inequality (25) turns into a nonlinear variational equation, which can be solved by standard nonlinear solution algorithms. We propose to use the regularized dissipation function This choice has been analyzed to work well for elasto-plasticity in [16]. Indeed, according to this reference, we expect the error due to regularization to be of order √ ε. For completeness we provide the variational equation to be solved in each time step. We denote the (continuous) derivative of the regularized dissipation functionφ ε with respect to the remanent polarization rate asφ ε . The updates ∆u ∈ V, ∆D ∈ D 0 and ∆P i ∈ P satisfy Another issue for regularization is the free energy function ψ, or its additive part ψ i . The formula proposed in (10) is unbounded as the saturation polarization is approached, This unbounded potential acts as a penalty term as the polarization approaches saturation, where also (ψ i ) (P i ) diverges to infinity. As analyzed in [36], this fact does not harm solvability of the update equation (27) itself, but it may harm convergence. Therefore, if convergence problems are observed, we propose to modify the energy densityψ i in a way that (ψ i ε ) (P i ) becomes large but stays bounded. We use the following regularized expression (ψ i ε ) ,

Gauss' law
Gauss' law of zero charges (1) implies that the dielectric displacement field is divergence free inside the domain Ω. So far, this condition has been included in the space of admissible updates D 0 without further ceremony. As Gauss' law is a linear and continuous restriction, the space D 0 is a linear and continuous subspace of D and H(div). However, as this condition is non-local, it will be difficult to implement a finite element space consisting of divergence-free finite element functions in practice. In the following paragraph, we discuss how to reformulate the condition in an equivalent way, such that finit element implementation will be possible. We stress that, both in theory in this section, as well as in the finite element implementation, the dielectric displacement field will be exactly divergence free in the strong sense.
We propose to add a Lagrangian multiplier ensuring the divergence-free condition. As it turns out that this multiplier resembles the electric potential, we denominate it as ϕ. In the minimization problem, we propose to use an augmented energy Ψ 0 , which assumes infinity if Gauss' law is violated. This augmented energy reads Note that we use ϕ ∈ Q := L 2 (Ω), which means that discontinuous electric potentials are admissible testing functions. Also, the electric potential does not satisfy any boundary conditions a priori. Of course, as the approach is consistent, the potential will assume boundary values matching the prescribed external work (6). Also, it is not necessary to implement an equi-potential condition on electrodes specifically. The regularized variational equation including Gauss' law is stated below. The updates ∆u ∈ V, ∆D ∈ D and ∆P i ∈ P as well as the (total) electric potential ϕ ∈ Q satisfy for all δu ∈ V, δD ∈ D, δP i ∈ P, δϕ ∈ Q. The augmented variational equation (32) is equivalent to (27). The dielectric displacement update ∆D is exactly divergence free and thus ∆D ∈ D 0 .

Solvability and stability
The issue of existence of solutions, also to the time-dependent variational inequality (25), has been treated by the authors in [36]. We recall that a unique solution to the update problem (25) exists, given the following essential assumption of stong monotonicity holds: for all u,ũ ∈ V, D,D ∈ D 0 and P i ,P i ∈ P as well as corresponding forces T,T, E,Ẽ andÊ,Ẽ there holds This is the case if the reversible part of energy Ψ r is strongly monotone with respect to the reversible parts S − S i (P i ) and D − P i , and the additive part Ψ i is strongly monotone with respect to P i . This is the case, as long as the linear moduli lead to a positive definite material tensor. This issue is addressed by [44]. Conditions ensuring positive definiteness of the additive energy partψ i can be found in [7].

Finite elements
In this section we propose a set of consistent finite elements to discretize the variational equation (32) above. The elements shall be of arbitrary order, by k ≥ 1 we denote the lowest overall polynomial order. In the following, we assume that T is a simplicial mesh of the volume Ω. The elements can be defined also on prismatic meshes, as is shown in our numerical results. Note that, in the context of our mixed elements, "consistent" means different notions of continuity for the different fields.

Displacement elements
The displacement element has to be consistent for H 1 , such that u ∈ V. This condition is satisfied when using standard nodal or hierarchical continuous elements. We propose to use hierarchical elements of arbitrary order. For simplicial elements, the finite element space V h is then given by For tetrahedral elements, this results in (k + 2)(k + 3)(k + 6)/2 degrees of freedom, which means 30 degrees of freedom for k = 1 (second order displacment elements) and 60 degrees of freedom for k = 2 (third order displacment element). A prismatic element sports 3(k + 1) 2 (k + 2)/2 degrees of freedom. For k = 1 this results in 54 degrees of freedom (second order displacment elements), for k = 2 we obtain 120 degrees of freedom (third order displacment elements).

Dielectric displacement and electric potential
The dielectric displacement elements are less standard. To be consistent, we need that D ∈ D, i.e. that D allows for a divergence in weak sense. This means that the elements need to be normal continuous. In other words, different moments of the normal component D · n on element interfaces are degrees of freedom instead of nodal values. Different families of such normal-continuous elements have been introduced in the literature. For a fundamental introduction and overview, we refer the interested reader to the monograph [6]. As a first proposition, we propose to use elements from the BDM family as introduced by [8] (BDM after Brezzi, Douglas and Marini, who first proposed corresponding elements in two dimensions, see [9]). These elements span the space However, as the dielectric displacement is divergence free, the number of unknowns can be further reduced. As already mentioned, the divergence-free condition is non-local, therefore it is not possible to restrict the space D h such that it becomes a subspace of D 0 in a simple, local manner. However, in context of incompressible Navier-Stokes equation, [31] have shown that it is possible to do so for the higher-order shape functions. We propose to use this reduced set, div D| T ∈ P 0 for T ∈ T , D · n cont., D · n = 0 on Γ ins } ⊂ D.
We see that this way not only the number of unknowns for D drops, but also less Lagrangian multipliers ϕ are needed. As div D is constant per element, it is enough to use piecewise constant ϕ, Then, the condition implies that div D = 0 in strong sense for D ∈ D h0 . The finite element space D h0 is intrinsically different from standard nodal finite element spaces well-known in mechanics. While in nodal finite element spaces such as V h degrees of freedom are function values evaluated at element nodes, a more abstract concept of degrees of freedom as introduced by [11] is needed in the context of H(div). In our case, degrees of freedom are associated to element interfaces (to D · n), and additional "element bubbles" that are associated to the element interior (similar to element-interior nodes). There are no degrees of freedom associated to nodes or element edges. Different from standard nodal elements as used for the displacements, the inter-element coupling acts only through degrees of freedom associated to element faces. Therefore, the coupling is less strong than for standard elements, and even for the same number of overall degrees of freedom, the system matrix is sparser and thereby easier to invert by a direct solver. As the system matrix is indefinite, it is important to choose a direct solver that can handle such systems.
We count the total number of degrees of freedom for the electric quantities per element. In all cases, one degree of freedom is used for the electric potential, while the number of dielectric displacement degrees of freedom varies depending on k and the element shape. For a tetrahedral element of order k, there are 2(k + 1)(k + 2) coupling and k(k − 1)(2k + 5)/6 interior degrees of freedom. For k = 1 this means 12 coupling degrees of freedom only, for k = 2 this adds up to 24 coupling and 3 interior degrees of freedom. For a prismatic element, we have (k + 1)(k + 2) + 3(k + 1) 2 coupling and (k + 2)(k + 1)(2k + 1)/2 + 1 interior degrees of freedom. For a first order element with k = 1 this adds up to 18 coupling and 10 inner degrees of freedom. For a second order element with k = 2 we get 39 coupling and 31 inner degrees of freedom.

Remanent polarization
There are no conditions on the remanent polarization other than P i ∈ P = [L 2 ] 3 . Neither in derivation nor in stability analysis any assumptions on the continuity of P i were made. Therefore, we propose the remanent polarization to be of the same order as the dielectric displacement and discontinuous,

Implementation
All finite elements described in this section are implemented in the open source software package Netgen/NGSolve available from ngsolve.org. Netgen/NGSolve sports a python interface, where one can symbolically enter variational equations or energies. The necessary derivations for iterative methods such as Newton's iteration are then performed automatically. As all finite elements, including the divergence-conforming ones, are predefined in this software package. The variational equation (32), or the minimization problem corresponding to (21) are entered symbolically, using the definitions (10) for the stored energy and (26) for the regularized dissipation functionφ ε . Poisson ratio ν 0.31 el. permittivity S 2.77 × 10 −8 F m −1 piezoelectric coupling d 31 2.74 × 10 −10 m V −1 piezoelectric coupling d 33 5.93 × 10 −10 m V −1 coercive electric field E 0 820 × 10 3 V m −1 saturation polarization P 0 0.24 C m −2 saturation strain S 0 9.3 × 10 −3 shape parameter m 1.4 hardening parameter h 0 714 × 10 3 m F −1

Non-proportional loading
As a first numerical example, we consider the repolarization of a piezoelectric cube. This benchmark problem was originally proposed by [19], who provide measurement results. In their experiment, a cube was cut from pre-polarized PZT-5H, such that the axes of the cube are at an angle θ compared to the original polarization direction. Then, the cube is electroded and re-polarized by applying an electric field aligned with the cube's z axis. The change in dielectric displacement at the electrodes was measured with growing prescribed electric field. This example has been widely used as a benchmark example for non-proportional loading e.g. by [34,28,1,43]. In our computations, we assume the cube to be of size 5 mm. We use the material parameters for PZT-5H as provided by [39], see Table 1. The material is pre-polarized at angle θ ∈ [0 • , 180 • ]. The regularization parameter for the dissipation function was chosen as ε = 10 −4 P 0 . The free energy function was not regularized.
Concerning this initial state, [43] point out that it is not consistent if θ = 0 • , 180 • , as then surface charges of size ρ P = P i · n occur at the non-electroded boundaries. In this work, we assume one of their remedy proposals, where these charges are modeled as non-moveable surface charges. This way, on the non-electroded boundaries we have constant surface charges of D · n = ρ P .
In a first computation, we perform repolarization in 20 load steps, where a total electric field of 2E 0 is applied. We discretize the cube using an unstructured tetrahedral mesh consisting of 745 elements. We choose the finite element order k = 1, which means second order displacement elements and first order dielectric/remanent polarization elements. In total, we arrive at 18733 degrees of freedom. In Figure 1 the change in dielectric displacements as measured from the finite element model is provided, together with the original data as taken from [19]. A good qualitative resemblance is observed.
In a second computation, we analyze how the loadstep size affects the results. Comparison to the original measurement data taken from [19].
We consider now the case of θ = 90 • . We observe that even large load increments lead to accurate results, see Figure 2. In Table 2, the relative error of the dielectric displacement difference measured at an applied electric field of 2E 0 is provided for different numbers of load steps. To this end, the result from applying 40 loadsteps was used as a reference. We see that using only one load step in this highly nonlinear problem led to a relative error below one percent. In this case, a total number of 35 Newton iterations with line search were needed.
number of loadsteps rel. error 20 0.0001241 10 0.0005419 5 0.0018539 1 0.0099080 Table 2: Relative error measured for the dielectric displacement difference at fully applied electric field 2E 0 for different load step sizes.

Hexahedron with cylindrical hole
This example is a benchmark example taken from [39]. A hexahedron of size 20 mm × 20 mm × 6 mm made of PZT-5H sports a cylindrical hole of diameter 4 mm. Two opposite faces of the hexahedron are electroded, and an electric potential is applied. Material properties are the same as in the previous example and collected in Table 1. The regularization parameter for the dissipation function was chosen as ε = 10 −4 P 0 . The free energy function was not regularized. Symmetry of the problem is used, such that the geometry is reduced to one eighth. In the original reference, the potential was applied in five load steps, and iteration counts were provided. We use a prismatic mesh with one prism over the height, which was also done in the original reference. We provide iteration counts for the first and second order method in Table 3. These counts are higher than those listed by [39], however we only needed load sub-steps in the second order method in the last load increment. Remanent polarization and strain at ∆V = 15 kV are depicted in   aligned at the top and bottom at an angle of 90 • . The computation of macroscopic material constants resembling the mircostructure has been addressed in the literature in different contributions. Analytic mixing rules for a simplified setup were provided by [13]. In [12], the authors address the fact that the geometric setup of the d 33 MFC leads to non-uniform polarization directions in three dimensions. However, they did not simulate the polarization process, but assumed the remanent polarization to be of uniform absolute value P 0 , but in direction of the electric field. [40] presented a first simulation of the polarization process, but using a different model from the one proposed in this work, and only in two space dimensions. Our contribution is now a complete simulation of the poling process, where a poling electric field of size 4E 0 is applied to the interdigitated electrodes. Afterwards, using this remanent polarization field, a small voltage is applied. Strain and dielectric displacement computed in this linear problem then lead to macroscopic d 33 , d 32 and T 33 constants. A sketch of the unit cell as used by [12] is provided in Figure 5. Conforming with [40], we use W = 413.90 µm for the width of an epoxy/PZT fibre set. The width of the PZT fibre w p is then defined by the volume ratio ρ of the fibre. We use ρ = 0.8, which implies w p = ρW = 331.12 µm. Height of fibre and electrode are also taken from [40], and set to h p = 206.40 µm and h e = 18 µm. To compare our results to the values provided by [12], we use their settings w e = h p for the width, and L/w e = 6 for the distance of the copper electrodes. This implies w e = 206.40 µm and L = 1238.40 µm.
Due to symmetry of the geometry and the applied load, which is an applied voltage to the electrodes, we may consider an eigth of the geometry only. Boundary conditions are defined in Figure 6. While the displacement boundary conditions are straightforward for this symmetric problem, we shortly comment on the electric boundary conditions. On the copper electrodes, we assume a constant electric potential of ϕ = ±∆V . This is realized in the following way in our mixed finite element formulation: dielectric displacements are not modelled in the electrode volume, but only in the PZT and epoxy parts. On the electrode surfaces, which are now surfaces of the electric domain, the potential boundary condition ϕ = ± 1 2 ∆V is applied via a surface integral in the external work statement (6). The symmetry plane between electrodes is considered as grounded due to symmetry of the electric potential, which means a (zero) volume integral is added in theory. On all other surfaces, surface charges D · n are set to zero. Poisson ratio ν 0.41 el. permittivity T 1900 0 piezoelectric coupling d 31 −1.85 × 10 −10 m V −1 piezoelectric coupling d 33 4.40 × 10 −10 m V −1 coercive electric field E 0 1150 × 10 3 V m −1 saturation polarization P 0 0.255 C m −2 saturation strain S 0 5.5 × 10 −3 hardening parameter h 0 1/(40 000 0 ) Table 4: Homogenization of macro fibre composite: material parameters for fibres made from PZT-5A.
The choice of material parameters was conducted in the following way: All constants that could be re-used from the linear model by [13,12] were used. This includes the linear properties of epoxy (Y E = 2.9 GPa, ν = 0.3, = 4.25 0 ) and copper (Y E = 117 GPa, ν = 0.31). In this example, we assume that the stiffness at constant electric field c E is constant and isotropic, and defined by Young's modulus Y E and Poisson ratio ν. The permittivity at constant strain S is assumed isotropic, but is not constant. In linear theory, it is smaller than the permittivity at constant stress T . In this contribution, it is assumed to be a convex combination of T = 1900 0 for zero polarization, and 900 0 for fully polarized material. Further, the additive part of the free energy functionψ i is chosen such that its derivative needed in implementations readŝ The saturation polarization P 0 and hardening parameter h 0 were chosen such that the hysteresis measurements depicted in [18, Figure 7(d)] could be reproduced. The remanent saturation strain S 0 was chosen in accordance with the measurements provided by [15]. The regularization parameter for the dissipation function was chosen as ε = 10 −4 P 0 , and to avoid convergence problems at the singularity at the electrode alsoψ i was regularized using the same parameter. All material parameters are also provided in Table 4. In Figure 7 the absolute value at full poling electric field 3E 0 is depicted. Figures 9 and 8 show the strain distributions of S yy and S zz , respectively. After the poling field was removed, a linear computation was performed to compute the homogenized material constants. To this end, an electric potential was applied to the electrodes, and the macroscopic material constants d 32 , d 33 and T 33 were computed by the following formulae provided by [12], The average values ofS yy andS zz were derived from the constant displacements arising at the boundaries, compare Figure 6,S yy = 2c y /W andS zz = 2c z /L. The average value of the electric field was chosen according to [12] asĒ z = ∆V /L. The average value of the dielectric displacement was chosen as the total charge at one electrode divided by the cross section area W H. Due to Gauss' law and the insulation boundary conditions on all outer surfaces, this value is equal to the average of D·n over the symmetric boundary at L/2. The following values were found in our computations, We note that the piezoelectric constants are lower than those estimated by [12]. There are several facts possibly contributing to this outcome: first, the absolute value of the polarization always stays below saturation polarization. It decreases slightly more as the poling electric field is reduced in the used material law, as the remanent polarization at zero voltage is below the saturation polarization also in one-dimensional settings. Moreover, the polarization is lower close to the fibre boundaries, where the interface to the epoxy matrix is located. And last, in the zones below the electrodes, the electric field and polarization are not oriented in z direction, but rather in vertical x direction. This last issue has also been discussed in the finite element simulation by [12].

Conclusion
In this contribution, we have presented energy-based model of the ferroelectric polarization process. In an energy-based setting, dielectric displacement and strain (or displacement) are the primary independent unknowns, while the remanent polarization vector is added as an internal variable. Electric field and stress are evaluated through post-processing utilizing the material law. The model is governed by two constitutive functions: the free energy function and the dissipation function. We provide standard choices from the literature for the free energy, and motivate our choice of dissipation function. As it is nondifferentiable, we proposed to regularize the problem. Then, a variational equation can be posed, which is subsequently discretized using conforming finite elements for each quantity. One section is devoted to the finite elements, as non-standard continuity is needed for the dielectric displacement to obtain a conforming method. We provide corresponding finite elements, that are chosen such that Gauss' law of zero charges is satisfied exactly. The discretized variational equations are solved for all unknowns at once in a single Newton iteration. We present numerical examples gained in the open source software package Netgen/NGSolve.