---
abstract: '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.' 
author:
- 'Astrid S. Pechstein, Martin Meindlhumer and Alexander Humer'
date: January 2020
title: 'High-order mixed finite elements for an energy-based model of the polarization process in ferroelectric materials'
---

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-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.  [@HuberFleckLandisMcMeeking:1998].  Finite element methods based on microscopic models date back to Chen and Lynch [@ChenLynch:1998].  With increasing computing power, methods using one finite element per grain have become feasible, we cite non-exhaustively [@ArockiarajanSivakumarSansour:2009; @LangeRicoeur:2015; @JayendiranGanapathiZineb:2016; @LobanovSemenov:2019]. 

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 [@HwangEtal:1995; @HughesWen:1995; @HegewaldEtal:2008].  Non-linear finite elements for the latter model have been proposed by [@KKaltenbacherHegewaldLerch:2010]. 

A thermodynamically consistent model was first presented in the series of papers by the group around Maugin [@BassiounyGhalebMaugin:1988a; @BassiounyGhalebMaugin:1988b; @BassiounyMaugin:1989a; @BassiounyMaugin:1989b].  Kamlah [@Kamlah:2001] 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 [@KamlahBoehle:2001].  More recently, Elhadrouz et al.  [@ElhadrouzZinebPatoor:2005] 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 [@McMeekingLandis:2002] as well as Schröder and Romanowski [@SchroederRomanowski:2005] proposed formulations where the remanent polarization strain depends one-to-one on the remanent polarization.  Nevertheless, these models are capable of mechanic depolarization, see the numerical examples from [@McMeekingLandis:2002].  Linnemann et al.  [@LinnemannKlinkelWagner:2009] added magnetostriction in their work.  Another approach with an independent polarization strain has been proposed by Landis [@Landis:2002].  Also Klinkel [@Klinkel:2006] uses two independent remanent quantities remanent electric field and remanent strain in his model.  Laxman et al.  [@LaxmanManiprakashArockiarajan:2018] additionally modeled creep behavior in their approach. 

Miehe et al.  [@MieheRosatoKiefer:2011] presented a variational framework for electro-magneto-mechanics.  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 [@SandsGuz:2013] describe a one-dimensional model based on the free energy.  Semenov et al.  [@SemenovLiskowskyBalke:2010] 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.  [@StarkNeumeisterBalke:2016I; @StarkNeumeisterBalke:2016II] 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 [@PechsteinMeindlhumerHumer:2019arxiv].  This concept has been used in the mathematical analysis of other nonlinear problems in mechanics, such as elasto-plasticity by Han and Reddy [@HanReddy:1999] or contact mechanics by Sofonea and Matei [@SofoneaMatei:2011].  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 [@MieheRosatoKiefer:2011] 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 $\Omega$, which is in total or in part made of ferroelectric material.  We use ${\mathbf{u}}: \Omega \times [0,T] \to \mathbb{R}^3$ to denote the displacement field, and use ${\mathbf{S}}= \frac12 (\nabla {\mathbf{u}}+ (\nabla {\mathbf{u}})^T)$ for the linearized strain tensor. 

In an electrostatic regime on a simply connected domain, the electric field ${\mathbf{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 $\varphi$ with ${\mathbf{E}}=  - \nabla \varphi$. The dielectric displacement vector ${\mathbf{D}}$ satisfies Gauss’ law of zero charge in the material, $$\begin{aligned}
 \label{eq:Gauss}
{\operatorname{div}}{\mathbf{D}}&= 0 &&\text{in } \Omega.\end{aligned}$$ Last, we use the remanent polarization ${\mathbf{P}^{i}}$ as an internal variable. 

The reversible material response is characterized by the internal energy $\Psi$ or its density $\psi$, which shall depend on strain, dielectric displacement as well as remanent polarization vector, $$\begin{aligned}
\psi({\mathbf{x}}, t) &= \hat \psi({\mathbf{S}}({\mathbf{x}},t), {\mathbf{D}}({\mathbf{x}},t), {\mathbf{P}^{i}}({\mathbf{x}},t)),\\
\Psi(t) &= \int_\Omega \psi({\mathbf{x}},t)\, d{\mathbf{x}}.\end{aligned}$$

Total stress ${\mathbf{T}}$ and electric field ${\mathbf{E}}$ are conjugate to strain ${\mathbf{S}}$ and dielectric displacement ${\mathbf{D}}$ in the following sense, $$\begin{aligned}
{\mathbf{T}}&= \frac{\partial \hat\psi}{\partial {\mathbf{S}}}, & {\mathbf{E}}&= \frac{\partial \hat \psi}{\partial {\mathbf{D}}}. 
\label{eq:defTE}\end{aligned}$$ Conjugate to the remanent polarization is the (negative) electric driving force $\hat {\mathbf{E}}$, $$\begin{aligned}
\hat {\mathbf{E}}&= -\frac{\partial \hat \psi}{\partial {\mathbf{P}^{i}}}. 
\label{eq:defhatE}\end{aligned}$$

We denote by $P^{ext} = \frac{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 ${\mathbf{f}}$ in $\Omega$, a surface load ${\mathbf{t}}$ on $\Gamma_{t}$ and an applied electric potential $V_0$ on $\Gamma_{el}$, the power from external sources is given by $$\begin{aligned}
 \label{eq:externalwork}
P^{ext} &= \int_\Omega {\mathbf{f}}\cdot \dot{\mathbf{u}}\, d{\mathbf{x}}+ \int_{\Gamma_{t}} {\mathbf{t}}\cdot \dot{\mathbf{u}}\, ds - \int_{\Gamma_{el}} V_0\, \dot{\mathbf{D}}\cdot {\mathbf{n}}\, ds.\end{aligned}$$

A reversible process is fully characterized by internal energy $\Psi$ 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.  [@MieheRosatoKiefer:2011], the dissipative behavior can be described by the dissipation function $\Phi$, 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, $$\begin{aligned}
\phi({\mathbf{x}},t) &= \hat \phi(\dot {\mathbf{P}^{i}}({\mathbf{x}},t)),\\
\Phi(t) &= \int_\Omega \phi({\mathbf{x}},t)\, d{\mathbf{x}}.\end{aligned}$$ 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)$, $$\begin{aligned}
\begin{split}
\Psi(t_1) - \Psi(t_0)& + \int_{t_0}^{t_1}\!\!  \Phi(t)\, dt - \int_{t_0}^{t_1}\!\!\!  P^{ext}(t)\, dt \to \min_{{\mathbf{u}}, {\mathbf{D}}}\\
\text{subject to } & \left\{\begin{array}{l}{\operatorname{div}}{\mathbf{D}}= 0 \text{ and}\\
\text{boundary conditions on } {\mathbf{u}}, {\mathbf{D}}\cdot{\mathbf{n}}.\end{array} \right. 
\end{split} \label{eq:minimization_cont}\end{aligned}$$ For a more detailed introduction, we refer to the extensive derivations originally provided by [@MieheRosatoKiefer:2011]. 

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 ${\mathbf{c}}^{D}, {\mathbf{h}}$ and ${\boldsymbol{\beta}}^{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$. $$\begin{aligned}
\begin{split}
\hat \psi =& \ \hat \psi^{r} +\hat  \psi^{i},\\
\hat \psi^r =&\ \frac12 ({\mathbf{S}}- {\mathbf{S}}^{i}({\mathbf{P}^{i}})):{\mathbf{c}}^{D}:({\mathbf{S}}- {\mathbf{S}}^{i}({\mathbf{P}^{i}}))\\
& -({\mathbf{S}}- {\mathbf{S}}^{i}({\mathbf{P}^{i}})):{\mathbf{h}}\cdot({\mathbf{D}}- {\mathbf{P}^{i}}) \\
& +\frac12 ({\mathbf{D}}- {\mathbf{P}^{i}})\cdot{\boldsymbol{\beta}}^{S}\cdot({\mathbf{D}}-{\mathbf{P}^{i}}),\\
\hat \psi^{i} =&\ \frac{h_0\, P_0^2}{(m-1)(m-2)} \left( 1-\frac{|{\mathbf{P}^{i}}|}{P_0} \right)^{2-m}\\
& - \frac{h_0\, P_0^2}{m-1}\frac{|{\mathbf{P}^{i}}|}{P_0}. 
\end{split} \label{eq:defpsi}\end{aligned}$$ We use a one-to-one relationship between remanent polarization and polarization strain as suggested by [@McMeekingLandis:2002], $$\begin{aligned}
S^{i}_{kl}({\mathbf{P}^{i}}) &= \frac{S_0}{2 P_0^2} (3 P^{i}_k P^{i}_l - \delta_{kl} P^{i}_j P^{i}_j).\end{aligned}$$ Note that this choice is volume preserving. 

Choice of dissipation function
------------------------------

Apart from the free energy density $\hat \psi$, the dissipation function $\hat \phi$ has to be specified to obtain a valid model through the minimization problem .  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 $$\begin{aligned}
D = \dot {\mathbf{P}^{i}}\cdot \hat {\mathbf{E}}\geq 0.\end{aligned}$$ The choice of dissipation function $\hat \phi$ 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 $\hat {\mathbf{E}}$ exceeds the coercive field $E_0$, $$\begin{aligned}
|\hat {\mathbf{E}}| - E_0 & \leq 0, & |\dot {\mathbf{P}^{i}}| &\geq 0, & (|\hat {\mathbf{E}}| - E_0)|\dot {\mathbf{P}^{i}}| = 0.\end{aligned}$$ The dissipation function is chosen such that dissipation is maximized over all possible values for $\hat {\mathbf{E}}$, $$\begin{aligned}
\hat \phi(\dot {\mathbf{P}^{i}}) = \sup_{\hat {\mathbf{E}}, |\hat {\mathbf{E}}|\leq E_0} \dot {\mathbf{P}^{i}}\cdot \hat {\mathbf{E}}.\end{aligned}$$ This supremum is of course taken for $\hat {\mathbf{E}}= E_0 {\mathbf{e}}_P$, with ${\mathbf{e}}_P = {\mathbf{P}^{i}}/|{\mathbf{P}^{i}}|$ the unit vector in direction of remanent polarization, and for $\hat {\mathbf{E}}= \mathbf{0}$ in case $|\dot {\mathbf{P}^{i}}| = 0$. Thus, the dissipation function can be simplified to read $$\begin{aligned}
\hat \phi(\dot {\mathbf{P}^{i}}) = E_0 |\dot {\mathbf{P}^{i}}|.\end{aligned}$$ Note that this dissipation function is positively homogeneous (compare ).  Such a dissipation function leads to rate-independent behavior, see @MieheRosatoKiefer:2011. 

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 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 $\Gamma_{fix}$ are included, such that we are searching for $$\begin{aligned}
{\mathbf{u}}\in {\mathbb{V}}&:= \{{\mathbf{v}}\in H^1(\Omega): {\mathbf{v}}= 0 \text{ on } \Gamma_{fix}\}.\end{aligned}$$ The choice of the dielectric displacement space is less standard.  We resort to the $H({\operatorname{div}})$ space of $L^2$ vector fields with weak divergence.  Gauss’ law is then included explicitely.  Additionally, insulation boundary conditions of vanishing surface charges ${\mathbf{D}}\cdot {\mathbf{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 that voltage boundary conditions on electrodes are natural.  Summing up, we search for $$\begin{aligned}
\begin{split}
{\mathbf{D}}\in {\mathbb{D}}_0 &:= \{{\mathbf{d}}\in {\mathbb{D}}: {\operatorname{div}}{\mathbf{d}}= 0\},\\
\text{with } {\mathbb{D}}& = \{{\mathbf{d}}\in [L^2]^3: {\operatorname{div}}{\mathbf{d}}\in L^2,\ {\mathbf{d}}\!\cdot\!{\mathbf{n}}= 0 \text{ on } \Gamma_{ins}\}\\
& \subset H({\operatorname{div}}). 
\end{split}\end{aligned}$$

The remanent polarization does not sport any smoothness, nor does it allow for boundary conditions.  Therefore we choose the vector-valued Lebesgue space $[L^2]^3$, $$\begin{aligned}
{\mathbf{P}^{i}}\in {\mathbb{P}}&:= [L^2]^3.\end{aligned}$$

We will now formulate the incremental variational problem for some small but finite time increment $\Delta t$, starting from time $t_0$ to $t_1 = t_0+\Delta t$. To this end, we assume ${\mathbf{u}}_0 = {\mathbf{u}}(t_0)$, ${\mathbf{D}}_0 = {\mathbf{D}}(t_0)$ and ${\mathbf{P}^{i}}_0 = {\mathbf{P}^{i}}(0)$ known.  We are interested in finding the increments $\Delta {\mathbf{u}}, \Delta {\mathbf{D}}$ and $\Delta {\mathbf{P}^{i}}$ such that ${\mathbf{u}}(t_1) = {\mathbf{u}}_0 + \Delta {\mathbf{u}}$, ${\mathbf{D}}(t_1) = {\mathbf{D}}_0 + \Delta {\mathbf{D}}$ and ${\mathbf{P}^{i}}(t_1) = {\mathbf{P}^{i}}_0 + \Delta {\mathbf{P}^{i}}$. Miehe et al.  [@MieheRosatoKiefer:2011] stated that these increments minimize the potential increment, $$\begin{aligned}
\begin{split}
&(\Delta {\mathbf{u}}, \Delta {\mathbf{D}}, \Delta {\mathbf{P}^{i}}) =\\
&= \arg\min \Big( \Psi(t_1) - \Psi(t_0) + \Delta t \Phi(\Delta {\mathbf{P}^{i}}/\Delta t) - W^{ext}_{t_0\to t_1}\Big)
\\
& \quad \text{subject to } \Delta {\mathbf{u}}\in {\mathbb{V}}, \Delta {\mathbf{D}}\in {\mathbb{D}}_0, \Delta {\mathbf{P}^{i}}\in {\mathbb{P}}. 
\end{split}
\label{eq:minimization_deltat_ratedep}\end{aligned}$$ In the special case of a rate-independent response, where the dissipation function is positively homogeneous, $$\begin{aligned}
 \label{eq:poshom}
\Phi(\alpha \dot {\mathbf{P}^{i}}) &= \alpha \Phi(\dot {\mathbf{P}^{i}}) & \text{for } \alpha > 0,\end{aligned}$$ the time step size cancels out in .  The parameter $t$ acts as pseudo-time then, the minimization problem reads $$\begin{aligned}
\begin{split}
&(\Delta {\mathbf{u}}, \Delta {\mathbf{D}}, \Delta {\mathbf{P}^{i}}) =\\
&= \arg\min \Big( \Psi(t_1) - \Psi(t_0) + \Phi(\Delta {\mathbf{P}^{i}}) - W^{ext}_{t_0\to t_1}\Big)
\\
&\quad \text{subject to }   \Delta {\mathbf{u}}\in {\mathbb{V}}, \Delta {\mathbf{D}}\in {\mathbb{D}}_0, \Delta {\mathbf{P}^{i}}\in {\mathbb{P}}. 
\end{split}
\label{eq:minimization_deltat}\end{aligned}$$

Variational inequality
----------------------

As the dissipation function is not differentiable at $\dot {\mathbf{P}^{i}}= \mathbf{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. [@HanReddy:1999]) or contact problems (see e.g. [@SofoneaMatei:2011]), a thorough analysis of the arising time-dependent variational inequality has been carried out in [@PechsteinMeindlhumerHumer:2019arxiv].  In this contribution, we state the variational inequality which is equivalent to the optimization problem .  We use the definitions for total stress ${\mathbf{T}}$ and electric field ${\mathbf{E}}$ and for the electric driving force $\hat {\mathbf{E}}$ evaluated at $t_1$, $$\begin{aligned}
{\mathbf{T}}_{t_1} &= \frac{\partial \hat\psi}{\partial {\mathbf{S}}}({\mathbf{S}}({\mathbf{u}}_0+\Delta {\mathbf{u}}),{\mathbf{D}}_0 + \Delta {\mathbf{D}}, {\mathbf{P}^{i}}_0 + \Delta {\mathbf{P}^{i}}),\\
{\mathbf{E}}_{t_1} &= \frac{\partial \hat\psi}{\partial {\mathbf{D}}}({\mathbf{S}}({\mathbf{u}}_0+\Delta {\mathbf{u}}),{\mathbf{D}}_0 + \Delta {\mathbf{D}}, {\mathbf{P}^{i}}_0 + \Delta {\mathbf{P}^{i}}),\\
\hat {\mathbf{E}}_{t_1} &= -\frac{\partial \hat\psi}{\partial {\mathbf{P}^{i}}}({\mathbf{S}}({\mathbf{u}}_0+\Delta {\mathbf{u}}),{\mathbf{D}}_0 + \Delta {\mathbf{D}}, {\mathbf{P}^{i}}_0 + \Delta {\mathbf{P}^{i}}).\end{aligned}$$ Then, $\Delta {\mathbf{u}}\in {\mathbb{V}}$, $\Delta {\mathbf{D}}\in {\mathbb{D}}_0$ and $\Delta {\mathbf{P}^{i}}\in {\mathbb{P}}$ satisfy the following variational inequality $$\begin{aligned}
\begin{split}
\int_\Omega \big( {\mathbf{T}}_{t_1}:({\mathbf{S}}(\tilde {\mathbf{u}})-{\mathbf{S}}(\Delta{\mathbf{u}})) + {\mathbf{E}}_{t_1} \cdot (\tilde {\mathbf{D}}- \Delta {\mathbf{D}}) &\\
- \hat {\mathbf{E}}_{t_1} \cdot ( \tilde {\mathbf{P}^{i}}- \Delta {\mathbf{P}^{i}})
+ \hat \phi(\tilde {\mathbf{P}^{i}}) - \hat \phi(\Delta {\mathbf{P}^{i}}) \big) d{\mathbf{x}}&\\
\geq  \tilde W^{ext}_{t_0\to t_1} - W^{ext}_{t_0\to t_1} &\\
\text{for all } \tilde {\mathbf{u}}\in {\mathbb{V}}, \tilde {\mathbf{D}}\in {\mathbb{D}}_0, \tilde {\mathbf{P}^{i}}\in {\mathbb{P}}.& 
\end{split} \label{eq:varin}\end{aligned}$$

Regularization
--------------

In the framework of elasto-plasticity, Han and Reddy [@HanReddy:1999] 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 turns into a nonlinear variational equation, which can be solved by standard nonlinear solution algorithms. 

We propose to use the regularized dissipation function $$\begin{aligned}
\hat \phi_\varepsilon &= \left\{ \begin{array}{ll}
    E_0\big(|\dot{\mathbf{P}^{i}}| - \frac{\varepsilon}{2}\big) & |\dot {\mathbf{P}^{i}}| \geq \varepsilon\\
    \frac1{2\varepsilon} E_0 |\dot {\mathbf{P}^{i}}|^2 & |\dot {\mathbf{P}^{i}}| < \varepsilon. 
    \end{array} \right. 
    \label{eq:defphieps}\end{aligned}$$ This choice has been analyzed to work well for elasto-plasticity in [@HanReddy:1999].  Indeed, according to this reference, we expect the error due to regularization to be of order $\sqrt{\varepsilon}$.

For completeness we provide the variational equation to be solved in each time step.  We denote the (continuous) derivative of the regularized dissipation function $\hat \phi_\varepsilon$ with respect to the remanent polarization rate as $\hat \phi_\varepsilon'$. The updates $\Delta {\mathbf{u}}\in {\mathbb{V}}$, $\Delta {\mathbf{D}}\in {\mathbb{D}}_0$ and $\Delta {\mathbf{P}^{i}}\in {\mathbb{P}}$ satisfy $$\begin{aligned}
\begin{split}
\int_\Omega \big( {\mathbf{T}}_{t_1}:\,\delta {\mathbf{S}}+ {\mathbf{E}}_{t_1} \cdot \delta {\mathbf{D}}- \hat {\mathbf{E}}_{t_1} \cdot \delta {\mathbf{P}^{i}}&\\
+ \hat \phi'_\varepsilon(\Delta {\mathbf{P}^{i}})\cdot \delta {\mathbf{P}^{i}}\big) &d{\mathbf{x}}=  \delta W^{ext}_{t_0\to t_1}  \\
\text{for all } \tilde {\mathbf{u}}\in {\mathbb{V}}, \tilde {\mathbf{D}}\in {\mathbb{D}}_0, \tilde {\mathbf{P}^{i}}\in {\mathbb{P}}.& 
\end{split} \label{eq:vareq}\end{aligned}$$

Another issue for regularization is the free energy function $\psi$, or its additive part $\psi^{i}$. The formula proposed in is unbounded as the saturation polarization is approached, $$\begin{aligned}
\hat \psi^{i}({\mathbf{P}^{i}}) \to \infty &\text{ and } (\hat \psi^{i})'({\mathbf{P}^{i}})\cdot \delta{\mathbf{P}^{i}}\to \infty \\
&\text{ as } |{\mathbf{P}^{i}}| \to P_0.\end{aligned}$$ This unbounded potential acts as a penalty term as the polarization approaches saturation, where also $ (\hat \psi^{i})'({\mathbf{P}^{i}})$ diverges to infinity.  As analyzed in [@PechsteinMeindlhumerHumer:2019arxiv], this fact does not harm solvability of the update equation itself, but it may harm convergence.  Therefore, if convergence problems are observed, we propose to modify the energy density $\hat \psi^{i}$ in a way that $ (\hat \psi^{i}_\varepsilon)'({\mathbf{P}^{i}})$ becomes large but stays bounded.  We use the following regularized expression $ (\hat \psi^{i}_\varepsilon)'$, $$\begin{aligned}
(\hat \psi^{i}_\varepsilon)'({\mathbf{P}^{i}}) &= \left\{ \begin{array}{l}
\left( \frac{h_0 P_0^m}{m-1}(P_0-|{\mathbf{P}^{i}}|)^{1-m} -1\right) {\mathbf{P}^{i}}/|{\mathbf{P}^{i}}| \\
\text{  if } |{\mathbf{P}^{i}}| \leq P_0-\varepsilon,\\
\left(\frac{h_0 P_0^m}{m-1}(\varepsilon)^{1-m} -1\right) {\mathbf{P}^{i}}/(P_0-\varepsilon) \\
\text{ if } |{\mathbf{P}^{i}}| > P_0-\varepsilon. 
\end{array} \right.  \end{aligned}$$

Gauss’ law
----------

Gauss’ law of zero charges implies that the dielectric displacement field is divergence free inside the domain $\Omega$. So far, this condition has been included in the space of admissible updates ${\mathbb{D}}_0$ without further ceremony.  As Gauss’ law is a linear and continuous restriction, the space ${\mathbb{D}}_0$ is a linear and continuous subspace of ${\mathbb{D}}$ and $H({\operatorname{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 $\varphi$. In the minimization problem, we propose to use an augmented energy $\Psi_0$, which assumes infinity if Gauss’ law is violated.  This augmented energy reads $$\begin{aligned}
\Psi_0 &= \Psi({\mathbf{S}}, {\mathbf{D}}, {\mathbf{P}^{i}}) - \inf_{\varphi \in {\mathbb{Q}}} \int_\Omega {\operatorname{div}}{\mathbf{D}}\, \varphi\, d{\mathbf{x}}.\end{aligned}$$ Note that we use $\varphi \in {\mathbb{Q}}:= L^2(\Omega)$, 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 .  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 $\Delta {\mathbf{u}}\in {\mathbb{V}}$, $\Delta {\mathbf{D}}\in {\mathbb{D}}$ and $\Delta {\mathbf{P}^{i}}\in {\mathbb{P}}$ as well as the (total) electric potential $\varphi \in {\mathbb{Q}}$ satisfy $$\begin{aligned}
\begin{split}
\int_\Omega \big( {\mathbf{T}}_{t_1}:\,\delta {\mathbf{S}}+ {\mathbf{E}}_{t_1} \cdot \delta {\mathbf{D}}- \hat {\mathbf{E}}_{t_1} \cdot \delta {\mathbf{P}^{i}}&\\
-{\operatorname{div}}(\delta {\mathbf{D}})\, \varphi - {\operatorname{div}}(\Delta {\mathbf{D}}) \, \delta \varphi&\\
+ \hat \phi'_\varepsilon(\Delta {\mathbf{P}^{i}})\cdot \delta {\mathbf{P}^{i}}& \big) d{\mathbf{x}}=  \delta W^{ext}_{t_0\to t_1}  \\
\text{for all } \delta {\mathbf{u}}\in {\mathbb{V}}, \delta {\mathbf{D}}\in {\mathbb{D}}, \delta {\mathbf{P}^{i}}\in {\mathbb{P}}, \delta \varphi \in {\mathbb{Q}}.& 
\end{split} \label{eq:vareq_lag}\end{aligned}$$

The augmented variational equation is equivalent to .  The dielectric displacement update $\Delta {\mathbf{D}}$ is exactly divergence free and thus $\Delta {\mathbf{D}}\in {\mathbb{D}}_0$.

Solvability and stability
-------------------------

The issue of existence of solutions, also to the time-dependent variational inequality , has been treated by the authors in @PechsteinMeindlhumerHumer:2019arxiv.  We recall that a unique solution to the update problem exists, given the following essential assumption of stong monotonicity holds: for all ${\mathbf{u}}, \tilde {\mathbf{u}}\in {\mathbb{V}}$, ${\mathbf{D}}, \tilde {\mathbf{D}}\in {\mathbb{D}}_0$ and ${\mathbf{P}^{i}}, \tilde {\mathbf{P}^{i}}\in {\mathbb{P}}$ as well as corresponding forces ${\mathbf{T}}, \tilde {\mathbf{T}}$, ${\mathbf{E}}, \tilde {\mathbf{E}}$ and $\hat {\mathbf{E}}, \tilde{\hat {\mathbf{E}}}$ there holds $$\begin{aligned}
\int_\Omega \big( &({\mathbf{T}}- \tilde {\mathbf{T}}):({\mathbf{S}}- \tilde {\mathbf{S}}) + ({\mathbf{E}}- \tilde {\mathbf{E}})\cdot({\mathbf{D}}- \tilde {\mathbf{D}}) \nonumber\\
& - (\hat {\mathbf{E}}- \tilde {\hat {\mathbf{E}}})\cdot({\mathbf{P}^{i}}- \tilde {\mathbf{P}^{i}}) \big)\, d{\mathbf{x}}\\
&\geq C \big( \|{\mathbf{u}}- \tilde {\mathbf{u}}\|_{\mathbb{V}}^2 + \|{\mathbf{D}}- \tilde {\mathbf{D}}\|_{\mathbb{D}}^2 + \|{\mathbf{P}^{i}}-\tilde {\mathbf{P}^{i}}\|_{\mathbb{P}}^2\big).\nonumber\end{aligned}$$ This is the case if the reversible part of energy $\Psi^{r}$ is strongly monotone with respect to the reversible parts ${\mathbf{S}}- {\mathbf{S}}^{i}({\mathbf{P}^{i}})$ and ${\mathbf{D}}- {\mathbf{P}^{i}}$, and the additive part $\Psi^{i}$ is strongly monotone with respect to ${\mathbf{P}^{i}}$. This is the case, as long as the linear moduli lead to a positive definite material tensor.  This issue is addressed by [@StarkNeumeisterBalke:2016a].  Conditions ensuring positive definiteness of the additive energy part $\hat \psi^{i}$ can be found in [@BotteroIdiart:2018]. 

Finite elements
===============

In this section we propose a set of consistent finite elements to discretize the variational equation above.  The elements shall be of arbitrary order, by $k \geq 1$ we denote the lowest overall polynomial order.  In the following, we assume that $\mathcal{T}$ is a simplicial mesh of the volume $\Omega$. 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 ${\mathbf{u}}\in {\mathbb{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 ${\mathbb{V}}_h$ is then given by $$\begin{aligned}
 {\mathbb{V}}_h &:= \{{\mathbf{v}}\in \mathcal{C}(\Omega): {\mathbf{v}}|_T \in [P^{k+1}]^3\text{ for } T \in \mathcal{T}\} \subset {\mathbb{V}}. 
 \end{aligned}$$ 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 ${\mathbf{D}}\in {\mathbb{D}}$, i.e. that ${\mathbf{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 ${\mathbf{D}}\cdot {\mathbf{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 [@BoffiBrezziFortin:2013].  As a first proposition, we propose to use elements from the BDM family as introduced by [@BrezziDouglasDuranFortin:1987] (BDM after Brezzi, Douglas and Marini, who first proposed corresponding elements in two dimensions, see [@BrezziDouglasMarini:1986]).  These elements span the space

$$\begin{aligned}
 \begin{split}
 {\mathbb{D}}_h:= \{ &{\mathbf{D}}\in [L^2(\Omega)]^3: {\mathbf{D}}|_T \in [P^k]^3\text{ for } T \in \mathcal{T},\\
 & {\mathbf{D}}\cdot {\mathbf{n}}\text{ cont.}, { \mathbf{D}}\cdot {\mathbf{n}}= 0 \text{ on }\Gamma_{ins}\} \subset {\mathbb{D}}. 
 \end{split}
 \end{aligned}$$

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 ${\mathbb{D}}_h$ such that it becomes a subspace of ${\mathbb{D}}_0$ in a simple, local manner.  However, in context of incompressible Navier-Stokes equation, [@LehrenfeldSchoeberl:2016] have shown that it is possible to do so for the higher-order shape functions.  We propose to use this reduced set, $$\begin{aligned}
\begin{split}
{\mathbb{D}}_{h0}:= \{ &{\mathbf{D}}\in [L^2(\Omega)]^3: {\mathbf{D}}|_T \in [P^k]^3\text{ for } T \in \mathcal{T},\\
& {\operatorname{div}}{\mathbf{D}}|_T \in P^0 \text{ for } T \in \mathcal{T},\\
& {\mathbf{D}}\cdot {\mathbf{n}}\text{ cont.}, { \mathbf{D}}\cdot {\mathbf{n}}= 0 \text{ on }\Gamma_{ins}\} \subset {\mathbb{D}}. 
\end{split}\end{aligned}$$

We see that this way not only the number of unknowns for ${\mathbf{D}}$ drops, but also less Lagrangian multipliers $\varphi$ are needed.  As ${\operatorname{div}}{\mathbf{D}}$ is constant per element, it is enough to use piecewise constant $\varphi$, $$\begin{aligned}
 {\mathbb{Q}}_h &:= \{ \varphi \in L^2(\Omega): \varphi|_T \in P^0\text{ for } T \in \mathcal{T}\}
 \end{aligned}$$ Then, the condition $$\begin{aligned}
 \int_\Omega {\operatorname{div}}{\mathbf{D}}\, \delta \varphi \, &= 0 && \text{for all } \delta \varphi \in {\mathbb{Q}}_h,
 \end{aligned}$$ implies that ${\operatorname{div}}{\mathbf{D}}= 0$ in strong sense for ${\mathbf{D}}\in {\mathbb{D}}_{h0}$.

The finite element space ${\mathbb{D}}_{h0}$ is intrinsically different from standard nodal finite element spaces well-known in mechanics.  While in nodal finite element spaces such as ${\mathbb{V}}_h$ degrees of freedom are function values evaluated at element nodes, a more abstract concept of degrees of freedom as introduced by [@Ciarlet:1978] is needed in the context of $H({\operatorname{div}})$. In our case, degrees of freedom are associated to element interfaces (to ${\mathbf{D}}\cdot {\mathbf{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 ${\mathbf{P}^{i}}\in {\mathbb{P}}= [L^2]^3$. Neither in derivation nor in stability analysis any assumptions on the continuity of ${\mathbf{P}^{i}}$ were made.  Therefore, we propose the remanent polarization to be of the same order as the dielectric displacement and discontinuous, $$\begin{aligned}
 {\mathbb{P}}_h &:= \{ {\mathbf{P}^{i}}\in [L^2(\Omega)]^3: {\mathbf{P}^{i}}|_T \in [P^k]^3\text{ for } T \in \mathcal{T}\}
 \end{aligned}$$

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 , or the minimization problem corresponding to are entered symbolically, using the definitions for the stored energy and for the regularized dissipation function $\hat \phi_\varepsilon$.

Numerical results
=================

Non-proportional loading
------------------------

As a first numerical example, we consider the repolarization of a piezoelectric cube.  This benchmark problem was originally proposed by [@HuberFleck:2001], 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 $\theta$ 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 [@McMeekingLandis:2002; @Landis:2002; @ArockiarajanSivakumarSansour:2009; @StarkNeumeisterBalke:2016II].  In our computations, we assume the cube to be of size $\SI{5}{\milli\meter}$. We use the material parameters for PZT-5H as provided by [@SemenovLiskowskyBalke:2010], see Table \[tab:materialPZT5H\].  The material is pre-polarized at angle $\theta \in [0^\circ, 180^\circ]$. The regularization parameter for the dissipation function was chosen as $\varepsilon=10^{-4}P_0$. The free energy function was not regularized. 

Concerning this initial state, [@StarkNeumeisterBalke:2016II] point out that it is not consistent if $\theta \neq 0^\circ, 180^\circ$, as then surface charges of size $\rho_{P} = {\mathbf{P}^{i}}\cdot{\mathbf{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 ${\mathbf{D}}\cdot {\mathbf{n}}= \rho_P$.

In a first computation, we perform repolarization in 20 load steps, where a total electric field of $2 E_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 \[fig:repoling\_theta\] the change in dielectric displacements as measured from the finite element model is provided, together with the original data as taken from [@HuberFleck:2001].  A good qualitative resemblance is observed. 

  --------------------------------- -----------------------------------------
  Young’s modulus $Y_E$             $\SI{61e9}{\newton\per\meter\squared}$
  Poisson ratio $\nu$               0.31
  el.  permittivity $\epsilon^{S}$   $\SI{2.77e-8}{\farad\per\meter}$
  piezoelectric coupling $d_{31}$   $\SI{2.74e-10}{\meter\per\volt}$
  piezoelectric coupling $d_{33}$   $\SI{5.93e-10}{\meter\per\volt}$
  coercive electric field $E_0$     $\SI{820e3}{\volt\per\meter}$
  saturation polarization $P_0$     $\SI{0.24}{\coulomb\per\meter\squared}$
  saturation strain $S_0$           $\SI{9.3e-3}{}$
  shape parameter $m$               1.4
  hardening parameter $h_0$         $\SI{714e3}{\meter\per\farad}$
  --------------------------------- -----------------------------------------

  : Material parameters for PZT-5H used for repolarization, [@SemenovLiskowskyBalke:2010].[]{data-label="tab:materialPZT5H"}

![ Repoling of piezoelectric cube: The change of dielectric displacement is measured with growing applied electric field for different angles $\theta$. Comparison to the original measurement data taken from [@HuberFleck:2001].[]{data-label="fig:repoling_theta"}](repolarization_theta){width="0.9\columnwidth"}

In a second computation, we analyze how the loadstep size affects the results.  We consider now the case of $\theta = 90^\circ$. We observe that even large load increments lead to accurate results, see Figure \[fig:repoling\_stepsize\].  In Table \[tab:relerrors\_repoling\], the relative error of the dielectric displacement difference measured at an applied electric field of $2 E_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. 

![ Repoling of piezoelectric cube: Comparison of different loadstep sizes for an initial polarization angle of $90^\circ$.[]{data-label="fig:repoling_stepsize"}](repolarization_stepsize_90deg){width="0.9\columnwidth"}

    number of loadsteps   rel.  error
  --------------------- ------------
                     20    0.0001241
                     10    0.0005419
                      5    0.0018539
                      1    0.0099080

  : Relative error measured for the dielectric displacement difference at fully applied electric field $2 E_0$ for different load step sizes.[]{data-label="tab:relerrors_repoling"}

Hexahedron with cylindrical hole
--------------------------------

This example is a benchmark example taken from [@SemenovLiskowskyBalke:2010].  A hexahedron of size $\SI{20}{\milli\meter}\times\SI{20}{\milli\meter}\times\SI{6}{\milli\meter}$ made of PZT-5H sports a cylindrical hole of diameter $\SI{4}{\milli\meter}$. 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 \[tab:materialPZT5H\].  The regularization parameter for the dissipation function was chosen as $\varepsilon=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 \[tab:itcount\].  These counts are higher than those listed by [@SemenovLiskowskyBalke:2010], however we only needed load sub-steps in the second order method in the last load increment.  Remanent polarization and strain at $\Delta V = \SI{15}{\kilo\volt}$ are depicted in Figure \[fig:15kv\]. 

                                 order $k=1$   order $k=2$
  ---------------------------- ------------- -------------
      $0\to\SI{5}{\kilo\volt}$             8            10
      $5\to\SI{8}{\kilo\volt}$            19            31
      $8\to\SI{9}{\kilo\volt}$            33            34
     $9\to\SI{10}{\kilo\volt}$            21            21
    $10\to\SI{15}{\kilo\volt}$            48     13+11+8+7

  : Hexahedron with cylindrical hole: Iteration counts in Newton’s method with linesearch for the five load steps.[]{data-label="tab:itcount"}

![ Hexahedron with cylindrical hole: Remanent polarization $|{\mathbf{P}^{i}}|$ and strain component $S_{zz}$ at $\Delta V = \SI{15}{\kilo\volt}$.[]{data-label="fig:15kv"}](cubewithhole_PI_k1 "fig:"){width="\columnwidth"} ![ Hexahedron with cylindrical hole: Remanent polarization $|{\mathbf{P}^{i}}|$ and strain component $S_{zz}$ at $\Delta V = \SI{15}{\kilo\volt}$.[]{data-label="fig:15kv"}](cubewithhole_PI_k2 "fig:"){width="\columnwidth"}

![ Hexahedron with cylindrical hole: Remanent polarization $|{\mathbf{P}^{i}}|$ and strain component $S_{zz}$ at $\Delta V = \SI{15}{\kilo\volt}$.[]{data-label="fig:15kv"}](cubewithhole_szz_k1 "fig:"){width="\columnwidth"} ![ Hexahedron with cylindrical hole: Remanent polarization $|{\mathbf{P}^{i}}|$ and strain component $S_{zz}$ at $\Delta V = \SI{15}{\kilo\volt}$.[]{data-label="fig:15kv"}](cubewithhole_szz_k2 "fig:"){width="\columnwidth"}

Poling of a $d_{33}$ macro fiber composite
------------------------------------------

As an advanced numerical example, we consider a $d_{33}$ macro fibre composite ($d_{33}$ MFC).  Such a $d_{33}$ MFC is commercially available e.g. from Smart Materials Corp., Sarasota, FL.  This MFC consists of thin (microscopic) fibres of PZT which are aligned in an epoxy matrix.  Interfingering copper electrodes are aligned at the top and bottom at an angle of $90^\circ$. 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 [@DeraemaekerEtal:2009].  In [@DeraemaekerNasser:2010], 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. [ @ShindoEtal:2011] 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 $4 E_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 $\epsilon^T_{33}$ constants. 

A sketch of the unit cell as used by [@DeraemaekerNasser:2010] is provided in Figure \[fig:MFC\_sketch\].  Conforming with [@ShindoEtal:2011], we use $W = \SI{413.90}{\micro\meter}$ for the width of an epoxy/PZT fibre set.  The width of the PZT fibre $w_p$ is then defined by the volume ratio $\rho$ of the fibre.  We use $\rho = 0.8$, which implies $w_p = \rho W = \SI{331.12}{\micro\meter}$. Height of fibre and electrode are also taken from [@ShindoEtal:2011], and set to $h_p = \SI{206.40}{\micro\meter}$ and $h_e = \SI{18}{\micro\meter}$. To compare our results to the values provided by [@DeraemaekerNasser:2010], 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 = \SI{206.40}{\micro\meter}$ and $L = \SI{1238.40}{\micro\meter}$.

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 \[fig:MFC\_bc\].  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 $\varphi = \pm \Delta 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 $\varphi = \pm\frac12 \Delta V$ is applied via a surface integral in the external work statement .  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 ${\mathbf{D}}\cdot{\mathbf{n}}$ are set to zero. 

The choice of material parameters was conducted in the following way: All constants that could be re-used from the linear model by [@DeraemaekerEtal:2009; @DeraemaekerNasser:2010] were used.  This includes the linear properties of epoxy ($Y_E = \SI{2.9}{\giga\pascal}$, $\nu = 0.3$, $\epsilon=4.25 \epsilon_0$) and copper ($Y_E = \SI{117}{\giga\pascal}$, $\nu = 0.31$).  In this example, we assume that the stiffness at constant electric field ${\mathbf{c}}^{E}$ is constant and isotropic, and defined by Young’s modulus $Y_E$ and Poisson ratio $\nu$. The permittivity at constant strain ${\boldsymbol{\epsilon}}^S$ is assumed isotropic, but is not constant.  In linear theory, it is smaller than the permittivity at constant stress ${\boldsymbol{\epsilon}}^T$. In this contribution, it is assumed to be a convex combination of $\epsilon^T = 1900 \epsilon_0$ for zero polarization, and $900 \epsilon_0$ for fully polarized material.  Further, the additive part of the free energy function $\hat \psi^{i}$ is chosen such that its derivative needed in implementations reads $$\begin{aligned}
\hat \psi^{i}({\mathbf{P}^{i}}) &= \frac{h_0 P_0}{2 |{\mathbf{P}^{i}}|} \left(\log\left(\frac{P_0+|{\mathbf{P}^{i}}|}{P_0-|{\mathbf{P}^{i}}|}\right) \right) {\mathbf{P}^{i}}.\end{aligned}$$ The saturation polarization $P_0$ and hardening parameter $h_0$ were chosen such that the hysteresis measurements depicted in [@Hooker:1998 Figure 7(d)] could be reproduced.  The remanent saturation strain $S_0$ was chosen in accordance with the measurements provided by [@FanStollLynch:1999].  The regularization parameter for the dissipation function was chosen as $\varepsilon=10^{-4}P_0$, and to avoid convergence problems at the singularity at the electrode also $\hat \psi^{i}$ was regularized using the same parameter.  All material parameters are also provided in Table \[tab:materialPZT5A\]. 

  --------------------------------- ------------------------------------------
  Young’s modulus $Y_E$             $\SI{48e9}{\newton\per\meter\squared}$
  Poisson ratio $\nu$               0.41
  el.  permittivity $\epsilon^{T}$   $1900 \epsilon_0$
  piezoelectric coupling $d_{31}$   $\SI{-1.85e-10}{\meter\per\volt}$
  piezoelectric coupling $d_{33}$   $\SI{4.40e-10}{\meter\per\volt}$
  coercive electric field $E_0$     $\SI{1150e3}{\volt\per\meter}$
  saturation polarization $P_0$     $\SI{0.255}{\coulomb\per\meter\squared}$
  saturation strain $S_0$           $\SI{5.5e-3}{}$
  hardening parameter $h_0$         $1/(40\,000\epsilon_0)$
  --------------------------------- ------------------------------------------

  : Homogenization of macro fibre composite: material parameters for fibres made from PZT-5A.[]{data-label="tab:materialPZT5A"}

![ Homogenization of macro fibre composite: schematic sketch of unit cell.[]{data-label="fig:MFC_sketch"}](MFCd33_sketch "fig:"){width="0.98\columnwidth"}\

![ Homogenization of macro fibre composite: symmetry boundary conditions posed on eigth of unit cell, top: displacement boundary conditions, bottom: electric boundary conditions.[]{data-label="fig:MFC_bc"}](MFCd33_bc "fig:"){width="0.98\columnwidth"}\

![ Homogenization of macro fibre composite: remanent polarization at full poling electric field $3E_0$. []{data-label="fig:MFC_PI"}](MFC_112_PI_fullE "fig:"){width="0.98\columnwidth"}\

![ Homogenization of macro fibre composite: longitudinal strain $S_{zz}$ at full poling electric field $3E_0$. []{data-label="fig:MFC_strainzz"}](MFC_112_Strainzz_fullE "fig:"){width="0.98\columnwidth"}\

![ Homogenization of macro fibre composite: transverse strain $S_{yy}$ at full poling electric field $3E_0$. []{data-label="fig:MFC_strainyy"}](MFC_112_Strainyy_fullE "fig:"){width="0.98\columnwidth"}\

In Figure \[fig:MFC\_PI\] the absolute value at full poling electric field $3E_0$ is depicted.  Figures \[fig:MFC\_strainyy\] and \[fig:MFC\_strainzz\] 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 $\epsilon^{T}_{33}$ were computed by the following formulae provided by [@DeraemaekerNasser:2010], $$\begin{aligned}
d_{32} &= \frac{\bar S_{yy}}{\bar E_z}, & d_{33} &= \frac{\bar S_{zz}}{\bar E_z}, & \epsilon^{T}_{33} &= \frac{\bar D_z}{\bar E_z}.\end{aligned}$$ The average values of $\bar S_{yy}$ and $\bar S_{zz}$ were derived from the constant displacements arising at the boundaries, compare Figure \[fig:MFC\_bc\], $\bar S_{yy} = 2c_y / W$ and $\bar S_{zz} = 2c_z / L$. The average value of the electric field was chosen according to [@DeraemaekerNasser:2010] as $\bar E_z = \Delta 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 ${\mathbf{D}}\cdot {\mathbf{n}}$ over the symmetric boundary at $L/2$. The following values were found in our computations, $$\begin{aligned}
d_{32} & =  \SI{-1.1554e-10}{\meter\per\volt},\\
d_{33} &=  \SI{3.1598e-10}{\meter\per\volt},\\
\epsilon^T_{33} &=  \SI{1.7736e-08}{\ampere\second\per\volt\per\meter} = 2003 \epsilon_0 .\end{aligned}$$ We note that the piezoelectric constants are lower than those estimated by [@DeraemaekerNasser:2010].  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 [@DeraemaekerNasser:2010]. 

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 non-differentiable, 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. 

Acknowledgements
================

Martin Meindlhumer acknowledges support of Johannes Kepler University Linz, Linz Institute of Technology (LIT).\
This work has been supported by the Linz Center of Mechatronics (LCM) in the framework of the Austrian COMET-K2 program. 

[10]{}

A Arockiarajan, SM Sivakumar, and C Sansour.  A thermodynamically motivated model for ferroelectric ceramics with grain boundary effects. ,  19(1):015008, 2009. 

E. Bassiouny, A.F. Ghaleb, and G.A. Maugin.  Thermodynamical formulation for coupled electromechanical hysteresis effects–[I]{}. [ B]{}asic equations. ,  26(12):1279–1295, 1988. 

E. Bassiouny, A.F. Ghaleb, and G.A. Maugin.  Thermodynamical formulation for coupled electromechanical hysteresis effects–[I]{}[I]{}. [ P]{}oling of ceramics. ,  26(12):1297–1306, 1988. 

E. Bassiouny and G.A. Maugin.  Thermodynamical formulation for coupled electromechanical hysteresis effects–[I]{}[I]{}[I]{}. [ P]{}arameter identification. ,  27(8):975–987, 1989. 

E. Bassiouny and G.A. Maugin.  Thermodynamical formulation for coupled electromechanical hysteresis effects–[I]{}[V]{}. [ C]{}ombined electromechanical loading. ,  27(8):989–1000, 1989. 

D. Boffi, F. Brezzi, and M. Fortin. ,  volume 44 of [*Springer Series in Computational Mathematics*]{}.  Springer, Heidelberg, 2013. 

C.J. Bottero and M.I. Idiart.  An evaluation of a class of phenomenological theories of ferroelectricity in polycrystalline ceramics. ,  113(1):13–22, 2018. 

Franco Brezzi, Jim Douglas, Ricardo Dur[á]{}n, and Michel Fortin.  Mixed finite elements for second order elliptic problems in three variables. ,  51(2):237–250, 1987. 

Franco Brezzi, Jim Douglas, and L.D. Marini.  Recent results on mixed finite element methods for second order elliptic problems.  In Dorodnitsyn Balakrishanan and Lions, editors, [*Vistas in Applied Math.,  Numerical Analysis, Atmospheric Sciences, Immunology.*]{}  Optimization Software Publications, New York, 1986. 

W Chen and CS Lynch.  A micro-electro-mechanical model for polarization switching of ferroelectric materials. ,  46(15):5303–5311, 1998. 

Philippe G. Ciarlet. .  North-Holland Publishing Co., Amsterdam-New York-Oxford, 1978.  Studies in Mathematics and its Applications, Vol.  4. 

Arnaud Deraemaeker and Houssein Nasser.  Numerical evaluation of the equivalent properties of macro fiber composite (mfc) transducers using periodic homogenization. ,  47(24):3272–3285, 2010. 

Arnaud Deraemaeker, Houssein Nasser, Ayech Benjeddou, and Andr[é]{} Preumont.  Mixing rules for the piezoelectric properties of macro fiber composites. ,  20(12):1475–1482, 2009. 

Mourad Elhadrouz, Tarak Ben Zineb, and Etienne Patoor.  Constitutive law for ferroelastic and ferroelectric piezoceramics. ,  16(3):221–236, 2005. 

J Fan, WA Stoll, and CS Lynch.  Nonlinear constitutive behavior of soft and hard pzt: Experiments and modeling. ,  47(17):4415–4425, 1999. 

W. Han and B.D. Reddy. ,  volume 9.  Springer Science & Business Media, 1999. 

Thomas Hegewald, Barbara Kaltenbacher, Manfred Kaltenbacher, and Reinhard Lerch.  Efficient modeling of ferroelectric behavior for the analysis of piezoceramic actuators. ,  19(10):1117–1129, 2008. 

Matthew W Hooker.  Properties of pzt-based piezoelectric ceramics between-150 and 250 c. 1998. 

J.E. Huber and N.A. Fleck.  Multi-axial electrical switching of a ferroelectric: theory versus experiment. ,  49(4):785–811, 2001. 

JE Huber, NA Fleck, CM Landis, and RM McMeeking.  A constitutive model for ferroelectric polycrystals. ,  47(8):1663–1697, 1999. 

Declan C Hughes and John T Wen.  Preisach modeling and compensation for smart material hysteresis.  In [*Active materials and smart structures*]{}, volume 2427, pages 50–64.  International Society for Optics and Photonics, 1995. 

SC Hwang, CS Lynch, and RM McMeeking.  Ferroelectric/ferroelastic interactions and a polarization switching model. ,  43(5):2073–2084, 1995. 

Raja Jayendiran, M Ganapathi, and T Ben Zineb.  Finite element analysis of switching domains using ferroelectric and ferroelastic micromechanical model for single crystal piezoceramics. ,  42(9):11224–11238, 2016. 

Manfred Kaltenbacher, Barbara Kaltenbacher, Thomas Hegewald, and Reinhard Lerch.  Finite element formulation for ferroelectric hysteresis of piezoelectric materials. ,  21(8):773–785, 2010. 

M. Kamlah.  Ferroelectric and ferroelastic piezoceramics–modeling of electromechanical hysteresis phenomena. ,  13(4):219–268, 2001. 

Marc Kamlah and Ulrich B[ö]{}hle.  Finite element analysis of piezoceramic components taking into account ferroelectric hysteresis behavior. ,  38(4):605–633, 2001. 

S. Klinkel.  A phenomenological constitutive model for ferroelastic and ferroelectric hysteresis effects in ferroelectric ceramics. ,  43(22-23):7197–7222, 2006. 

C.M. Landis.  Fully coupled, multi-axial, symmetric constitutive laws for polycrystalline ferroelectric ceramics. ,  50(1):127–152, 2002. 

Stephan Lange and Andreas Ricoeur.  A condensed microelectromechanical approach for modeling tetragonal ferroelectrics. ,  54:100 – 110, 2015. 

L.V. Laxman, S. Maniprakash, and A. Arockiarajan.  A phenomenological model for nonlinear hysteresis and creep behaviour of ferroelectric materials. ,  229(9):3853–3867, Sep 2018. 

Christoph Lehrenfeld and Joachim Sch[ö]{}berl.  High order exactly divergence-free hybrid discontinuous galerkin methods for unsteady incompressible flows. ,  307:339 – 361, 2016. 

K Linnemann, S Klinkel, and W Wagner.  A constitutive model for magnetostrictive and piezoelectric materials. ,  46(5):1149–1166, 2009. 

SM Lobanov and AS Semenov.  Finite-element modeling of ferroelectric material behavior at morphotropic phase boundaries between tetragonal, rhombohedric and orthorhombic phases.  In [*Journal of Physics: Conference Series*]{}, volume 1236, page 012062.  IOP Publishing, 2019. 

R.M. McMeeking and C.M. Landis.  A phenomenological multi-axial constitutive law for switching in polycrystalline ferroelectric ceramics. ,  40(14):1553–1577, 2002. 

C Miehe, D Rosato, and B Kiefer.  Variational principles in dissipative electro-magneto-mechanics: A framework for the macro-modeling of functional materials. ,  86(10):1225–1276, 2011. 

Astrid S. Pechstein, Martin Meindlhumer, and Alexander Humer.  The polarization process of ferroelectric materials analyzed in the framework of variational inequalities, 2019. 

C. M. Sands and I. A. Guz.  Unidimensional model of polarisation changes in piezoelectric ceramics based on the principle of maximum entropy production. ,  78(1):249–259, Feb 2013. 

J[ö]{}rg Schr[ö]{}der and Holger Romanowski.  A thermodynamically consistent mesoscopic model for transversely isotropic ferroelectric ceramics in a coordinate-invariant setting. ,  74(11-12):863–877, 2005. 

A.S. Semenov, A.C. Liskowsky, and H. Balke.  Return mapping algorithms and consistent tangent operators in ferroelectroelasticity. ,  81(10):1298–1340, 2010. 

Yasuhide Shindo, Fumio Narita, Koji Sato, and Tomo Takeda.  Nonlinear electromechanical fields and localized polarization switching of piezoelectric macrofiber composites. ,  6(7):1089–1102, 2011. 

M. Sofonea and A. Matei.  History-dependent quasi-variational inequalities arising in contact mechanics. ,  22(5):471–491, 2011. 

S Stark, P Neumeister, and H Balke.  A hybrid phenomenological model for ferroelectroelastic ceramics.  part i: Single phased materials. ,  95:774–804, 2016. 

S Stark, P Neumeister, and H Balke.  A hybrid phenomenological model for ferroelectroelastic ceramics.  part ii: Morphotropic pzt ceramics. ,  95:805–826, 2016. 

S. Stark, P. Neumeister, and H. Balke.  Some aspects of macroscopic phenomenological material models for ferroelectroelastic ceramics. ,  80:359 – 367, 2016.