In recent years dielectric elastomers (DE) have become popular for the usage in actuators. The efficient transformation of electrical energy in mechanical energy and the ability to maintain large strains makes them very attractive. A Constitutive model, which describes the specific material behavior, is introduced by Dorfmann and Ogden [1], Steigmann [2] and the references therein. The numerical treatment in the context of the finite element method is presented by Vu et al. [3], who provide a brick element formulation, which incorporates the nonlinear constitutive model. In the present work a shell finite element formulation for DE is presented. It is motivated by the fact that the most DE devices are thin structures, which have a very high length to thickness ratio. The electromechanical coupling is considered in the body force and the couple density, see e.g. Eringen and Maugin [4] or Müller et al. [5]. The angular momentum equation is fulfilled by assuming a Maxwell stress tensor. The linear momentum equation is approximately fulfilled by the finite element method. A mixed solid shell element formulation is introduced. It incorporates specific interpolations for the displacements, strains, stresses as well as for the electric potential, electric field and dielectric displacements. As nodal degrees of freedom three displacements and the electric potential are assumed. The mixed formulation allows for a consistent finite element approximation to avoid electromechanical locking effects. The element formulation is able to simulate large deformations. Some numerical examples show the applicability of the proposed solid shell element.

Let Φ be a deformation that point maps ⃗⃗ of the reference configuration to ⃗ ⃗ of the current configuration at time , see Figure 1. The deformation gradient is declared as the tangent to Φ and given by ⃗ ⃗ ⃗⃗ .

(

With this definition for the deformation gradient the right Cauchy-Green deformation tensor reads, .

(

Followed by the Green-Lagrange strain tensor defined as A potential is used to describe the electric field ⃗ ⃗ . With satisfying the Laplace's equation ⃗ ⃗ reads

Applying the pull-back operation the physical electric field ⃗ ⃗ ⃗ is observed. The displacement vector ⃗ ⃗ is defined by

Boundary conditions for ⃗ ⃗ and are given on and , respectively.

To get a macroscopic representation of force and couple density a close look on the microscopic level is presented. Therefore the physical model will start with a particle of the current configuration to deduce electric force in the presence of electromagnetic matter and in absence of a magnetic field. Relative to the mass center eccentric by ⃗ not to be confused with the convective coordinates are the point charges and the point masses , as shown in Figure 2. In an electric field a force is acting on each point charge called Lorentz force and it is defined by ⃗ . The gravitation field ⃗ ⃗ is acting on each point mass producing the Newton force given by ⃗ ⃗ . Summing over leads to the resultant force on the particle:

The position of a point within the particle can be described by ⃗ ⃗ ⃗ ⃗ ⃗ . Assuming the gravitation to be constant it follows: ⃗ ⃗ ( ⃗ ⃗ ) ⃗ ⃗ ( ⃗ ⃗ ). This assumption and expanding the external field in a Taylor series ⃗ ( ⃗ ⃗ ) ⃗ ( ⃗ ⃗ ) , ⃗ -⃗ by neglecting higher order terms results in

After defining the resultant mass ∑ , charge ∑ , and polarization ⃗ ⃗ ∑ ⃗ of the particle and performing a simple space averaging ∑ , ∑ , ⃗ ⃗ ∑ ⃗ ⃗ leads to the macroscopic body force density

where ⃗ combines the electric contribution to the body force density. Deriving the corresponding couple density with the same arguments the resultant couple of the particle with respect to the mass center reads

Here it is assumed that the mass dipole moment . In quasi static processes the integral form of the balance of linear momentum is given as

where ⃗ ⃗ is the body force density (8) and the traction vector on . With the Cauchy stress tensor the traction is determined by a linear map of the normal vector ⃗ ⃗ , according to Cauchy's stress theorem. Considering conservation of mass, applying the divergence theorem and the localization theorem, the field equation along with the jump condition are observed as

With as a prescribed traction on the boundary . The boundaries with a given traction and a given displacement satisfy ⋃ and ⋂ . The global form of the balance of angular momentum reads

where ⃗ is the couple density Using the integral theorem, considering conservation of mass and the linear momentum balance along with the localization theorem results in the field equation

For dielectric materials the conservation of charge in integral form with the surface charge density on reads

The dielectric displacement vector is denoted by ⃗ ⃗ and determined by Gauss' law ⃗ ⃗ ⃗ ⃗

. Applying the divergence theorem results in the field equation along with the jump condition

Here, is a prescribed surface charge. For the total boundaries with a given electric potential and a given surface charge it holds ⋃ and ⋂ . With the constitutive equations in matter the dielectric displacements are given as ⃗ ⃗ ⃗ ⃗ ⃗ , where is the permittivity in vacuum and the polarization depends on the considered material. An electric stress tensor is introduced as

such that ⃗ and ⃗ . It is remarked that ⃗ ⃗ ⃗ ⃗ is also known as Maxwell stress tensor.

With  it follows that the total stress tensor has to be symmetric. The remaining field equations and boundary conditions are

Considering the 2 nd Piola-Kirchhoff stress tensor and , the pull-back of the dielectric displacement ⃗⃗ ⃗ ⃗ and the transformation of the densities by leads to the material description

Where ⃗ ⃗ is the traction with respect to the reference configuration.

Introducing the energy function

where is a function of to fulfill material objectivity, here an Ogden-type material is chosen, and the susceptibility of the material is denoted by , the total stress and the dielectric displacements are derived as

VI. WEAK FORMULATION In this section a mixed variational formulation is introduced. Let the set { ⃗ ⃗ , ( )-⃗ ⃗ ⃗ ⃗ } be the space of admissible displacement variations and { , ( )-} be the space of admissible electric potential variations. Further let * , ( )-+ , * , ( )-+ the spaces of admissible variations of the total stresses and strains and { ⃗⃗ , ( )-} , { ⃗ ⃗ , ( )-} the spaces of admissible variations of the dielectric displacement and the electric field. Since the variations are arbitrary the field equations ( 24)-( 27), the constitutive equations ( 29), (20), and the kinetic field equations ( 3), ( 4) are rewritten as

Applying integration by part, using the divergence theorem and considering the boundary conditions the weak formulation reads The gradient fields are defined with respect to the curvilinear coordinates . The constitutive equations  and  will be given with respect to a local orthonormal coordinate system ⃗ . This necessitates a transformation of the strains and the electric field. Introducing the transformation matrix

,

where denotes the partial derivative of the shape function with respect to the curvilinear coordinates. The matrix at the node is determined as , -.

The physical stresses and dielectric displacements are derived from the potential function and are arranged in the vector

Here, ⃗ are independently assumed quantities for strain and electric field components and are approximated with the following interpolations, see Klinkel and Wagner [7]:

, ⃗⃗ , and

. The matrices and ⃗ are given as

 Quantities, which are evaluated at the element center are denoted with the index and is a identity matrix. The transformation matrix is obtained by  considering and . The matrices and ⃗ are defined as

The approximation of the independent field is defined as

] and ⃗⃗ .

Here the matrix is equivalent to of , where instead of ( ) the transformation matrix ( ) is used. The interpolation ⃗ ⃗ is identical to ⃗ .

The weak formulation of  will be approximated on element level as following:

with ⃗ ⃗ , ⃗ ⃗ and ⃗ ⃗ contains the components of , , and ⃗⃗ accordingly to the vector notation. The weak form is solved iteratively by employing Newton-Raphson's method. This requires the linearization

Considering the above interpolations  and  one obtains the following matrices:

where is the matrix of [7] and vectors

with and . In Eq. ( 59) the body and surface loads are determined by ⃗ ⃗ , ⃗ ⃗ and ̃ 0 ⃗ ⃗ 1.

Having in mind that (50) is solved iteratively the following approximation on element level is obtained

.

Taking into account that the finite element interpolations for the fields ⃗ ⃗ and ⃗ are discontinuous across the element boundaries a condensation on element level yields the element stiffness matrix and the right hand side vector

To analyze the locking modes of the solid shell element an eigenvalue problem is solved. In this first example a cube form DE material with an edge length of is examined. The cube is zero-stress supported. The data for the Ogden-type material are according to [8] , , , and . The relative permittivity is given as ( ) . In Table 1 it is shown that there is only eigenvalue number 18 much greater than zero. This eigenvalue accounts for the volumetric locking mode which arises by incompressible materials like DE. At this point no other locking modes are observed. The result is also the same when taking an irregular non-cube element.

The second example is presented to demonstrate the valid reflection of the electro-mechanical coupling phenomenon of DEs. Therefore a square plate with the dimension of is investigated. It is clamped on one side and consists of two layers with the thickness of each, see Figure 3. The dataset of the material is the same as given in the first example. To get a bending answer of the thin structure either the upper or lower layer is loaded by an electric field applied in thickness direction of the structure. Due to the electromechanical coupling the loaded layer responds with an inplane extension. Because of the eccentricity the bending effect of the whole structure occurs. Figure 4 shows several deformed configurations. In the simulation the surface charge is increased from to and decreases back to again. Then the other layer is loaded in the same way. The tip deflection versus the applied voltage is shown in Figure 5.

The presented element formulation is based on a mixed variational approach. It results in an independent interpolation of the displacement, electric potential, strains, electric field, stresses and dielectric displacements. The element possesses only four nodal degrees of freedom, displacements and the electric potential. It allows a consistent approximation of the electromechanical coupled problem. The governing field equations and boundary conditions are presented. The formulation accounts for geometric and material nonlinear behavior. Further tasks will be to embed the material formulation into an improved stabilized element formulation. For more accurate simulations other material models for the mechanical part should be discussed.