Nonlinear electromechanical coupling in ferroelectric materials: large deformation and hysteresis

Smart materials respond to external stimuli, e.g., electric fields, which enables their use as sensors and actuators. The electromechanical coupling of the direct and converse piezoelectric effects, for instance, is used for both actuation and sensing in diverse engineering applications. The response of ferroelectric materials depends on their state of remanent polarization and the presence of an external electric field. To extend the operational range of sensors and actuators, an accurate understanding of the evolution of the material’s state of polarization is imperative, which requires both physical and geometric nonlinearities to be taken into account. Moreover, polymeric smart materials like PVDF allow significantly larger deformation as compared to conventional piezoelectric ceramics. The electromechanical coupling in piezoelectric materials manifests in ferroelectric and ferroelastic hystereses, which are related to both reversible and irreversible processes. Focusing on the latter, we transfer phenomenological models for domain switching in ferroelectric materials to the geometrically nonlinear regime. For this purpose, we follow related concepts of geometrically nonlinear elastoplasticity, where the concept of a multiplicative decomposition of the deformation gradient plays a key role. Accordingly, an additional deformation path that describes the evolution of the poled state from the unpoled referential configuration is introduced. The constitutive response of the material to mechanical and electrical loads is discussed, and dissipative internal forces that drive the evolution of the remanent polarization are derived within a thermodynamical framework and the principle of maximum dissipation.

comparatively high coupling coefficients. These polycrystalline ceramics are characterized by a relatively high brittleness, which is why they do not admit large strains. For this reason, their constitutive models are usually formulated within the range of small (linearized) strains, which gives Voigt's linearized theory of piezoelectricity for the reversible response, see, e.g., Eringen and Maugin [8] and Yang [39]. The response of ferroelectric materials is intrinsically related to their state of remanent polarization, which is assumed to be constant for most applications. Electromechanical loadings beyond so-called coercive electric fields or coercive stresses, however, induce an irreversible re-orientation of dipoles in the material and therefore a change in the remanent polarization, see, e.g., the review article of Kamlah [11]. These mechanisms are referred to as switching and result in ferroelectric and ferroelastic hystereses. A first thermodynamically consistent macroscopic model for the irreversible response of ferroelectrics was presented by Bassiouny et al [1,2] and Bassiouny and Maugin [3,4]. The mathematical modeling and the computational approximation of switching phenomena in ferroelectric materials share many aspects with well-established approaches in elastoplasticity, see, e.g., Landis [13], McMeeking and Landis [18] and Klinkel [12]. Recently, novel mixed finite-element formulations for reversible and irreversible problems were proposed by Pechstein et al. [24][25][26] and Meindlhumer et al. [21].
Unlike ceramic materials, ferroelectric polymers, e.g., polyvinylidene (di)fluoride (PVDF), can be subjected to comparatively large strains well beyond the linear range. Based on the modeling of the reversible electromechanically coupled response presented in [10,33], we propose a first attempt toward the continuum modeling of the ferroelectric hysteresis in the geometrically nonlinear range of deformation. For this purpose, we repeatedly adopt the concept of a multiplicative decomposition of the deformation gradient, which has successfully been applied to the modeling of various inelastic processes, see Lubarda [15] for a review. The idea of multiple multiplicative decompositions was discussed by Meggyes [20]. In the present approach, the deformation gradient is decomposed into three parts: The irreversible part governs the evolution of a poled (intermediate) configuration from the unpoled reference configuration. The second part describes the reversible local response of a poled volume element subjected to an electric field. Under the electric field, it deforms into a second intermediate configuration, which can be regarded as mechanically destressed, i.e., the current configuration is relieved from mechanical stresses. The third part accounts for the elastic stress response of the material. The modeling of the dissipative response follows the concepts of finite-strain elastoplasticity [30][31][32], where ideas of damage modeling are included to incorporate changes in the material parameters due to poling [14]. The evolution equation for the dissipative internal forces follows as an associated flow rule from the constrained optimization problem that is obtained from the principle of maximum dissipation [32].
The present paper starts with a discussion of the fundamental relations of nonlinear electromechanics. We first consider the reversible case, for which we present both a spatial and a material description of the governing equations. Then, we discuss the case of irreversible constitutive processes resulting from a remanent polarization of the material. After a classical approach using the principle of maximum dissipation, we introduce the known concept of a dissipation function. Finally, the specific constitutive framework is introduced based on a multiplicative decomposition of the deformation gradient tensor into reversible and irreversible parts and on a corresponding additive decomposition of the free energy. A brief review of the numerical implementation and the presentation of results for two example problems complete the present paper.

Reversible response
Continuum mechanics is based on the idealization of material bodies as a more or less continuous distribution of matter. In problems of large (finite) deformations, we are used to a material description of the governing equations, in which we refer to a reference configuration B (0) of the body under consideration. The availability of a referential state allows us to relate physical quantities that represent volume and surface densities to given volumes and surfaces in the reference configuration. In electrostatics, on the contrary, the deformation of the domain of interest typically plays a tangential role, if any, and the domain's geometry is assumed to be known. Eventually, we intend to bring together continuum mechanics and electrostatics in strongly deforming bodies. For this reason, we choose to first adopt a spatial representation, in which we refer to a body's current configuration B, which is-figuratively and literally-where physics actually takes place.
To begin with, we recall the fundamental relations of reversible electromechanics, in which the state of remanent polarization is assumed constant, see, e.g., [8,39] for a more detailed exposition. In the quasistatic case, the spatial description of the local balance of linear momentum reads div σ + f + f E = 0. (1) The Cauchy stress tensor is denoted by σ , and f are mechanical body forces. A nonconducting dielectric body additionally experiences electrical body forces f E in the presence of a nonhomogeneous electric field e, which can be represented by means of the divergence of the electrostatic stress tensor σ E , In the above relation, the (macroscopic) polarization vector has been introduced as p. Within the spatial description of continuum mechanics, derivatives and differential operators are to be understood with respect to a set of appropriate spatial coordinates. Note that the electrostatic stresses are not uniquely defined since any divergence-free tensor can be added without changing the electrostatic body forces f E . By direct computation, we can convince ourselves that the common definition of σ E as the sum of the polarization stresses σ P , which, in turn, are defined by the tensorial product of the electric field and the polarization, and Maxwell's stress tensor σ M of classical electrostatics, satisfies the requirement (2): The vacuum permitivity is denoted by 0 . To obtain the above representation of the electrostatic stresses, we have already made use of the constitutive equation among the electrical fields, i.e., which relates the dielectric displacement d, the electric field, and the polarization vector. In a polarized dielectric, an electric field also induces an electric (body) couple c E , which enters the balance of angular momentum as an additional source term. By allowing body couples, we leave the field of classical continuum mechanics, in which all torques are assumed to be moments of forces [38]. Instead, we consider dielectric materials as a case of polar continua, a concept introduced in the pioneering work of the Cosserat brothers, see, e.g., Maugin [17]. Unlike Cauchy's classical concept, not only the position but also the orientation of the material points (i.e., the polarization) needs to be accounted for in polar continua. From the local form of the balance of angular momentum, we eventually find that the sum of (mechanical) Cauchy stresses and polarization stresses is symmetric, rather than the mechanical stresses alone. In what follows, we will primarily consider the total (Cauchy) stresses, which are defined as the sum of the mechanical stresses and the electrostatic stresses, In view of the symmetry condition (8) and the definition of Maxwell's stresses (4), we note that the total stress tensor must be symmetric, i.e., σ tot = (σ tot ) T . Turning our attention to the balance laws for the electrical quantities, we first recall Gauss' law of conservation of charge, whose local form reduces to div d = 0 ( 1 0 ) in the case of nonconducting continua, for which the free charge density vanishes. In quasistatic problems, Faraday's law of induction in local form, implies that the electric field is conservative, i.e., the electric field can be represented by the (negative) gradient of the scalar potential φ. The mechanical and the electrical domains are brought together in the balance of energy, whose local form, assuming isothermal conditions, is given by where the polarization per unit mass π has been introduced. The velocity field is denoted by v; ρ is the mass density per unit (current) volume. The spatial velocity gradient can be expressed in terms of the deformation gradient tensor F as grad v =Ḟ · F −1 . Material time derivatives are indicated by superimposed dots throughout the present paper. The specific internal energy e is a thermodynamic potential in the deformation gradient F and the polarization p as independent fields, i.e., e = e (F, p). In other words, the electric field plays the role of a generalized force that drives changes in the state of polarization. By a Legendre transformation, the electric field is introduced as independent variable in the potential ψ = ψ(F, e), We obtain the corresponding balance law from (12) as Upon the Legendre transformation (13), the polarization has become the generalized driving force for changes in the electric field. We want to mention here that the terminology for thermodynamic potentials is not consistent in the literature. Eringen and Maugin [8] and Yang [39] refer to ψ as (generalized) Helmholtz free energy, whereas Miehe et al. [22] denote a potential in the deformation gradient and the electric field as mixed energyenthalpy function. Miehe et al. [22], in turn, use the term free energy function in their energy-based formulation for a thermodynamical potential that is a function of the deformation gradient and the dielectric displacement. We will adopt the terms free energy or simply energy in what follows. The balance law (14) suggests the definition of the electromechanical power density p per unit current volume as In nondissipative, reversible processes, the electromechanical power represents the source term in the balance of free energy, which then simplifies to Comparing the balance law (14) with the generic rate of the free energyψ(F, e), we can identify the constitutive equations for the Cauchy stress tensor and the polarization vector as Not least from a computational point of view, we introduce a total stress formulation by augmenting the free energy ψ with the free-space (or vacuum) energy ψ aug , The material description of the total free energy ψ tot will subsequently serve as basis for the variational formulation of the governing equations used in the finite-element approximation. The rate of the (specific) augmented free energy follows from the above relation aṡ Substituting the fundamental relation (14), we obtain the balance law for the augmented free energy as or, in terms of the total stress tensor (9), The constitutive equations for the total stresses and dielectric displacement again follow from a comparison of the total power, and the material time derivative of ψ tot = ψ tot (F, e) as

Material description of the governing equations
In solid mechanics, we typically make use of the fact that a referential state of the body under consideration is readily available. Material points are identified by their position vector X in the reference configuration, whose components with respect to some spatially fixed Cartesian frame (e 1 , e 2 , e 3 ) are denoted by X = X 1 e 1 + X 2 e 2 + X 3 e 3 . Though the electric field is not intrinsically tied to matter, we introduce a material description, in which we follow the evolution of physical quantities along the path of a material point. For this purpose, we transform the global (integral) balance laws such that integration is performed in the reference configuration B (0) rather than in the current configuration B, see, e.g., [6,7]. In Gauss' law of conservation of charge, we use the relation among material surface elements in the reference configuration and the current configuration da = J F −T · da (0) with J = det F, which is also known as Nanson's formula, The above equation suggests the definition of the material dielectric displacement as The transformation of material line elements dx = F · dX allows us to rewrite the global form of Faraday's law of induction as We introduce the integrand of the line integral in the reference configuration as the material electric field Using the definitions of the material dielectric displacement (26) and the material electric field (28), the corresponding local balance laws retain the same structure as in the spatial representation, i.e., Gauss' law reads Div and the material electric field is related to the electric potential by Note that derivatives with respect to spatial coordinates are replaced by derivatives in material coordinates, which, in the present coordinate-free representation, is indicated by capitalized differential operators. For the material polarization vector, we assume the same transformation as for the dielectric displacement, In terms of their material counterparts, the relation among the electrical fields (6) follows as Before proceeding with mechanical quantities, we spend a few remarks on the material description of electrical fields. As opposed to the spatial (physical) electrical fields, the material counterparts are invariant with respect to superimposed rigid-body motions. Let Q denote an orthogonal second-order tensor, which rotates the material body from the current configuration along with the electric fields, i.e.,F = Q · F, and, hence, From definitions (26) and (28), we immediately find that the material fields remain preserved upon the rotation, Rotation-invariant quantities facilitate the formulation of constitutive equations that satisfy the principle of material frame indifference. We note that other choices for rotation-invariant electric fields have been used in the material modeling. McMeeking and Landis [19], for instance, used the rotation tensor R in the polar decomposition of the deformation gradient F = R · U to define a rotation-invariant polarization as Π = p · R.
Regarding mechanical quantities, we recall the relation between the mass densities in the reference configuration and in the current configuration, In finite-deformation problems, constitutive equations are preferably formulated in the right Cauchy-Green tensor C or Green's strain tensor E, which are related to the deformation gradient by The second Piola-Kirchhoff stress tensor S is the conjugate stress measure, which is obtained from a pullback transformation of the Cauchy stresses σ to the reference configuration: In terms of the conjugate pairs of mechanical and electrical variables, the electromechanical power density per unit volume in the reference configuration follows from (15) as where the material polarization stresses (3), have been introduced using the conventional pullback of Cauchy-type stress measures (37). In reversible, isothermal processes, the balance of specific free energy is given by where we have introduced superscripts to clearly indicate the change in the independent fields in the material description, ψ = ψ (0) (C, e (0) ). Using the definition of the material electric field (28) and the mass density in the reference configuration (35), the free-space contribution to the total free energy (18) is rewritten as Its rate can be expressed in terms of the material Maxwell stresses S M , which are related to their spatial counterpart by With the material representation of the polarization stresses (39), the local balance of the total free energy ψ (0) tot follows in terms of the total stresses S tot = S + S P + S M as As in the spatial description, the source terms on the right-hand side of the above equation represent the total electromechanical power, whose material representation is

Dissipative response
In our summary of the equations that govern the reversible response in nonlinear electromechanics, the state of remanent polarization p i has been assumed to be constant. Processes that change the state of remanent polarization are irreversible by definition, which is why we have to abandon the idea of a constant remanent polarization p i . Unlike reversible processes, in which a formulation in terms of potentials naturally complies with the second law of thermodynamics, the constraint on the dissipated energy requires our attention in irreversible processes. For this reason, we recall the Clausius-Duhem inequality, which states that the dissipation density (per unit referential volume) d (0) is always positive, i.e., In irreversible processes, the total free energy can no longer be a function of the right Cauchy-Green tensor C and the material electric field e (0) only as it is in the reversible case. Instead, we assume ψ (0) tot to also depend on the remanent material polarization p (0) i , which can be regarded as an internal variable, Substituting the total power (45) and the generic rate of the total free energy (47) into the Clausius-Duhem inequality (46), we can identify the constitutive equations for the total stresses S tot and the material dielectric displacement and also the material representation of the dissipation and dissipative driving forcesê (0) : The dissipative driving forcesê (0) are conjugate to the remanent polarization p (0) i and therefore have the physical dimension of an electric field. To further specify the irreversible response, we adopt the principle of maximum dissipation, which is a fundamental concept in rate-independent plasticity, see, e.g., Simo and Hughes [32]. We expect the dissipation to take a maximum for a givenṗ (0) i , i.e., the dissipative driving forceŝ e (0) are determined such that In the above relation, E denotes the admissible domain, which is defined by a switching function f , The switching function represents a threshold for the dissipative driving forces in analogy to the yield function in plasticity. As long as we stay in the interior of the admissible range, the state of remanent polarization p (0) i is preserved. On the boundary of the admissible range ∂E, which can be regarded as yield surface where f (ê (0) ) = 0 holds, the remanent polarization may change depending on the direction of loading. The constrained maximization problem given by (51) can be solved by introducing a Lagrange multiplier λ, i.e., As shown by Simo and Hughes [32], the principle of maximum dissipation implies an associated flow rule that governs the (nonsmooth) evolution of the remanent polarization aṡ and the classical Karush-Kuhn-Tucker conditions for loading and unloading: From the algorithmic point of view, return-mapping methods are often used to construct (approximate) solutions, see, e.g., Simo and Hughes [32] for problems in plasticity and the formulations by Semenov [29] and [12] for the polarization of ferroelectric materials.
In the present formulation, we abandon the concept of Lagrange multipliers in favor of the notion of a dissipation function ϕ. Following the reasoning of Miehe et al. [22], we define the dissipation function such that its derivative with respect to the rate of the remanent polarization gives the driving forces: As the dissipation is generally nonsmooth, the derivative is to be understood as the subgradient of ϕ, which represents a set. In what follows, we will regularize the dissipation function appropriately such that the subgradient reduces to the conventional derivative, see [9,24]. To a priori fulfill the Clausius-Duhem inequality of positive dissipation (46), which, in terms of the dissipation function, transforms into we require ϕ to be normalized and positive, and convex, i.e., for all α ∈ [0, 1]. We further assume ϕ to be a positively homogeneous function of degree 1, which characterizes a dissipative response that is rate independent. Taking the derivative of the above relation with respect to α, we recover the dissipation (50), In other words, the dissipation density d (0) and the dissipation function ϕ are identical in the rate-independent case. As an example, we consider the classical switching condition that bounds the norm of the driving force by a threshold, which is referred to as coercive field strength e c in the context of ferroelectrics: By direct computation, we can convince ourselves that the dissipation function satisfies the principle of maximum dissipation (51). Its derivative with respect to the rate of the remanent polarization gives a rate-independent driving force (forṗ i = 0) in the direction ofṗ i , the norm of the driving force coincides with the coercive field ê (0) = e c . The above result for the driving force maximizes the dissipation (51) within the admissible range E, since we obtain For the numerical analysis, we regularize the dissipation function such that its derivative always exists. In the context of plasticity, Han and Reddy [9] proposed and analyzed a regularization of the form which they showed to result in errors of the order √ ε.

Multiplicative decomposition
The modeling of the constitutive response, in particular in coupled problems of multiple physical domains, undoubtedly is a challenging task and, quite commonly, subject of debates. Only few nonlinear problems allow us to motivate constitutive equations on the macroscopic continuum scale from microscopic considerations. Exemplarily, we mention material models of rubber, which can be derived from models of polymer chains [37]. Alternatively, if no assumptions on the structure of material models can be a priori made, the governing thermodynamic potentials are expanded into power series in the independent physical fields. To account for symmetries of a material, formulations in terms of invariants of the independent fields are often adopted in the modeling. The series coefficients represent the material parameters and need to be identified from experimental data.
The key idea of the multiplicative decomposition is to decouple the constitutive response in finitedeformation problems that are characterized by multiple underlying processes, e.g., thermoelasticity or elastoplasticity [15]. By assuming appropriate models for the decoupled subproblems, we endow constitutive models for the coupled response with a suitable structure, by which the description of the coupling is ideally simplified.
From a phenomenological point of view, we observe that ferroelectric materials deform in the course of the poling process, in which the state of remanent polarization is irreversibly altered. To understand the underlying mechanism, recall that ferroelectric materials are characterized by the presence of microscopic structures (unit cells) that posses a spontaneous polarization, i.e., the centers of positive and negative charges do not coincide and an electric dipole is formed, see, e.g., Kamlah [11] for a detailed exposition on ceramics. The unit cells may form so-called domains, which represent regions of equal spontaneous polarization. Stepping up from the microscale, physical properties on the macroscopic level can be regarded as averages of microscopic quantities, see Eringen and Maugin [8]. The remanent polarization p i represents the macroscopic average of the spontaneous polarization of unit cells and domains, respectively. In the unpoled state of a ferroelectric material, the polarization of microstructures is randomly oriented, which is why p i vanishes. Under sufficiently large electric fields, domains of ferroelectric materials switch their orientation, i.e., the spontaneous polarization becomes aligned with the applied field and maintains that orientation upon unloading. The re-orientation of the spontaneous polarization shifts the centers of positive and negative charges. The microscopic deformation of unit cells induced by the poling process macroscopically manifests as a remanent state of strain. In view of coupling between remanent polarization and deformation, mechanical loadings can also induce domains to irreversibly switch their state in ferroelectric materials. Though the underlying physical mechanism is the same as for electric loadings, the term ferroelasticity is often used to describe this characteristic property of ferroelectrics.
From a mechanical perspective, the poling process bears some resemblance to phenomenological models of elastoplasticity, in which irreversible deformation may occur if mechanical loads reach (or exceed) a threshold value, i.e., the yield stress. The similarities between poling and plasticity are reflected in both the structure of the constitutive models and numerical methods that have been proposed for the construction of approximate solutions. Leaving the range of infinitesimally small deformations, we also take the fundamental ideas of nonlinear elastoplasticity as inspiration in what follows. If deformations are large, the conventional additive decomposition of stresses and strain measures into reversible and irreversible parts loses some of its physical underpinning. Instead, we decompose the deformation gradient F into a product of an irreversible part F i and a reversible part F r : By the multiplicative decomposition, the notion of an additional intermediate configuration B (1) is introduced in the deformation path. In the present context of ferroelectric materials, the irreversible part F i locally describes the evolution from the unpoled reference configuration into a poled intermediate configuration, which depends on the state of remanent polarization. The second part F r , in turn, represents the reversible deformation induced by mechanical or electric loadings. Owing to the electromechanical coupling in ferroelectric materials, the reversible part F r can again be split into two parts, i.e., an inelastic (electrical) part F r,e and an elastic (mechanical) part F r,m , The electrical part F r,e governs the local response of a poled material due to the presence of an electric field, whereas the current state of stress is determined from the mechanical part F r,m . In the present approach, the total deformation gradient F is therefore composed of three parts, which implies that a second intermediate configuration B (2) is introduced into the deformation path, see Fig. 1.
For later use, we additionally define the electrical part of the deformation gradient as We can associate deformation and strain measures with each part of the deformation gradient. The right Cauchy-Green tensor related to the irreversible electrical part F i is naturally defined by We will later also make use of the right Cauchy-Green tensor of the mechanical part of the deformation, where the relation F r,m = F · F −1 e has been used. The above representation in terms of C also reveals an intuitive interpretation of C m : According to its definition, C m represents a material deformation measure in the reversible mechanical deformation path, but, at the same time, it can be regarded as a spatial deformation measure in the electrical parts of the deformation, which follows from a push-forward of C by means of F e .  Fig. 1 The multiplicative decomposition of the deformation gradient splits the deformation path into three parts: The irreversible part governs the evolution of a poled but unloaded first intermediate configuration B (1) from the unpoled reference configuration B (0) . The reversible electrical part describes the local deformation induced by reversible electromechanical coupling, which results in a second, mechanically unloaded intermediate configuration B (2) The decomposition of the total deformation into three parts is also reflected in the (total) free energy, which is additively composed of three contributions that correspond to the respective deformation paths: i.e., we have an irreversible part denoted by ψ i , a reversible electrical part ψ r,e and a reversible mechanical part ψ r,m . The fundamental idea of the multiplicative decomposition is to separately regard each of the three (sub)paths of the total deformation and to formulate appropriate relations for the material response.

Irreversible response
To further specify the proposed constitutive model, we first consider the irreversible part in the deformation. From the kinematic point of view, we note two characteristic features that are to be captured by the present approach: First of all, the irreversible deformation induced by the remanent polarization is anisotropic. Polycrystalline ferroelectrics usually have one principal direction that coincides with the remanent polarization, whereas polymeric materials may also be orthotropic due to the stretching that is applied during poling. For the case of transversally isotropic materials, the anisotropy can be captured by a second-order structural tensor, which is determined by the material remanent polarization in the reference configuration, Secondly, the irreversible deformation is volume-preserving, i.e., the condition J i = det F i = 1 must hold. Unlike plasticity, where plastic deformation accumulates until fracture, the number of switchable unit cells or domains of a ferroelectric material is finite. Therefore, the remanent polarization in ferroelectrics eventually saturates well before electric breakdown occurs, and so does the irreversible deformation. To describe the degree of saturation, we introduce the ratio between the material remanent polarization p (0) i and the saturation polarization p sat : In view of the isochoric nature of the poling-induced deformation, we assume a stretching-type deformation with a saturation strain S sat in analogy to what Zäh and Miehe [40] proposed for the electrostrictive response of graft polymers, The parameter γ is proportional to the saturation strain S sat and the square of the degree of polarization by As volumes are preserved, the mass density in the reference configuration is also preserved under the irreversible deformation, i.e., ρ (0) = ρ (1) . For the irreversible part of the free energy, we translate the phenomenological model used by Landis [13] and McMeeking and Landis [18] to nonlinear problems: If the remanent polarization saturates, i.e., β → 1, the energy tends to infinity, which prohibits any remanent polarization beyond p sat . The parameter c in the above relation represents a kinematic hardening modulus that has the physical dimension of inverse permittivity.

Reversible electrical response
As a next step, we consider the reversible response of material volume elements to electrical loads, i.e., an (externally) applied electric field. For the corresponding deformation path, the second intermediate configuration can be viewed as target configuration, whereas the first intermediate configuration represents the referential state. In the constitutive relations, we therefore use the electric field e (1) , which is related to its spatial and material counterparts by To describe the electromechanical coupling generically, we could assume a series expansion of the electrical part of the deformation gradient F r,e in the electric field e (1) as F r,e = I + e (1) · d (1) + e (1) e (1) : where the third-order tensor d (1) describes the piezoelectric coupling and b (1) is the fourth-order tensor of electrostrictive moduli. As the state of remanent polarization may change, however, the material parameters can be no longer assumed to be constant in general. We recall that the piezoelectric effect is only present in poled ferroelectrics, and also electrostriction is influenced by the remanent polarization. In what follows, our focus lies on a conventional piezoelectric coupling, which is why we neglect second-order electrostrictive effects. The piezoelectric coupling is kinematically described by three Biot-type strain tensors, The individual parts of the piezoelectric coupling, which are assumed to scale linearly in the degree of polarization β, correspond to the d 33 , d 31 , and d 15 -effects of Voigt's linear theory of piezoelectricity. In terms of the electric field e (1) and the normalized remanent polarization n (1) i in the first intermediate configuration, the three parts of the deformation are given by (1) , (1) , where (·) s denotes the symmetric part of a second-order tensor. We have also introduced the structural tensor in the first intermediate configuration B (1) as The structural tensor is pushed forward from and pulled back to the reference configuration, respectively, by means of the transformations where i is the square of the irreversible stretch in the direction of the remanent polarization. Within the framework of the multiplicative decomposition, the couplings among the individual processes of the deformation path are described kinematically. The (reversible) electromechanical coupling has already been taken care of in terms of F r,e . Having specified the constitutive response from the kinematic point of view, the thermodynamic relations remain to be provided in terms of an appropriate potential ψ (0) r,e . For this purpose, we disregard any coupling and assume the material to respond to a loading by an electric field as a (deformable) dielectric. Accordingly, the reversible electrical part of the internal energy is quadratic in the reversible polarization, where denotes the second-order permittivity tensor. The permittivity characterizes the capability of a material to (reversibly) polarize under an applied electric field. By means of a Legendre transformation, we replace the dielectric displacement d by the electric field e as independent field. Together with the constitutive relation of dielectrics d = · e + p i , we obtain the spatial description of the specific free energy as The permittivity of ferroelectric materials generally depends on both the current deformation and the state of remanent polarization, and, unlike dielectric elastomers [33,34], is not necessarily isotropic, see, e.g., Kamlah [11]. As we are interested in the behavior of the material (rather than the vacuum part), we introduce the second-order susceptibility tensor χ as In conventional polycrystalline piezoelectric materials, a transversally isotropic behavior is observed. The anisotropy is induced by the remanent polarization, which introduces a preferred direction in the material. For this reason, we assume the susceptibility to be composed of an isotropic part and an anisotropic part Δχ that results from the remanent polarization, In view of what is known from deformable dielectrics [19,33], the isotropic part has been assumed as χ(ρ) = J −1 χ 0 , where χ 0 is a material constant. From the spatial representation (86), the specific free energy corresponding to the reversible electrical deformation path follows upon substitution of the material quantities as a sum of an isotropic part and an anisotropic part, which are given by r,e,aniso = − In the isotropic part, we have recovered the free-space contribution of the free energy (41). The anisotropic part of the susceptibility has been pulled back to the reference configuration by means of the transformation which naturally follows from the material representation of the free energy. Exemplarily, we assume a transversal isotropy that depends linearly on the degree of remanent polarization β, i.e., where χ and χ ⊥ are material constants that represent the maximum change in susceptibility in the poling direction and in the direction perpendicular to it. Note that the anisotropic part of the susceptibility only depends on the (material) remanent polarization but not on the current state of deformation.

Reversible mechanical response
The elastic response is governed by the reversible mechanical part of the free energy ψ r,m . We can choose an appropriate strain energy function from finite-strain elasticity, which is typically formulated in terms of Cauchy-Green or Green-type measures of deformation. In view of the multiplicative decomposition, the strain energy refers to the deformation measures of the second intermediate configuration, which represents the referential state for the reversible mechanical path. Assuming a transversally isotropic behavior, the elastic strain energy ψ r,m depends on the structural tensor a (2) in the second intermediate configuration B (2) , i.e., a (2) = n (2) i n (2) i , In analogy to the structural tensor in the first intermediate configuration (84), its counterpart in the second intermediate configuration a (2) follows upon a push-forward from the reference configuration by where C e,ii = n i denotes the square of the electrically induced stretch in the direction of the remanent polarization.
As mechanical deformation measures C m or E m = (C m − I)/2 are used in the strain energy function, one is tempted to object that ferroelectric and piezoelectric couplings are missing. Analogously to the reversible electric path, we assume an uncoupled elastic stress response since any coupling is represented kinematically in terms of F i and F r,e , respectively. Exemplarily, we mention the St. Venant-Kirchhoff material of finitestrain elasticity, which can be regarded as a generalization of Hooke's law of linear elasticity to geometrically nonlinear problems, ρ (2) ψ (2) r,m = where c denotes the fourth-order tensor of elastic moduli. In our example problems, we will use a neo-Hookean material model, see Sect. 5.

Constitutive equations
Having specified the kinematics of the nonmechanical deformation and the thermodynamic potentials of each deformation path, we evaluate the general constitutive equations (49) and (50) in the present multiplicative setting and discuss the implications in what follows. For the total stresses, the corresponding second Piola-Kirchhoff stress tensor follows from the derivative of the respective free energy functions as If the nonisotropic part of the material susceptibility does not depend on the deformation as in example (92), the above relation is evaluated as where we have used the chain rule of differential calculus and the derivative of C m with respect to C according to (71). The irreversible part of the free energy ψ (0) i does not depend on the state of deformation, which is why it does not enter the above relations. From the free-space energy, we obtain the (material) Maxwell stresses that contribute to the (total) stress response. The constitutive relation for the total stresses (97) suggests the definition of the elastic stresses in the second intermediate configuration as the derivative of the reversible mechanical part of the free energy, i.e., The corresponding stresses in the reference configuration are obtained upon a pullback from the second intermediate configuration by means of the electrical part of the deformation gradient F e , Note that the irreversible part ψ (0) i is not only independent of the deformation, but also independent of the material electric field e (0) , which is why the material dielectric displacement field follows from the derivatives of the reversible parts of the (total) free energy ψ ∂F r,e ∂ e (1) In view of the decomposition of the susceptibility into isotropic and anisotropic parts (88), we split the material polarization vector, which follows from the dielectric displacement by discarding the free-space contribution, into four separate parts, The three parts in addition to p (0) i are identified from the constitutive equation (100) and the susceptibility (88) as The first term p (0) elast depends on the elastic stresses and represents the electromechanical coupling, which is kinematically described by the reversible electrical part of the deformation gradient F r,e . The second term p (0) iso describes the isotropic dielectric response of the material under an electric field, and the polarization induced by the anisotropy of the susceptibility is denoted by p (0) aniso . The decomposition of the polarization allows us to gain more insight into the stress response (97). Upon substitution of the elastic stresses (99) along with the isotropic polarization (102), the constitutive equation for the total stresses follows as From the definition of the total stress tensor (9) and the polarization stresses (39), we recognize that the mechanical stresses are related to the elastic stresses S elast by The above result illustrates that the mechanical stresses are composed of a symmetric elastic part and the polarization stresses induced by the electromechanical coupling, an anisotropic electric susceptibility and the remanent polarization. We continue with the dissipative driving forces, which are obtained from the derivative of the total free energy with respect to the material remanent polarization p (0) i . The irreversible part ψ (0) i naturally depends on the remanent polarization. We have assumed the isotropic part of the susceptibility (88) to describe the dielectric response of an unpoled material. Therefore, only the anisotropic part of reversible electric free energy (90) enters the constitutive equation of the dissipative driving force. As the material parameters that govern the elastic response of the material also depend on the state of remanent polarization in general, ψ (0) r,m also contributes to the driving forces, which eventually follow aŝ Within the framework of the multiplicative decomposition, the concept of destressing is introduced to illustrate the notion of intermediate configurations, see Fig. 1. The second intermediate configuration B (2) is taken by a body if no mechanical deformation occurs, i.e., F r,m = I and C m = I hold. In this case, the symmetric Piola-Kirchhoff stresses in the second intermediate configuration (98) vanish identically, S (2) elast = 0. In other words, the second intermediate configuration follows upon elastic unloading of the current configuration. The total stresses are only composed of the Maxwell stresses and the polarization stresses induced by the isotropic dielectric response. As the electromechanical coupling polarization also vanishes, p (0) elast = 0, the mechanical stresses in the second intermediate configuration are the negative of the sum of the polarization stresses that correspond to the dielectric anisotropy and the polarization stresses from the remanent polarization.
We recover the first intermediate configuration B (1) in the deformation path if the body is also unloaded electrically, i.e., e (0) = 0. As a consequence, both total and mechanical stresses vanish in this configuration, and only the remanent polarization remains from the electrical point of view. Accordingly, we can regard the first intermediate configuration as a fully unloaded but poled configuration. As the individual parts of the deformation gradient describe the local response of material volume elements in the respective configurations, the intermediate configurations are incompatible as soon as inhomogeneities occur due to external loadings or kinematic boundary conditions. In fact, incompatibilities are the sources for a nontrivial elastic stress response S elast and the corresponding polarization p (0) elast .

Algorithmic treatment
In view of both geometric and physical nonlinearities considered within the proposed approach, we need to make use of some numerical method to construct approximate solutions. In what follows, we develop the incremental variational formulation of the governing equations, which provides the basis for an algorithmic treatment by means of the finite-element method. For this purpose, we first introduce the electromechanical power P as the integral of the material power density p where the rate of the (total) free energy and the dissipation are integrals of the respective densities (44) and (50): The time integral of the total electromechanical power over an interval Δt = t n+1 − t n gives the corresponding work performed in that period. Subtracting the power of the external forces P ext , for which we assume a potential W ext to exist, we obtain the incremental work functional The work of the external loadings is composed of a mechanical part W m,ext and an electrical part W e,ext , The mechanical part, in turn, comprises body forces f (0) and surface tractions t (0) , which are to be understood as total stresses. The work of electrical loads depends on the density of free charges σ (0) e per unit surface in the reference configuration B (0) , which represents the outer normal component of the dielectric displacement, i.e., σ (0) e = d (0) · n (0) , In the time-discrete setting, we approximate the energy dissipated during the time period Δt by By requiring the variation in the independent fields of the work functional to vanish, i.e., δΠ = 0, we obtain the variational formulation of the governing equations, which generally serves as basis for a finiteelement analysis. We use the open-source, multipurpose finite-element code Netgen/NGSolve (www.ngsolve. org), which relieves us from the burden of computing the variations in the independent fields and a consistent linearization by hand. Owing to automatic differentiation, we only have to specify the work functional (108) symbolically. We rely on Newton's method to solve the nonlinear system of equations iteratively for the increments in the unknown fields, which are the displacement field u, the electric potential φ, and the material remanent polarization p (0) i . Due to the introduction of the dissipation function, we can solve for all fields at once without resorting to some kind of return-mapping algorithm. In our example problems, we use hexahedral elements to spatially discretize the bodies. Second-order hierarchical polynomials are used as shape functions for the displacement field and for the electric potential. The components of the remanent polarization are interpolated by means of linear polynomials.

Uniaxial poling
By means of the first example, we want to check the consistency of the proposed nonlinear formulation with results obtained from geometrically linear models. For this purpose, we study the uniaxial poling process of a ferroelectric material under a homogeneous electric field. In all our examples, we use material properties that are derived from the dimensionless parameters used by McMeeking and Landis [18], see Table 1. To further connect to Voigt's linear theory of piezoelectricity, we recall that the second intermediate configuration B (2) , which evolves from the first intermediate configuration B (1) under electric loads, can be regarded as mechanically unloaded. In the mechanical part of the deformation path, which leads to the current configuration B, the electric field is assumed to be constant. Therefore, the material parameters of the present approach correspond to the (tangent) moduli of the linear theory, where the electric field and the mechanical stresses are used as independent variables. The elastic behavior is assumed to be isotropic and independent of the remanent polarization, and no anisotropic parts develop in the electric susceptibility. In our examples, we use the neo-Hookean model for compressible materials, for which the strain energy is given by ρ (2) ψ (2) r,m =    see, e.g., [5]. In the above relation, λ denotes the first Lamé parameter and μ is the shear modulus. We consider a cube-shaped specimen with 1 mm side length, for which all rigid-body degrees of freedom are constrained. The edges of the cube are aligned with the coordinate axes (x 1 , x 2 , x 3 ) of the global frame. The bottom surface parallel to the (x 1 , x 2 )-plane is grounded; on the top surface, an electric potential V 0 is applied such that a homogeneous electric field in x 3 -direction is induced. By a cyclic electric loading with a maximum voltage amplitude of V 0 = ∓ 1700 V in steps of 10 V, we obtain the dielectric hysteresis illustrated in Fig. 2a.
The corresponding butterfly hysteresis, which describes the dependency of the strain E 33 in the poling direction as a function of the applied electric field e 3 , is shown in Fig. 2b. Additionally, we also provide the curves for the reversible electrical strain due to the piezoelectric coupling E r,e,33 and the poling-induced irreversible strain E i,33 . Our geometrically nonlinear approach based on the multiplicative decomposition is able to reproduce what is known from the geometrically linear case, cf. McMeeking and Landis [18].

Poling of a bent cantilever
As geometric nonlinearities are negligible in the previous benchmark example, we want to study a situation in which large deformations need to be considered. For this purpose, we investigate the poling of a ferroelectric cantilever, which is first deflected by a tip force in transverse direction, before it is poled by applying an electric voltage V 0 between its top and bottom surfaces. The cantilever, which has the dimensions L × W × H = 80 mm × 20 mm × 0.4 mm, is aligned with the global x 1 -axis in its reference configuration; x 3 is the thickness direction. The material parameters are taken from Table 1. The cantilever is discretized into hexahedral elements using six equidistant layers over the thickness. The in-plane discretization is shown in Fig. 4. For the sake of computational efficiency, symmetry conditions at the plane X 2 = 0 are exploited in the numerical analysis. The transverse tip force of F = −0.8 N is realized as a uniform surface traction acting on the front face Fig. 3 Deformed configuration of a ferroelectric cantilever subjected to a transverse tip load: the top surface is grounded; the voltage V 0 at the bottom surface is gradually increased in the poling process, which is initiated from the deformed equilibrium configuration. The magnitude of the remanent polarization p (0) i is illustrated at a poling voltage of V 0 = 420 V (X 1 = L) of the cantilever in the reference configuration. At the clamped end, any displacement is prohibited, i.e., u(X 1 = 0, X 2 , X 3 ) = 0. Under the imposed load, the tip deflection of the center plane is obtained as u 3 (X 1 = L , X 2 = 0, X 3 = H/2) = −32.84 mm. We assume the top surface to be grounded; the voltage V 0 at the bottom surface is gradually increased in a second loading step. The mechanically deformed, poled configuration (V 0 = 420 V) of the cantilever is illustrated in Fig. 3, which also shows the magnitude of the (material) remanent polarization p (0) i . In a homogeneous poling problem as our first example, we would expect the material to polarize starting from a voltage of V 0 = e c H = 400 V for the given dimensions. The present problem of a bent cantilever, however, is inhomogeneous due to the mechanical load and the imposed boundary conditions, which result in an inhomogeneous state of deformation. We observe remanent polarization to occur well before the threshold of the homogeneous problem.  (Fig. 4a), the top surface (X 3 = H ) does not exhibit any significant remanent polarization. The bottom surface (X 3 = 0), by contrast, has already been poled irreversibly to some extent. The remanent polarization starts from the clamped end of the cantilever (X 1 = 0). Directly at the clamping, however, poling is obstructed by the kinematic boundary conditions on the displacement, which constrain the deformation induced by the poling process. At V 0 = 400 V, we find regions with the highest remanent polarization in the center of both the top and the bottom surfaces of the cantilever, see Fig. 4b. Toward its free edges, no significant remanent polarization has yet occurred. When we increase the applied voltage to V 0 = 410 V (Fig. 4c), the polarization increases substantially and extends further toward the front (X 1 = 0) and lateral (X 2 = ±W/2) faces of the cantilever. At the lateral faces of the cantilever, the remanent polarization diminishes toward the clamped end.
To gain additional insight, we investigate the evolution of the (material) remanent polarization p i,3 within three cross sections along the cantilever's axis, i.e., cross sections close to the clamped end (X 1 = 0.1L), in the middle (X 1 = L/2), and close to the free end (X 1 = 0.9L). Figure 5a-c shows the remanent polarization as a function of the applied voltage V 0 at the (vertical) plane of symmetry (X 2 = 0) and three (horizontal) planes parallel to the transducer, i.e., the bottom surface (X 3 = 0), the center plane (X 3 = H/2) and the top surface (X 3 = H ). Apparently, a notable remanent polarization first occurs at the bottom surface (X 3 = 0) close to the clamped end (X 1 = 0.1L) of the ferroelectric cantilever at an applied voltage of approximately V 0 = 370 V, see Fig. 5a. At the center plane (X 3 = H/2), an increase in the polarization is observed at about V 0 = 392.5 V, whereas no considerable polarization occurs at the top surface (X 3 = H ) before a voltage of V 0 = 400 V is reached. At V 0 = 400 V, the rate of change in the remanent polarization (in terms of applied voltage) increases significantly in all three planes of the cantilever. In the cross section at X 1 = L/2, we observe a similar behavior: The remanent polarization starts to grow from the bottom surface (X 3 = 0) toward . At V 0 = 410 V, the remanent polarization has already extended throughout both surfaces (c). The gridlines correspond to the spatial discretization of the finite-element mesh. Note that only half of the structure is illustrated in view of the symmetry at X 2 = 0 the top surface (X 3 = H ). As compared to the cross section close to the clamped end, the onset of poling is delayed to higher voltages in the bottom plane and the center plane of the cantilever. A substantial increase in the rate of polarization, however, already occurs at a lower voltage than at the clamped end throughout all three planes over the height of the cross section. Moving further toward the free end (X 1 = 0.9L), the poling process occurs at even higher voltages. The remanent polarization starts to grow almost simultaneously in all three planes of the cross section at approximately V 0 = 400 V. Figure 6 shows the transverse tip displacement u 3 (X 1 = L , X 2 = 0, X 3 = H/2) in the course of poling. First, the deflection increases slightly with the onset of poling at the bottom parts of the cantilever. The polarization induces tensile strains along the poling axis, i.e., in the thickness direction, and compressive strains in the plane perpendicular to it, see (75)-(82). Increasing the applied voltage, the tip deflection of the cantilever reduces as the poled region grows. i,3 within the plane of symmetry X 2 = 0. Values for cross sections at X 1 = 0.1L, X 1 = L/2 and X 1 = 0.9L and three planes at the bottom (X 3 = 0), center (X 3 = H/2), and top surface (X 3 = H ) of the cantilever are shown Fig. 6 Evolution of the transverse tip deflection u 3 (X 1 = L , X 2 = 0, X 3 = H/2) during the poling process: The deflection of the cantilever first increases slightly before it decreases as the remanent polarization and the poled regions grow From a phenomenological point of view, the observed behavior is consistent with what is referred to as ferroelasticity, i.e., a change of remanent polarization (and deformation) under mechanical loadings. Polycrystalline ferroelectrics can be depolarized under sufficiently strong compressive loads on the one hand, and under tensile loads perpendicular to the remanent polarization on the other hand [11]. Conversely, poling is supported if a material is compressed by a load perpendicular to the poling axis, which corresponds to what we observe in the present results. In the static equilibrium configuration, from which the poling process is started, the bottom surface of the cantilever is in compression due to the bending moment about the x 2 -axis, which is induced by the transverse tip load. The bending moment has a maximum at the clamped end, which is where remanent polarization occurs first, and it decreases toward the free end of the cantilever. Therefore, the onset of poling is delayed to higher voltages further along the axis of the cantilever. Over the height of the cross sections, the compressive (axial) normal stresses are largest at the bottom surface, which is also consistent with the results we obtained in our example. We note that ferroelasticity is not fully included in the proposed constitutive model. Irreversible deformations under mechanical loads in the macroscopically unpoled state, which are observed in polycrystalline ferroelectrics, cannot be captured by our approach, see, e.g., Landis [13] for a discussion.

Conclusions
In the present paper, we have proposed a constitutive framework for the phenomenological modeling of the irreversible response of ferroelectric materials in the geometrically nonlinear domain. For this purpose, we have used the concept of the multiplicative decomposition of the deformation gradient. The deformation gradient has been split into three parts that are related to the state of remanent polarization, the reversible electrical response, and the elastic stress response, respectively. We have adopted a material description of the governing equations, in which we have referred to an undeformed and unpoled reference configuration of the body under consideration. The constitutive equations for the (second Piola-Kirchhoff) stresses, the (material) dielectric displacement, and the dissipative driving forces follow from standard arguments of thermodynamics. In analogy to the theory of plasticity, the evolution equation for the remanent polarization has been obtained from the principle of maximum dissipation. Moreover, the concept of using a dissipation function to satisfy the principle of maximum dissipation has been introduced, which has significantly simplified the algorithmic treatment of the governing equations. Finally, the proposed physically and geometrically nonlinear formulation has first been validated against existing results for small deformations, and eventually results for an example problem with large deformations have been presented.