Nonlinear Finite Element Simulation of Thin Dielectric Elastomer Structures Sandro Zwecker #1, Sven Klinkel #2, R. Müller*3 # Statik und Dynamik der Tragwerke, TU Kaiserslautern Paul-Ehrlich-Str. 14, 67663 Kaiserslautern, Germany 1 zwecker@rhrk.uni-kl.de 2 klinkel@rhrk.uni-kl.de * Lehrstuhl für Technische Mechanik, TU Kaiserslautern Postfach 30 49, 67653 Kaiserslautern, Germany 3 ram@rhrk.uni-kl.de Abstract— To simulate the behavior of thin dielectric elastomer for the electric potential, electric field and dielectric structures a finite solid shell element formulation is presented. displacements. As nodal degrees of freedom three Dielectric elastomers belong to the group of electroactive displacements and the electric potential are assumed. The polymers and their use as actuators is caused by the efficient mixed formulation allows for a consistent finite element coupling between electrical energy input and mechanical energy approximation to avoid electromechanical locking effects. The output. Also the large elongation strain of 120-380% of the dielectric elastomer actuators and their light weight are element formulation is able to simulate large deformations. advantages that make the material very attractive. Regarding the Some numerical examples show the applicability of the electro-mechanical coupling a constitutive model is expounded. proposed solid shell element. For the definition of an electric stress tensor and a total stress tensor the electrical body force and couple are considered in the II. KINEMATICS balance of linear momentum and angular momentum, Let Φ be a deformation that point maps ⃗⃗ of the reference respectively. The governing constitutive equations are derived and incorporated in a solid shell element formulation based on a configuration to ⃗ of the current configuration at time , Hu-Washizu mixed variational principle considering six fields: see Figure 1. The deformation gradient is declared as the displacements, electric potential, strains, electric field, tangent to Φ and given by mechanical stresses, and dielectric displacements. This ⃗ formulation allows large deformations and accounts for physical ⃗⃗ . (1) nonlinearities to capture the main characteristics of dielectric elastomers. With this definition for the deformation gradient the right Cauchy-Green deformation tensor reads, I. INTRODUCTION . (2) In recent years dielectric elastomers (DE) have become popular for the usage in actuators. The efficient Followed by the Green-Lagrange strain tensor defined as transformation of electrical energy in mechanical energy and the ability to maintain large strains makes them very attractive. ( ). (3) 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 Fig. 1 Reference and current configuration with position vectors ⃗⃗ and ⃗ 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 To describe the shell formulation convective coordinates element formulation is introduced. It incorporates specific are introduced, where and are the in-plane coordinates interpolations for the displacements, strains, stresses as well as and is the coordinate in thickness direction. The covariant tangent vectors are obtained in the reference and current ⃗ (⃗ ) ⃗ (⃗ ) , ⃗- ⃗ by neglecting higher ⃗⃗ ⃗ configuration as ⃗⃗ , ⃗⃗ respectively. The order terms results in contravariant basis vectors are defined in a standard manner ⃗ (∑ )⃗ (∑ )⃗ , ⃗ - (∑ ⃗ ). by ⃗⃗ ⃗⃗ and ⃗⃗ ⃗⃗ , where is the Kronecker- (7) delta. With this convective description the deformation gradient also reads ⃗⃗ ⃗⃗ . After defining the resultant mass ∑ , charge ∑ , and polarization ⃗ ∑ ⃗ A potential is used to describe the electric field ⃗ . With of the particle and performing a simple space averaging ∑ , satisfying the Laplace’s equation ⃗ reads ∑ ,⃗ ∑ ⃗ leads to the macroscopic body ⃗ ⃗⃗ . (4) force density ⃗ ⃗ ⃗ ⃗ , ⃗- ⃗ ⃗ , (8) Applying the pull-back operation the physical electric field ⃗ ⃗ is observed. The displacement vector ⃗ is defined where ⃗ combines the electric contribution to the body by force density. Deriving the corresponding couple density with ⃗⃗ . the same arguments the resultant couple of the particle ⃗ ⃗ (5) with respect to the mass center reads Boundary conditions for ⃗ and are given on and , respectively. ⃗ ∑ ⃗ ⃗ ∑ ⃗ ⃗ . (9) Here it is assumed that the mass dipole moment III. FORCE AND COUPLE DENSITY ∑ ⃗ ⃗ . With the resultant polarization ⃗ To get a macroscopic representation of force and couple ⃗ and a space averaging it follows the macroscopic ∑ density a close look on the microscopic level is presented. couple density Therefore the physical model will start with a particle of the current configuration to deduce electric force in the ⃗ ⃗ ⃗ . (10) presence of electromagnetic matter and in absence of a magnetic field. IV. BALANCE LAWS AND STRESS TENSOR The balance laws incorporating the electric force, couple, and power densities are summarized. After that the global integral forms and the local field equations are presented. The Cauchy stress tensor and an electric stress tensor are introduced. Conservation of mass ∫ , with is assumed. Localization of the material description results in . In quasi static processes the integral form of the balance of linear momentum is given as ∫ ⃗ ∫ , (11) Fig. 2 Physical model for the microscopic description; left: the current configuration with a particle ,right: the zoomed inner structure of the where ⃗ is the body force density (8) and the traction particle. 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 Relative to the mass center eccentric by ⃗ not to be conservation of mass, applying the divergence theorem and confused with the convective coordinates are the point the localization theorem, the field equation along with the charges and the point masses , as shown in Figure 2. jump condition are observed as In an electric field a force is acting on each point charge called Lorentz force and it is defined by ⃗ . The ⃗ ⃗ in , (12) gravitation field ⃗ is acting on each point mass producing the Newton force given by ⃗ . Summing over leads to the ⃗ on . (13) resultant force on the particle: With as a prescribed traction on the boundary . The ⃗ ∑ ⃗ ⃗ . (6) boundaries with a given traction and a given displacement satisfy ⋃ and ⋂ . The global The position of a point within the particle can be form of the balance of angular momentum reads described by ⃗ ⃗ ⃗ . Assuming the gravitation to be constant it follows: ⃗ (⃗ ) ⃗ (⃗ ). This assumption and ∫ ⃗ ⃗ ⃗ ∫ ⃗ , (14) expanding the external field in a Taylor series where ⃗ is the couple density Using the integral V. CONSTITUTIVE EQUATIONS theorem, considering conservation of mass and the linear momentum balance along with the localization theorem Introducing the energy function results in the field equation (⃗ ⃗) [(⃗ ⃗) ] , ⃗ ⃗ in . (15) (28) For dielectric materials the conservation of charge in where is a function of to fulfill material objectivity, integral form with the surface charge density on reads here an Ogden-type material is chosen, and the susceptibility of the material is denoted by , the total stress and the ∫ ∫ . (16) dielectric displacements are derived as The dielectric displacement vector is denoted by ⃗ and , (29) determined by Gauss’ law ⃗ ⃗ . Applying the divergence theorem results in the field equation along with the ⃗⃗ . (30) ⃗⃗ jump condition ⃗ VI. WEAK FORMULATION in , (17) In this section a mixed variational formulation is introduced. ⃗ ⃗ on . (18) Let the set { ⃗ , ( )- ⃗ ⃗ } be the space of Here, is a prescribed surface charge. For the total admissible displacement variations and { , ( )- boundaries with a given electric potential and a given surface } be the space of admissible electric potential charge it holds ⋃ and ⋂ . variations. Further let * , ( )- + , * With the constitutive equations in matter the dielectric , ( )- + the spaces of admissible variations of the total displacements are given as ⃗ ⃗ ⃗ , where is the stresses and strains and { ⃗⃗ , ( )- } , { ⃗ , ( )- } permittivity in vacuum and the polarization depends on the the spaces of admissible variations of the dielectric considered material. An electric stress tensor is introduced as displacement and the electric field. Since the variations are arbitrary the field equations (24)-(27), the constitutive ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ , (19) equations (29), (20), and the kinetic field equations (3), (4) are rewritten 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 ∫ ( ( )) ,  ( ) ⃗ ⃗ in ,  ∫ . ⃗⃗ ⃗⃗ / ⃗ ,  ⃗ on ,  ⃗ in ,  ∫ ( ( )) ( ) ,  ⃗ ⃗  on . ∫ (⃗ ) ⃗⃗ .  nd Considering the 2 Piola-Kirchhoff stress tensor Applying integration by part, using the divergence theorem and and , the pull-back of the dielectric displacement considering the boundary conditions the weak formulation ⃗⃗ ⃗ and the transformation of the densities by leads reads to the material description , ( )- ⃗ ⃗ in ,  ∫ ( ) ( ) ⃗⃗ ⃗⃗ ⃗ ⃗ ⃗ ∫ ⃗ ⃗ on ,  ∫ ( ) ( ) ⃗⃗ ⃗⃗ in ,  ∫ ⃗ ( ) ⃗ ⃗⃗ ⃗⃗ ⃗⃗ ⃗⃗ on .  ( ) ⃗⃗ ⃗ ,  Where ⃗ is the traction with respect to the reference configuration. with ⃗ , ( ) , and ⃗ ⃗ ⃗ ,  . with , - and [ ]. VII. FINITE ELEMENT APPROXIMATION The matrix is defined with some ANS interpolations in [6] In this section a solid shell element is introduced. The finite and reads element approximation is constructed in the sense that the whole domain is divided in element domains with ⋃ , where is the total number of elements. The ⃗⃗ geometry, displacements, and electric potential are ⃗⃗ approximated as ⃗⃗ ∑ ⃗⃗ , ⃗ ∑ ⃗ , and ∑ ( )( ) ⃗⃗ , ∑ with the same interpolation functions ⃗⃗ ⃗⃗ ( ). ⃗⃗ ⃗⃗ / ( ). ⃗⃗ ⃗⃗ / ( )( )( ), at [ ( ). ⃗⃗ ⃗⃗ / ( ). ⃗⃗ ⃗⃗ /] nodes . The vectors ⃗⃗ and ⃗ contain the nodal coordinates and displacements, respectively. Arranging in  the matrix , - with , - , the virtual quantities are where denotes the partial derivative of the shape function interpolated as with respect to the curvilinear coordinates. The matrix at the node is determined as ⃗ ⃗ [ ] ⃗ ,  , -.  where ⃗ [⃗ ⃗ ⃗ ⃗ ] is the vector of nodal The physical stresses and dielectric displacements are derived ⃗ ,⃗ - from the potential function and are arranged in the vector degrees of freedom ⃗ ⃗ . Accordingly, ⃗ is the vector of the virtual values. 0 ⃗ 1 . Here, ⃗ are independently ⃗ ⃗ ⃗ The gradient fields are defined with respect to the curvilinear assumed quantities for strain and electric field components and coordinates . The constitutive equations  and  will are approximated with the following interpolations, see Klinkel be given with respect to a local orthonormal coordinate system and Wagner [7]: ⃗ . This necessitates a transformation of the strains and the ⃗ ⃗⃗ ⃗⃗ ,  electric field. Introducing the transformation matrix ( ) ( ) ( ) ( ) ( ) ( ) with [ ], [ ] , ⃗⃗ , and ⃗ ⃗ ( ) ( ) ( ) , ⃗⃗ . The matrices and ⃗ are given as [ ]  ( ) , where defined with , and ⃗⃗ ⃗ [ ] ⃗⃗ ⃗⃗ ⃗⃗ and ⃗ , ⃗ , ⃗ ⃗ ⃗ . The Jacobian  ‖⃗⃗ ‖ ‖⃗⃗ ⃗⃗ ‖ matrix is denoted by ⃗⃗ ⃗ . For the sake of a compact notation the contravariant components of the virtual strain ⃗ [ ]. tensor and the virtual electric field vector are arranged in a generalized vector  ⃗⃗ ⃗⃗ ⃗⃗ ⃗⃗ 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 on element level of the virtual gradient  field ⃗ reads ⃗ ∫ ⃗ ,  ⃗ [ ].  ⃗ ∫ ⃗̃ ∫ ̃ ,  The approximation of the independent field is defined as 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 with [ ] 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:  ⃗ 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 ∫ ⃗ [ ] ∫ ⃗ ⃗ ,  ∫ ⃗ ⃗ ⃗ ⃗ ⃗ ,  ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ , with ⃗ , ⃗ - and ⃗ contains the components of , with ( ) and ⃗ ⃗ , and ⃗⃗ accordingly to the vector notation. The weak form is ( ) ⃗ . After assembly over all elements solved iteratively by employing Newton-Raphson’s method. ⋃ , ⃗ ⋃ ⃗ and ⃗ ⋃ ⃗ one This requires the linearization obtains , - ∫ ⃗ ⃗ ∫ ⃗ ⃗ ⃗ ⃗ ,  ∫ ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ ⃗ with the unknown incremental nodal displacements and the electric potential. The update of the internal degrees of +∫ ⃗ ⃗   freedom reads Considering the above interpolations  and  one obtains ⃗⃗ ⃗ ⃗ ,  the following matrices: ⃗⃗ ( ) ⃗ ( ) ⃗⃗ ,  ∫ ( ) ⃗ ⃗ ,  ⃗⃗ ⃗⃗ ⃗ .  ∫ ( ) ,  VIII. NUMERICAL EXAMPLES Embedding the solid shell element formulation in a ∫ ,  modified version of the program FEAP demonstrates the capturing of the characteristics of the DE material behavior. ∫ ,  The ability of the shell element is shown in the following examples. where is the matrix of [7] and vectors A. Eigenvalue Problem 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] , , In the simulation the surface charge is increased from to , and . and decreases back to again. Then the other The relative permittivity is given as ( ) . layer is loaded in the same way. The tip deflection versus the In Table 1 it is shown that there is only eigenvalue number applied voltage is shown in Figure 5. 18 much greater than zero. This eigenvalue accounts for the 80 volumetric locking mode which arises by incompressible 60 materials like DE. At this point no other locking modes are 40 observed. The result is also the same when taking an irregular Tip deflection [mm] non-cube element. 20 TABLE I 0 0 1000 2000 3000 4000 5000 6000 7000 EIGENVALUES OF A DIELECTRIC ELASTOMER CUBE -20 No. Eigenvalue No. Eigenvalue -40 1 0.11E-02 10 0.50E-01 2 0.31E-02 11 0.62E-01 -60 3 0.72E-02 12 0.75E-01 -80 4 0.11E-01 13 0.85E-01 Load j [V] 5 0.18E-01 14 0.87E-01 6 0.20E-01 15 0.97E-01 Fig. 5 Tip deflection versus voltage relations 7 0.21E-01 16 0.10E+00 8 0.27E-01 17 0.10E+00 IX. SUMMARY AND OUTLOOK 9 0.34E-01 18 0.57E+01 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, B. Bending Actuator stresses and dielectric displacements. The element possesses The second example is presented to demonstrate the valid only four nodal degrees of freedom, displacements and the reflection of the electro-mechanical coupling phenomenon of electric potential. It allows a consistent approximation of the DEs. Therefore a square plate with the dimension of electromechanical coupled problem. The governing field is investigated. It is clamped on one side equations and boundary conditions are presented. The and consists of two layers with the thickness of each, formulation accounts for geometric and material nonlinear see Figure 3. The dataset of the material is the same as given behavior. Further tasks will be to embed the material in the first example. formulation into an improved stabilized element formulation. For more accurate simulations other material models for the mechanical part should be discussed. ACKNOWLEDGMENT We want to thank the Landesforschungszentrum CM² for financially supporting our project Mikro- & Strukturmechanik zur Analyse des nichtlinearen Deformationsverhaltens von dielektrischen und porösen Elastomeren (MSADEL). Fig. 3 Dimensions of the bending actuator REFERENCES To get a bending answer of the thin structure either the [1] A. Dorfmann and R.W.Ogden, “Nonlinear electroelasticity”, Acta Mechanica 174, pp. 167-183, 2005. upper or lower layer is loaded by an electric field applied in [2] D.J. Steigmann, “Equilibrium theory for magnetic elastomers and thickness direction of the structure. Due to the electro- magnetoelastic membranes”, Int. J. Nonlinear Mech. 39(7), pp. 1193- mechanical coupling the loaded layer responds with an in- 1216, 2004. plane extension. Because of the eccentricity the bending effect [3] D.K. Vu, P. Steinmann and G. Possart, ”Numerical modelling of non- linear electroelasticity”, Int. J. Numer. Meth. Engng 70, pp.685-704, of the whole structure occurs. Figure 4 shows several 2007. deformed configurations. [4] A.C. Eringen and G.A. Maugin, Electrodynamics of continua I, Foundations and Solid Media, New York : Springer-Verlag, 1990. [5] R. Müller, B.X. Xu, D. Gross, M. Lyschik, D. Schrade, S. Klinkel, „Deformable dielectrics – optimization of heterogeneities“, International Journal of Engineering Science 48, pp. 647–657, 2010. [6] S. Klinkel, “A thermodynamic consistent 1d model for ferroelastic and ferroelectric hysteresis effects in piezoceramics”, Communications in Numerical Methods in Engineering, 22(7), pp.727-739, 2006. [7] S. Klinkel and W. Wagner, “A piezoelectric solid shell element based on a mixed variational formulation for geometrically linear and nonlinear applications”, Computers and Structures 86, pp. 38-46, 2008. [8] M.T. Wissler, Modeling dielectric elastomer actuators, PhD. Thesis, Fig. 4 Bending actuator - several deformed configurations Swiss Federal Institute of Technology Zürich, 2007.