A complete direct approach to nonlinear modeling of dielectric elastomer plates

In this paper, we present a complete direct approach to nonlinear modeling of thin plates, which are made of incompressible dielectric elastomers. In particular, the dielectric elastomers are assumed to exhibit a neo-Hookean elastic behavior, and the effect of electrostatic forces is incorporated by the purely electrical contribution to the augmented Helmholtz free energy. Our approach does not involve any extraction-type procedure from the three-dimensional energy to derive the plate augmented free energy, but directly postulates the form of this energy for the structural plate problem treated in this paper. Results computed within the framework of this novel approach are compared to results available in the literature as well as to our own three-dimensional finite element solutions. A very good agreement is found.

to piezoelectric ceramics. For three-dimensional Eulerian and Lagrangian formulations, we refer to [20,21]. Given the typical applications, devices made of dielectric elastomers are mostly put into practice in the form of thin films, such that a structural mechanics approach is well suited, which motivates our modeling approach for dielectric elastomer plates as electro-elastic material surfaces with mechanical and electrical degrees of freedom. In general problems of dielectric elastomer actuators, numerical methods, such as the finite element method, are applied implementing solid elements for general three-dimensional problems [21][22][23][24] or solid shell elements to account for the typical thinness of the dielectric elastomer actuators, as developed in [25]. Within the three-dimensional modeling framework for dielectric elastomers, the works [26,27] are just a few examples of the exhaustive literature. Multiple extensions-e.g., to electro-viscoelasticity treated in [28] or electrostriction in [29]-have been reported as well.
Concerning the modeling of elastic plates and shells in a geometrically nonlinear framework, we refer the reader to the above listed papers by Eliseev as well as to [30][31][32]. Shell-type finite element models for electroactive polymers were investigated in [33], electro-elastic coupling of dielectric elastomers in combination with studies on diverse failure phenomena, e.g., pull-in instability and the formation of wrinkles, was addressed in [34], and further contributions to these topics can be found in [35,36]. With respect to our own contributions to the field, we only mention recent works related to modeling of electro-active plates and shells as electroelastic material surfaces. In [37] we have studied piezoelectric shells, and in [38] dielectric elastomer plates were investigated. Two-dimensional constitutive relations for the plate were derived by numerical integration of a three-dimensional augmented free energy through the plate thickness imposing a plane stress assumption and an a priori assumption concerning the distribution of the strain through the thickness of the plate in [38]; such an approach has also been used in [39] for hyperelastic shells.
The present paper is structured as follows. After a brief summary of the three-dimensional theory of dielectric elastomers, the concept of electro-elastic plates as material surfaces is developed in detail. We start with the material independent equations and focus on the direct derivation of a plate augmented free energy for the electro-elastic material surface afterward. This derivation is based on the polar decomposition of the surface deformation gradient tensor, which enables an additive decomposition of the augmented free energy into a membrane part, a bending part, and an electrical part. The latter is shown to account for the electrostatic forces responsible for the actuation behavior of the plate. A numerical implementation completes the theoretical part of the paper. Results computed within the framework of this novel approach are compared to results available in the literature as well as to our own three-dimensional finite element solutions. A very good agreement is found.

Three-dimensional formulation
This Section presents a brief introduction to the three-dimensional theory of electro-elasticity; for more detailed presentations we refer to [20,21,38]. We start with the energy balance in which we have the Helmholtz free energy per unit mass Ψ 3 = Ψ 3 (C 3 , E 3 ) with the right Cauchy-Green tensor C 3 = F T 3 ·F 3 , the material electric field vector E 3 , and the mass density ρ 0 in the reference configuration. S me 3 and P 3 are the-in general non-symmetric-mechanical second Piola-Kirchhoff stress tensor and the material polarization vector. The notions of material electric field vector and material polarization vector refer to the transformation of the corresponding spatial entities e 3 and p 3 by means of E 3 = F T 3 · e 3 and P 3 = J 3 F −1 3 · p 3 , see, e.g., [20] for details. Moreover, the deformation gradient tensors F 3 = (∇ 3 r 3 ) T and J 3 = detF 3 were introduced with the invariant differential operator ∇ 3 in the reference configuration and the position vector r 3 in the actual configuration. S pol 3 = P 3 E 3 · C −1 3 is the so-called second Piola-Kirchhoff polarization stress tensor, which is not symmetric in general, whereas the symmetric tensor S S 3 = S me 3 + S pol 3 is denoted as the symmetric second Piola-Kirchhoff stress tensor. Classical arguments from thermodynamics for reversible electro-elasticity render the constitutive relations as In nonlinear electro-elasticity, the Helmholtz free energy is augmented with an electrical term accounting for the contribution to the free energy stemming from the permittivity in vacuum,

Constitutive modeling for isotropic dielectric elastomers
For isotropic dielectric elastomers, we assume the Helmholtz free energy as the sum of a purely mechanical part depending only on C 3 and an electrical part, which results in Accordingly, the material electric polarization is computed solely from the electrical part of the Helmholtz free energy, which we introduce in analogy to the augmentation term as with a deformation dependent spatial susceptibilityχ(C 3 ). In the simplest setting of a material, which remains isotropic when strained, see [26], this dependency isχ(C 3 ) = χ J −1 3 with the constant material susceptibility χ . In this case, the material polarization is such that the polarization stress is symmetric, and it can be computed solely from the electrical part of the Helmholtz free energy. Combining the electrical part of the Helmholtz free energy with the augmentation part to a total electrical contribution to the augmented free energy ρ 0 Ψ el,tot with the permittivity ε = ε 0 (J 3 + χ), from which the electrostatic stress is computed as Finally, the mechanical second Piola-Kirchhoff stress tensor turns out to be symmetric, and it can be computed from the mechanical part of the free energy only. The latter part of the free energy is assumed as any isotropic hyperelastic strain energy function in the form in which I C 3 = trC 3 , I I C 3 = C 3 · · C 3 and I I I C 3 = detC 3 are the invariants of C 3 . In this paper, we will focus on incompressible isotropic dielectric elastomers with J 3 = 1. In this case, the total stress tensor and the material electric displacement vector are in which p is a Lagrange multiplier to account for the incompressibility condition. In particular, a simple incompressible neo-Hookean material law will be considered in the numerical examples, for which the augmented free energy with J 3 = 1 is with the permittivity ε = ε 0 (1 + χ) and the Lamé parameter μ as the only material parameters.

Electro-elastic plates as material surfaces
In this Section, we discuss the governing equations of a thin plate modeled as a two-dimensional electro-elastic material surface with mechanical and electrical degrees of freedom. For details concerning these equations, we refer the reader to [37,38]. In particular, we consider the plate as a two-dimensional continuum with five mechanical degrees of freedom, three translations and two rotations. This resembles the notion of a single rigid director d with d · d = 1 attached to each particle of the plate, which results in a Reissner-type theory; see [40]. Furthermore, we take the director to correspond to the unit normal vector of the material surface, d = n; then, we obtain a Kirchhoff-Love theory, see [41]. Concerning the electrical degrees of freedom, we use only the dominant one, i.e., the electric potential difference V .

Strain measures
The material surface is plane in its reference configuration, and it is denoted as reference surface. In its deformed or actual configuration, the material surface is denoted as actual surface. The first metric tensor of the plane reference surface A = I is the two-dimensional identity tensor, and the second metric tensor is zero, B = 0.
For the actual surface, the first and second metric tensors are a and b. The reference configuration and the actual configuration of the material surface are related to each other by means of a deformation gradient tensor F = (∇ 0 r) T with the differential operator ∇ 0 of the reference surface and the position vector r of points of the material surface in the actual configuration. Also note that the second metric tensor of the actual surface is b = b T = −∇n with the unit normal vector n of the actual surface, and a = ∇r = I 3 − nn with the three-dimensional identity tensor I 3 holds for the first metric tensor. The differential operators ∇ 0 and ∇ with ∇ 0 = F T · ∇ refer to either the reference surface or the actual surface. We introduce two tensor-valued Green strain measures for the material surface, which are defined as the difference between the two metric tensors in the two configurations, with the proper transformation by means of F, Both strain measures remain constant if and only if the motion of the material surface is a rigid body motion, see [7,42] for a discussion.

Material independent equations
In analogy to the three-dimensional formulation, a variational principle is introduced for the electro-elastic material surface as with the area A 0 of the reference surface. Note that the contribution from the external electric charge σ per unit area of the reference surface, σ δV with the voltage V , has migrated into the area integral. p and m are purely mechanical external forces and moments per unit area of the reference surface, and δ A i is the virtual work of the internal forces, moments, and charges. External forces and moments acting at the boundary, which involve both, mechanical and electrical contributions, are accounted for by means of the virtual work δ A e,b . Due to the fact that external forces and moments, p and m are assumed purely mechanical, ponderomotive forces must be accounted for by means of δ A i . As we have already discussed above, the two strain measures ε and κ remain constant, if and only if the motion of the material surface is a rigid body motion. Therefore, the virtual work of the internal forces, moments, and charges must vanish, if the virtual motion is a rigid body motion with the voltage fixed; hence, must hold. δr and δn are formally assumed to be independent, such that the constraint that the variation of the unit normal vector must lie in the tangential plane, δn + ∇δr · n = 0, is accounted for by means of a vector valued Lagrange multiplier q. Furthermore, two symmetric tensor-valued Lagrange multipliers τ and μ and a scalar-valued Lagrange multiplier q are introduced to account for the constraints δε = 0, δκ = 0, and δV = 0, such that the virtual work δ A i in the variational principle is formally replaced with −τ · · δε − μ · · δκ + qδV + q · (δn + ∇δr · n) .
Then, the variational principle is rewritten to Running through some straight forward mathematical procedures, we obtain the mechanical equilibrium conditions within A 0 as For the sake of brevity, we omit the specific form of the dynamic boundary conditions, as they will not be needed in this paper. The Lagrange multipliers τ and μ are identified as total second Piola-Kirchhoff stress measures, q as the transverse shear force vector and q as the internal charge per unit reference area. Moreover, total external and internal charges are introduced, such that holds in addition to the mechanical equilibrium conditions. Besides the derivation of the above equations, we can also conclude on the fact that the negative virtual work of the internal sources can be written as the variation of an energy function; in particular, we have with the augmented free energy Ω = Ω(ε, κ, V ) per unit reference area of the plate. Formally computing the variation of this energy as one finds the constitutive relations of the electro-elastic material surface as

Constitutive modeling for isotropic dielectric elastomer plates
The present approach to the modeling of plates as electro-elastic material surfaces involves two strain measures, ε and κ. Talking about hyperelastic plates, we rather work with the right Cauchy-Green tensor C than with ε; hence, we recall the definitions Next, we introduce the polar decomposition of the deformation gradient tensor as F = R · U, in which R is an orthogonal tensor with n = R · N and U = U T is the symmetric stretch tensor of the material surface, see, e.g., [41]; N is the constant unit normal vector of the plane reference surface. Then, the two strain measures may as well be written in the form In the case of a purely elastic problem, the principle of material frame indifference has been used in the literature (see [41]) to show that the strain energy W per unit reference area of the material surface must have the form Therefore, W = W (C, κ) satisfies this condition, and in analogy holds for the electro-elastic problem. For the purpose of constitutive modeling, we involve the alternative second strain measure which we denote as curvature tensor. Then, Ω = Ω(C, K, V ) is a proper form for the augmented free energy as well, for which we will discuss the specific form in the next Section.

Augmented free energy
The bottleneck in our formulation is the specific form of the augmented free energy Ω = Ω(C, K, V ); in particular, when it involves hyperelastic material laws. A straightforward approach is to involve the threedimensional formulation as presented above. We briefly sketch this approach for an incompressible neo-Hookean material. The starting point is the three-dimensional augmented free energy given in Eq. (21) together with incompressibility condition detC 3 = 1 or equivalently Then, a plane stress assumption for the total second Piola-Kirchhoff stress tensor S 3 = S 2 enables to compute the Lagrange multiplier p. Eventually, one finds the plane stress augmented free energy as in which the material electric field vector has been assumed as E 3 = E N. The plane part S 2 of the total second Piola-Kirchhoff stress tensor and the thickness component D of the material electric displacement vector follow as Finally, the plane part of the right Cauchy-Green tensor C 2 is approximated by means of the strain measures ε and κ as C 2 = 2(ε + Z κ) + I, where Z denotes the thickness coordinate; then, numerical integration of Ω 2 through the thickness finds the plate augmented free energy as Ω = Ω(ε, κ, V ). This approach has been widely used in the literature, see, e.g., [39] or [43], as well as in our own previous work [38]. We note that, according to [43], it is obvious that the models constructed in 2D manifold form, i.e., directly for the material surface, are much more efficient than models requiring numerical integration. Therefore, we derive the augmented free energy Ω = Ω(C, K, V ) directly for the material surface without the need to involve the three-dimensional formulation. This latter aspect represents the main scientific novelty of this paper. In the following, we will only concern ourselves with thin plates with thickness h made of isotropic and incompressible hyperelastic materials, for which the material is assumed homogenous through the thickness and obeys a neo-Hookean law. In order to introduce the specific form of the plate augmented free energy, we start with an additive decomposition of Ω = Ω(C, K, V ) into a mechanical part and a total electrical part, Note that the total electrical part is composed of a purely electric part of the free energy plus an augmentation term to account for the contribution from vacuum.

Strain energy
We start with a discussion of the mechanical part, i.e., the strain energy, for which the polar decomposition F = R · U enables a further additive decomposition into a membrane and a bending part, The membrane part is taken in analogy to a plane stress incompressible neo-Hookean strain energy as Here, A = Y h(1 − ν 2 ) −1 is the membrane stiffness well known from linear plate theory. With ν = 0.5 we have Young's modulus as Y = 3μ, and A = 4μh holds. η 0 is the mass per unit reference area and w m the membrane energy per unit mass. In order to introduce the bending energy, we recall the polar decomposition F = R · U, by means of which the deformation gradient tensor is multiplicatively decomposed into a plane part U and an orthogonal part R. As the first part does not contribute to the bending deformation, the bending energy must not directly depend on the stretch tensor U, and the corresponding intermediate configuration must be free of bending stresses. Therefore, the bending energy should be formulated in terms of the curvature tensor only, and be referred to the mass η per unit area in the intermediate configuration with a plate stiffnessD accounting for the thickness of the plate in the intermediate configuration. Hence, we write the bending energy per unit area in the intermediate configuration in analogy to the one of an isotropic incompressible Kirchhoff plate as Here, we have taken the standpoint that the curvature tensor is of an order of smallness that justifies the use of an incompressible Saint Venant-Kirchhoff strain energy. With the area change from the reference configuration to the intermediate configuration J = detU = detF, η = J −1 η 0 holds; obviously, η also represents the mass per unit actual area. Noting in addition the three-dimensional incompressibility condition for the bending energy.
Total electrical energy In a homogenous plate with a homogenous electric field through the thickness, the contribution from the total electrical free energy is written in analogy to a capacitor as 2ηΨ el,tot = −cV 2 , with the voltage V and the capacityc and mass η per unit actual area; the latter is identical to the mass per unit area in the intermediate configuration emerging by means of U.c is related to the capacity per unit reference area c byc = J c. Therefore, we have for the electrically homogenous case. In such a case, the total electrical free energy only contributes to the membrane part of the augmented free energy, which only depends on C. In the more general scenario of a non-homogenous electric field through the thickness of the plate, the total electrical free energy must also contribute to the bending part of the augmented free energy; therefore, it must depend also on K. In order to take this contribution into account, we propose to extend the energy w b defined in Eq. (42) with an electrically motivated source term. Then, the total bending part of the augmented free energy is which formally constitutes the classical bending energy of an isotropic Kirchhoff plate with a so-called Eigenspannungsquellem * ; the second term represents the bending part of the total electrical free energy.
Here, we have naturally extended our previous standpoint that the curvature tensor is of an order of smallness that justifies the use of an incompressible Saint Venant-Kirchhoff strain energy to the electrical energy. Keeping in mind thatm * is an Eigenspannungsquelle per unit area in the intermediate configuration, it relates to the Eigenspannungsquelle per unit area in the reference configuration m * by virtue ofm * = J m * . Therefore, holds. The first part is the bending energy from above, whereas the second part can be combined with electrical energy introduced for a homogenous electric field to As m * must be proportional to V 2 and the capacity per unit reference area c is a constant, we write m * = cαV 2 and obtain the final form of the total electrical contribution to the augmented free energy as in which α is a parameter characterizing the non-homogeneity of the electric field through the thickness. Summarizing our result, the augmented free energy of an isotropic, elastically homogeneous and incompressible hyperelastic neo-Hookean electro-elastic plate is proposed as Finally, we note that the relations trK = tr(U −1 ·κ ·U −1 ) = tr(C −1 ·κ), detK = det(U −1 ·κ ·U −1 ) = det(C −1 ·κ) and C = I + 2ε enable us to write the augmented free energy as Ω = Ω(ε, κ, V ), which closes the theory of electro-elastic material surfaces for the specific case under consideration in this paper.

Constitutive relations revisited
Although the theory of electro-elastic material surfaces for the specific case under consideration is complete, we take the time to discuss the constitutive relations in some more detail. Recalling the augmented free energy as Ω = Ω(C, K, V ), we introduce alternative stress measures such that the variation of Ω is In order to relate n with τ and m with μ, we keep the two identities in mind, and re-formulate the variation of the augmented free energy to with the Biot stress measuresn = sym (n · U) andτ = sym (τ · U). This establishes the two relations The alternative stress measures n and m are particularly useful to discuss the constitutive processes involved in the present formulation. We start this discussion by recalling the polar decomposition F = R · U, by means of which a plane intermediate configuration resulting from the stretch tensor U only is introduced, see Fig. 1.
In the reference configuration, the material surface is plane and undeformed with the strain measures C = I and K = 0, in the intermediate configuration it experiences only the right Cauchy-Green tensor C = U 2 with K = 0, and the actual configuration conserves the right Cauchy-Green tensor C = U 2 , and due to the non-trivial R e f e r e n c e c o n fi g u r a t io n

Intermediate configuration
Actual configuration Besides this purely geometrical look at the polar decomposition, more insight into the constitutive processes is gained, if the augmented free energy is additively decomposed into a membrane part and a bending part, Ω(C, , see again Fig. 1. We now proceed with discussing a proper decomposition of the stress measures n and m as well as of the internal charge q. For that sake we take the standpoint that the voltage is present in all three configurations, the reference one, the intermediate one, and the actual configuration. Therefore, it is near at hand to introduce the charge q 0 = cV , which would be the resulting charge in case the electro-elastic material surface were assumed as rigid.
Internal charge The total internal charge is with Ω el,tot = Ω el,m + Ω el,b as defined in Eq. (55). We immediately conclude that hence, q 0 can be assigned to the reference configuration with C = I and K = 0. In the intermediate configuration with C = U 2 and K = 0, the internal charge changes to and it can be computed from the membrane part Ω el,m of the total electrical contribution to the augmented free energy. This change of the charge is related to the fact that due to U the effective capacity is changed, and hence, the charge per unit area in the reference configuration must change as well. Therefore, the internal charge q in the actual configuration is in which the additional contribution q b results from the bending part Ω el,b of the total electrical contribution to the augmented free energy. Moreover, q b = 0, if the geometry parameter α characterizing the inhomogeneity of the electric field through the thickness is zero.
Electrostatic stress measures We proceed in our discussion with the electrostatic stress measures n es and m es , which result from the total electrical contribution to the augmented free energy as We start with n es and find which can be either computed from the membrane part Ω el,m of the total electrical contribution to the augmented free energy or from the internal charge q m in the intermediate configuration. Moving further to the actual configuration, an additional contribution to n es comes into the picture, such that n es = n es,m + n es,b with which can be computed from either Ω el,b or q b . Now, we turn our attention to the second electrostatic stress measure m es . We first recall the definition of the Eigenspannungsquelle m * per unit area in the reference configuration as m * = cαV 2 = αq 0 V . Then, we compute m es from Ω el,b as As m * has been introduced as an Eigenspannungsquelle, it is present in any of the three configurations as long as a voltage is applied, and as long as α is not trivial. This fact follows from the classical interpretation of an Eigenspannungsquelle as a source of self-stress acting in an elastic background structure, see the original contributions [44,45]. The electrostatic stress measure m es resulting from m * is m es = −m * I in the reference configuration and it is m es = m es,m = −m * I(detC) in both, the intermediate and the actual configuration.

Elastic stress measures
The elastic stress measures n e and m e follow from the purely mechanical part of the augmented free energy, the strain energy, which we decomposed into a membrane part and a bending part, W = W m (C)+ W b (C, K). Accordingly, the reference configuration is free from elastic stress measures n e = 0 and m e = 0. We compute the overall membrane elastic stress measure n e as and m e = 0. Hence, the intermediate configuration is free of any bending stresses; this is not surprising, as this fact was actually used to derive the bending part of the strain energy. The bending elastic stress measure m e then evolves only for a non-trivial K, such that in the actual configuration we have for the bending stress measure. In addition a further contribution to the membrane stress measure shows up, because our approach to constitutive modeling involved a plate stiffness, which depends on the strain measure C. This contribution is such that n e = n e,m + n e,b is composed of two contributions-one from W m and one from W b , whereas m e = m e,b only involves the contribution from W b . We summarize our results and represent them graphically in Fig. 2, by assigning specific parts of the stress measures n and m as well as of the internal charge q to the three configurations. Actual configuration

Small strain regime
Finally, we approximate the augmented free energy in the small strain regime, in which we set C = I + 2λε with λ as a formal small parameter; also κ is formally replaced by λκ. Concerning the voltage V , we assume its square to be of order λ. Then, we expand the augmented free energy in the vicinity of λ = 0 up to terms of order λ 2 , and setting λ → 1 we find This resembles the physically linear constitutive theory for isotropic and incompressible dielectric elastomers, insofar as the resulting constitutive relations for the total stress measures are linear in the strain measures. On the other hand, the internal charge q is a nonlinear function of the voltage and the strain measures, as a bias voltage is required to enable a sensing effect in dielectric elastomers. The mechanical part of the approximated augmented free energy corresponds to the Koiter energy [46] for such a material, and the electrical part is a straightforward extension toward electro-elasticity. Besides this interpretation, the approximated augmented free energy also enables us to identify the material parameters of the nonlinear theory. For that sake, we recall the linear theory of Kirchhoff plates with eigenstrains, see [47]; for the case of an isotropic incompressible material obeying Hooke's law the strain energy reads in which the actuation enters classically by means of Eigenspannungsquellen τ * and μ * . Comparing the approximated augmented free energy with the linear strain energy, we conclude on τ * = cV 2 and μ * = αcV 2 = ατ * . We discuss two scenarios of importance for the present paper in detail; a single-layer plate with a homogenous electric field and a bi-morph plate with a homogenous electric field only present in one of the two identical layers.
Single-layer plate For a single-layer plate of thickness h the extensional and the plate stiffness for an incompressible and isotropic material are defined as as ν = 0.5 and Y = 3μ hold. This confirms our original definition of the stiffnesses for the nonlinear theory. The Eigenspannungsquellen τ * and μ * in such a single-layer scenario are known as which confirms the meaning of c as the capacity per unit reference area. Moreover, we have for the nonlinear theory.
Bi-morph plate For a bi-morph plate of thickness h with two identical layers, the stiffnesses need no further discussion. Concerning τ * and μ * we assume one layer-specifically, the top one-to be inactive, and a homogenous electric field in the bottom layer. Then, from which the capacity c and the parameter α follow as With the material parameters identified, one can compute solutions for the nonlinear electro-elastic plate problem.

Plate finite elements
In this Section, we briefly introduce the implementation of finite elements for the numerical solution of the electro-mechanically coupled plate problem. The present version of the classical Kirchhoff-Love theory of plates requires in general C 1 continuity in the approximation of the deformed surface, which we achieve using a four-node finite element with the following approximation scheme: 16 shape functions for each spatial component of the position vector exactly represent any bi-cubic polynomial; the element thus has 48 mechanical degrees of freedom. Despite inherent restrictions concerning the topology of the mesh and connections between shell segments, the present finite element has a relatively broad spectrum of potential applications with respect to both, research and development. For more details concerning this shell element, we refer to [32].

Implementation
The starting point for the implementation is the functional Σ = Σ(ε, κ, V, p), which follows from the principle of virtual work, see Eq. (23), in which we assume no external forces and moments acting at the boundary, no external moments applied in the domain of the plate, and the voltage to be prescribed. Hence, the electrical contribution to the surface tractions at the vertical portions of the boundary of the plate is neglected; this is justified, as has been shown in [38]. Then, the functional can be written as with the total plate augmented free energy and with the potential energy Σ ext of the external force loadings p. The plate is modeled as a continuum of material normal vectors, and its configuration is defined by the field r(q 1 , q 2 ) ≡ r(q α ); −1 ≤ q α ≤ 1 are two material coordinates on a finite element. The unit normal vector results from the natural basis: while the differential operator on the surface features the co-basis: We denote the four nodes of the finite element as i, j, k, and l and write the position vector as This approximation features 12 nodal degrees of freedom: position vector r m , its derivatives with respect to the local coordinates on the element r m 1 and r m 2 , and the mixed second order derivative r m 12 . The conditions of smooth coupling with the approximation on the neighboring elements lead to a unique set of the 16 bi-cubic shape functions S m,n (q α ). The element itself is isoparametric: the reference geometry R is also approximated by means of Eq. (85) and is thus C 1 continuous. The validity of the presented approximation requires that the coordinate lines q α = const are continuous across the element boundaries. The element degrees of freedom are collected in a global vector, which contains all mechanical degrees of freedom, but no electrical degrees of freedom as we only consider problems with applied voltages.
Finally, we seek for a stationary value of the total energy functional, Σ = Σ Ω +Σ ext . Here, Σ = Σ(r(q α )) depends on the field r(q α ), which is approximated using Eq. (85). To compute the total energy functional numerical integration over the elements, with 3 × 3 integration points per element, is used. This results in nonlinear algebraic equations, from which equilibrium solutions are computed numerically by employing Newton's method for seeking the stationary points of Σ.

Boundary conditions
If an edge is free from kinematic constraints, then the external force factors acting on that edge need to be accounted for; as we are only considering problems without external loadings at such edges this is trivial. At a simply supported edge, the material points are fixed by appropriate penalty terms for the nodal positions r m and the derivatives r m α (α corresponds to the direction along the edge). If the edge is clamped, then the direction of the normal vector n needs to be additionally constrained. For a straight edge n = const, and the constraint will be fulfilled exactly, if we demand N · r m β = 0 and N · r m 12 = 0, in which β corresponds to the direction of the outer normal to the boundary of the domain.

Solid finite elements
As we will be using the three-dimensional formulation to validate the plate formulation, solid finite elements are implemented for the incompressible neo-Hookean dielectric elastomer with the augmented free energy as given in Eq. (21). To account for the incompressibility of the material response, a mixed formulation, in which the displacement field u 3 from the referential position R 3 to the actual position r 3 , i.e., u 3 = r 3 − R 3 , and the pressure p serve as variables (up form), see, e.g., Zienkiewicz et al. [48]. In the electric domain, we employ the electric potential ϕ 3 as independent variable, from which the material electric field is obtained as its material gradient, E 3 = −∇ 3 ϕ 3 .

Implementation
For the numerical analysis, we use the open-source multi-purpose finite element code Netgen/NGSolve. 1 Due to their thin-walled nature, the structures in the subsequent problems are discretized with prism elements, which are naturally aligned with the thickness direction. Both the components of the displacement field and the electric potential are approximated by means of hierarchical quadratic shape functions, whereas the pressure field is linearly interpolated. To reduce the number of unknowns, we employ symmetry conditions wherever applicable. Throughout the subsequent examples, a discretization of the thickness direction into eight equidistant layers of elements has proven to be sufficient. Regarding the (unstructured) triangular in-plane discretization, the maximum length of an element edge is restricted to be not larger than a tenth of the smaller in-plane dimension of the respective structure (with symmetry already accounted for). We use Newton's method to solve the nonlinear boundary value problem that is obtained by requiring stationarity of the augmented free energy in Eq. (21).

Boundary conditions
For the three-dimensional model, simply supported boundary conditions are realized by constraining the displacements along the respective edge on the center plane of the structure. In the proposed plate formulation, the displacement of the center plane is prohibited at clamped boundaries, whereas a deformation in thickness direction is not constrained. For this reason, the displacement in the thickness direction is only constrained at the edge on the center plane of the structure. Any displacement parallel to the center plane, however, is prohibited unless otherwise stated.

Validation
In this Section, we study example problems involving an incompressible and isotropic dielectric elastomer, for which the Lamé parameter μ was computed from the material parameters of a three parameter Ogden material used by Klinkel et al [25], see Table 1; hence, 2μ = μ 1 α 1 + μ 2 α 2 + μ 3 α 3 holds.

Stability of a single-layer dielectric elastomer plate
We consider a single-layer rectangular plate with dimension a × b × h = 4 mm × 2 mm × 0.01 mm, which is clamped at x = 0 and x = a, and free at the other two edges at y = ±b/2. At the two clamped edges the displacement in the y-direction parallel to the clamped edges is not constrained. First, we compute the critical buckling voltage at which the plane configuration loses stability and show the results in Table 2 for the present theory (ShellFE3) and the three-dimensional formulation (Netgen/NGSolve) in comparison with the result given in [25] for the Ogden-type material. A very good agreement is observed. In [25], hexahedral solid shell elements with eight nodes and four nodal degrees of freedom are used, which are the three displacements and the electric potential; the problem at hand was solved with 12 × 10 × 1 elements in [25]. Converged solutions with 16 × 8 elements in ShellFE3 and a total of 4872 elements in Netgen/NGSolve were computed. Next the voltage V is increased starting with 0V, and the non-dimensional deflection w mid / h of the center point of the plate is presented in Fig. 3. An initial imperfection by means of a pre-deformation with a corresponding non-dimensional center point deflection w mid,initial / h = 0.01 is taken into account, and w mid / h is shown relative to w mid,initial / h. In the three-dimensional model, an imperfection is introduced by means of a transverse surface load of 1 × 10 −6 Nm −2 , which is imposed on the center surface. In the left Figure of Fig. 3 we show the present results (solid line) and the ones computed with Netgen/NGSolve (solid circular markers), which are in very good agreement up to a voltage of V = 107.3 V, at which the volume elements in Netgen/NGSolve failed to converge. In contrast, in the present formulation no convergence problems have been observed up to a critical voltage of V critical = 153.5 V, which corresponds to the electro-mechanical breakdown voltage V breakdown = E breakdown h = 0.687h √ μ/ε = 153.215 V, see [38] for details. Numerical results reported by Pechstein [24] using three-dimensional Tangential Displacement Normal Normal Stress (TDNNS) mixed finite elements are shown in the right Figure of Fig. 3 up to a voltage of V = 150 V; in these results only a quarter of the problem was discretized with 458 tetrahedral elements, and a nearly incompressible neo-Hookean material with the Lamé parameters μ = 20,698 Pa and λ = 100 MPa was assumed. Again, a very good agreement with the present results is observed. Figure 4 shows the non-dimensional center point  deflection in the vicinity of the critical buckling voltage V crit,buckling , at which the supercritical pitchfork bifurcation occurs, for a significantly smaller imperfection w mid,initial / h = 0.0001. The deviation between the results computed with ShellFE3 and the ones with Netgen/NGSolve are due to the fact that we have not changed the transverse surface load of 1 × 10 −6 Nm −2 in Netgen/NGSolve, which represents a much larger imperfection in comparison with w mid,initial / h = 0.0001. Finally, the center point deflection for V = 125 V computed with the present theory is compared to the results reported by Pechstein [24] and Klinkel et al. [25] in Table 3.

Bending of a bi-layer dielectric elastomer plate
As a second example, we are studying a rectangular plate with dimension a×b×h = 100 mm×50 mm×1 mm made of two perfectly bonded dielectric elastomer layers. The electrode at the bonded interface is grounded, and a voltage can be applied at the outer electrodes. In particular, the voltage is applied at the lower electrode, and the upper one is grounded as well; this configuration will result in a bending actuator. The plate is clamped at x = 0 and free at the other edges. We use 16 × 8 elements in the plate formulation and a total of 9272 elements in the three-dimensional formulation for computing converged solutions. At increasing voltage, the actuator bends out of plane, see Fig. 5 for deformed configurations at different voltages. The voltage was first increased with a step size of ΔV = 10 V, which was refined to ΔV = 1 V for voltages above 4000 V. Results   [25] reported a value of V = 4200 V for the Ogden-type material using 32 × 18 × 2 elements; the deviation of −3.71% is comparable to the deviation we observed for the center point deflection in the first example, see Table 3.

A non-symmetric stability problem
In the last example, we study the non-symmetric problem of a bi-morph plate, with simultaneously acting electrical and mechanical loads. The plate is composed of two identical and perfectly bonded layers, whereby only one layer is electrically actuated. The plate has the dimension a × b × h = 4 mm × 2 mm × 0.02 mm with simply supported boundary conditions at all four edges. Each individual layer has the thickness 0.01 mm. A voltage V is applied to the outer electrode at the bottom layer while the interface is grounded and the top layer is short circuited. In addition, a spatially constant mechanical pressure p = (ρgh)λ with the mass density ρ = 1100 kg m −3 and the load factor λ is applied in the negative thickness direction pointing upwards. Converged solutions are computed with 16 × 8 elements in ShellFE3 and a total of 3432 elements in Netgen/NGSolve. Two load scenarios are investigated. Either the load factor λ = (0, 1, 2)λ crit is held constant and the voltage V is varied, or the load factor λ is varied for constant voltages V = (0, 1, 2)V crit . The values λ crit and V crit have the following specific meaning. For values (λ, V ) < (λ crit , V crit ) only one stable equilibrium path exists, see the top row in Fig. 7 and for values (λ, V ) > (λ crit , V crit ) we obtained two stable solution branches, see the bottom row of Fig. 7. For the critical values (λ crit , V crit ) = (0.0122, 6.41 V) a single equilibrium path with a saddle point is found, see the center row of Fig. 7.
First, we discuss the stability behavior for a constant voltage V > V crit and for a varying mechanical pressure in more detail; see also the right graph in the bottom row of Fig. 7. Due to the applied voltage, the plate deforms into a configuration with center point deflections w mid > 0. Applying a mechanical pressure in negative thickness direction w mid decreases continuously until snapping to a configuration with negative values of w mid occurs. Upon decreasing the pressure back-snapping occurs at a pressure value smaller than the value for snapping. This behavior occurs for any values V > V crit with the critical pressure values depending intrinsically on the magnitude of the applied voltage. Next, the pressure load factor λ is held constant while the voltage V is varied; in this case bi-stable solutions exist for λ > λ crit , see the left graph in the bottom row of Fig. 7. Again snapping upon loading and back-snapping upon unloading is observed. Figure 7 shows results computed with ShellFE3 (solid lines) and Netgen/NGSolve (solid circular markers) with a very good agreement between them. A practical application of the bi-stable behavior might be, for example, a voltage sensitive switching device, where the necessary pressure level for snapping can be adjusted by choosing the thickness of the dielectric elastomer switch properly.
Finally, we note that a structural mechanics framework-as developed in this paper-is especially valuable when it comes to the analysis of stability problems using a semi-analytical approach to account for the dominating nonlinear terms, for example, studying the effect of geometric nonlinearities by simplified kinematics, e.g., using a von Kármán-type theory, combined with a low-order Ritz approximation. Examples of this approach for the case of piezoelectric or thermoelastic plates can be found in [49,50] or [51] identifying that in similar examples not only buckling and snap-through / snap-back, but also snap buckling behavior can be observed.

Conclusions
In this paper we presented a complete direct approach for modeling geometrically and physically nonlinear dielectric elastomer plates as two-dimensional electro-elastic material surfaces. In contrast to common approaches, our formulation does not rely on the three-dimensional theory, but it directly develops the twodimensional one. In particular, the polar decomposition of the surface deformation gradient tensor was used to postulate the augmented free energy of the material surface.
Numerical results computed with plate finite elements for the present novel theory were compared to 3D finite elements implemented in Netgen/NGSolve, which are based on the three-dimensional counterpart to our theory. A very good agreement between the results was found in all example problems.