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.