An electromechanically coupled beam model for dielectric elastomer actuators

In this work, the Cosserat formulation of geometrically exact beam dynamics is extended by adding the electric potential as an additional degree of freedom to account for the electromechanical coupling in the Dielectric Elastomer Actuators (DEAs). To be able to generate complex beam deformations via dielectric actuator, a linear distribution of electric potential on the beam cross section is proposed. Based on this electric potential, the electric field and the strain-like electrical variable are defined for the beam, where the strain-like electrical variable is work-conjugated to the electric displacement. The electromechanically coupled strain energy for the beam is derived consistently from continuum electromechanics, which leads to the direct application of the material models in the continuum to the beam model. The electromechanically coupled problem in beam dynamics is first spatially semidiscretized by 1D finite elements and then solved via variational time integration. By applying different electrical boundary conditions, different deformations of the beam are obtained in the numerical examples, including contraction, shear, bending and torsion. The damping effect induced by the viscosity as well as the total energy of the beam are evaluated. The deformations of the electromechanically coupled beam model are compared with the results of the 3D finite element model, where a good agreement of the deformations in the beam model and that in the 3D finite element model is observed. However, less degrees of freedom are required to resolve the complex deformations in the beam model.

compliant electrodes. When an external electric field is applied to the DEA, the dielectric material will be polarized, resulting in electrostatic pressure, see the models in [Pelrine et al., 1998] , [Wissler and Mazza, 2007], [Suo et al., 2008] and [Schlögl and Leyendecker, 2017] for instance. Due to the contractive pressure, the contraction of the DEA will be induced such that it can be applied as an actuator. The deformation behavior of the DEA is governed by the electromechanical coupling in the dielectric material.
For the general investigation of the electromechanical coupling behavior, much effort has been made to address the nonlinear electroelasticity in the past years, see e.g. the theory of interaction of electromagnetic and elastic fields in deformable continua in [Pao, 1978], the nonlinear electroelasticity formulation for the finite deformation in [Dorfmann and Ogden, 2005] and the variational formulations of the electroand magneto-elastostatics in [Vu et al., 2007]. Additionally, material models of the dielectric elastomers have been investigated, see e.g. [Zhao et al., 2007], [Vu et al., 2007] and [Suo, 2010]. In [Khan et al., 2013], a viscoelastic effect is introduced to account for the damped dynamic behavior in silicon based dielectric elastomers. A viscoelastic 3D finite element model of the DEA is developed by [Schlögl and Leyendecker, 2016b] for the dynamic analysis using a structure preserving time integration scheme. This model is extended to flexible multibody system dynamics in [Schlögl and Leyendecker, 2016a].
The finite element models introduced above provide a powerful and accurate tool for solving the electromechanical coupling problem in DEA. However, huge amounts of degrees of freedom are required in large 3D finite element models, which leads to inefficient computation and difficulties for the optimal control of DEA. Especially for long thin artificial muscles, 3D finite element model is more expensive in computational cost than beam model. Further more, the coupling between 3D finite element models and rigid body is possible but has to be specially addressed in multibody system. The geometrically exact beam performs well concerning the tradeoff between computational cost and accuracy for the simulation of slender structures like the stacked DEA. By assigning the rotational degrees of freedom to points in continuum, the Cosserat formulation [Cosserat and Cosserat, 1909] of geometrically exact beam closes the gap between classical continuum mechanics and rigid (multi-)body dynamics, which leads to a consistent formulation of flexible multibody systems. The fundamental formulations on geometrically exact beam can be found in [Simo, 1985] and [Antman, 2005] for instance. The time integration of constrained geometrically nonlinear beam dynamics has been discussed by many authors, see e.g. the energy conserving/decaying algorithms in [Armero and Romero, 2001], the energy-momentum scheme with null space method in [Betsch and Leyendecker, 2006] and the variational integrators in [Leyendecker et al., 2008]. To account for the electric field in a beam, the focus has been put on the piezoelectric effect, see [Krommer and Irschik, 2002], [Tadmor and Kósa, 2003] and [Schoeftner and Buchberger, 2012]. An analytical model for dielectric elastomer based microbeam is discussed in [Feng et al., 2011]. However, the electromechanically coupled problem of dielectric elastomers in geometrically exact beam is still not given.
The objective of this work is to develop an electromechanically coupled beam model for the simulation of stacked dielectric elastomer actuators, where the Cosserat formulation of geometrically exact beam dynamics is extended by adding the electric potential as the additional degree of freedom. A linear distribution of electric potential on the beam cross section is proposed to generate different beam deformations including contraction, shear, bending and torsion. The electric field in the beam is computed from the gradient of the electric potential. Based on the formulation of deformation gradient and electric field in the beam, the electromechanically coupled strain energy function for the beam is derived consistently from the strain energy function in continuum electromechanics, which leads to the direct application of the material models in continuum electromechanics to the beam model. The viscoelastic effect is taken into account in the non-conservative force term. The electromechanically coupled problem in beam dynamics is first semidiscretized with 1D spatial finite elements and then solved via variational time integration. By applying different electrical boundary conditions, different deformations of the beam are obtained in the numerical examples, which are compared with the results of the 3D finite element model. This paper is structured as follows: In Section 2 and 3, the governing equations for electromechanical coupling in continuum electromechanics and geometrically exact beam are presented, respectively. Then the formulation of kinematic variables including the deformation gradient, the electric potential and the electric field are derived for the beam in Section 4. Section 5 presents the consistent derivation of a strain energy function for the beam from continuum electromechanics. In Section 6, the electromechanical coupling problem is solved within the variational time integration scheme with null space projection.
The numerical examples of the developed model are presented in Section 7, followed by the conclusions in Section 8.

Governing equations for electromechanical coupling in continuum electromechanics
The finite deformation of a dielectric elastic solid occupying the domain B ⊂ R 3 is distinguished by the initial and current configurations. The boundary of the solid ∂B is composed by the Dirichlet type sections ∂ u B and ∂ φ B, and the Neumann type sections ∂ σ B and ∂ D B. By denoting the position of a material point in the initial configuration with X ∈ B 0 , the position of the material point in the current configuration at time t is given by where u is the displacement. The deformations within the body induced by the electric field satisfy the balance law of momentum and the Maxwell equations.

Balance of linear and angular momentum
The local balance law of linear momentum in the dynamic process is given by subject to the Dirichlet and Neumann boundary conditions where P is the first Piola-Kirchhoff stress tensor, ρ 0 is the mass density in initial configuration,b is the body force vector,ü is the acceleration,ū is the prescribed displacement andT is the prescribed traction. The local balance of angular momentum reads in which F is the deformation gradient defined as F = ∂x(X, t)/∂X.

Maxwell equations
By neglecting the magnetic field, the Maxwell equations are given by subject to the Dirichlet and Neumann boundary conditions with E e the electric field, D the electric displacement in the initial configuration, φ the electric potential, φ the prescribed electric potential, N the outward unit normal vector andQ the prescribed charges per unit area on the boundary ∂ D B. The Eq. (6) 1 leads to the definition of the electric field as the gradient of a scalar electric potential

Electromechanical coupling
When the external electric field is imposed in the body of dielectric elastomer, the contractive pressure will be induced due to the polarization effects and thus the deformation of the body will be generated. The coupling effect between the electric field and the mechanical deformation is described by the strain energy function Ω(F, E e ) of the dielectric material in the constitutive equations For the dielectric materials, the electromechanical coupling can be described by the strain energy function with the additive form with Ω m (F) referring to the purely mechanical behavior, Ω em (F, E e ) referring to the electomechanical coupling and Ω e (E e ) referring to the pure electric behavior. Accordingly, the first Piola-Kirchhoff stress can be written as two parts

Governing equations for electromechanical coupling in geometrically exact beam
In this work, the formulation of the electromechanical coupling problem presented above is extended to the geometrically exact beam. The deformation state of an initially straight beam over time can be distinguished by the initial configuration and the current configuration, as shown in Fig. 1.
In the Cosserat formulation of geometrically exact beam, the placement of a material point in the current configuration of the beam is given by x(X, t) where s ∈ [0, L] ⊂ R denotes the arc-length of the line of the centroids ϕ(s, 0) ∈ R 3 in the initial configuration, d i (s, t) is an orthonormal triad at s with the directors d k (s, t), k = 1, 2 spanning a principle basis of the cross section, and X k are the curvilinear coordinates on the cross section. Based on these assumptions, the governing equations for electromechanically coupled beam model can be consistently derived from the equations of continuum electromechanics.

Electromechanically coupled beam model in variational setting
The governing equations of the beam are usually derived by integrating Eq.
(2) in continuum mechanics over the beam cross section, see e.g. [Simo, 1985]. However, the Neumann boundary condition on the beam is not clearly described. To obtain a consistent derivation of beam equations from continuum mechanics considering the Neumann boundary conditions in Eq. (4) and (8), in this work, we formulate the governing equations for the beam from the variational setting. The variational setting is formulated according to the Lagrange-d'Alembert principle with S the action and W ext the external non-conservative work contributed by the body force and surface traction. According to the constitutive law in Eq. (10), the variation of the action can be formulated for the beam as where the formula −∂ F Ω : δF = (∇ X · ∂ F Ω) · δx − ∇ X · (∂ F Ω · δx) and the divergence theorem have been applied. Due to the assumption of geometrically exact beam, the volume integral in Eq. (15) is split into the curve integral over beam center line and the area integral over beam cross section Σ. At the same time, the surface integral in Eq. (15) is split into the curve integral over beam center line and the curve integral over the lateral contour of beam cross section. By using of P = t k ⊗ d k (s, 0) + t s ⊗ d 3 (s, 0) and D = D k d k (s, 0) + D s d 3 (s, 0), the force f , the torque m and the electric displacement d e for the beam are defined as with r = x − ϕ. For the sake of later use, the electric displacement vector d e for the beam is defined as d e = d e 1 d e 2 d e s T with the components d e k = Σ D k dA. In this case, the divergence terms in Eq. (15) can be formulated as In the surface integral of Eq. (15), the unit normal vector is in the plane of cross section, i.e. N = n k d k (s, 0)(k = 1, 2). Thus, by applying the divergence theorem, we have ∂Σ P · Ndl = f n , ∂Σ r × (P · N)dl = m n and ∂Σ D · Ndl = d en .
The variation of the position field is given by with δη = 1 2 d i × δd i , see the details in [Eugster et al., 2014]. The time derivatives of the position field are given byẋ with ω the spatial angular velocity. Using Eq. (19) -(21), the varitation of the action in Eq. (15) reads in which I = Σ ρ x 2 − x ⊗ x dA is the spatial mass moment of inertia tensor. It can be observed that f n and m n in the first curve integral of Eq. (22) come from the divergence term of the action in Eq. (15), f n and m n in the second curve integral of Eq. (22) come from the integration over a boundary surface in Eq. (15). Due to their opposite signs, they will be offset respectively. By considering the Neumann boundary conditions in Eq. (4) and (8), the variation of the external work is written as where the prescribed body force, lateral traction, torque and lateral charge are defined respectively as Combing Eq. (22) and (23), the Lagrange-d'Alembert principle for the beam is given by

Balance of linear and angular momentum
The requirement of stationary in Eq. (25) in terms of the translation field ϕ leads to the balance of linear momentum for the beam The requirement of stationary in Eq. (25) in terms of the rotation field η leads to the balance of the angular momentum for the beam

Maxwell equation
The requirement of stationary in Eq. (25) in terms of the electric potential φ leads to the Maxwell equation for the beam

Electromechanical coupling
The force f , the torque m and the electric displacement vector d e in the beam governing equations are related with the kinematic variables via the constitutive equations where Ω b is the strain energy per unit arc-length in a beam, ε is the strain-like electrical variable conjugated with the electric displacement d e of the beam, γ and κ are the beam strain measures conjugated with f and m, respectively. The electromechanical coupling can be specified by the strain energy function in the additive form Due to the coupled term in strain energy, the force and torque are contributed by two parts accordingly, the electric part and the mechanical part, i.e.

Mechanical and electrical kinematics in the beam
To compute the force f , the torque m and the electric displacement d e in the governing equations of the electromechanically coupled beam model in Section 3, the conjugated kinematic variables for mechanical and electrical parts have to be formulated. The mechanical strain measures γ and κ for the geometrically exact beam have been defined and related with the deformation gradient in the literature, see e.g. [Auricchio et al., 2008]. However, the electric potential φ as well as the strain-like electrical variable ε are still not given for the beam.

Mechanical kinematics in the beam
By setting the origin of the global Cartesian coordinate system O to one end of the beam as shown in Fig. 1, the location of the node ϕ(s, 0) of the straight beam can be rewritten in terms of the director d 3 (s, 0) and the arc length s Thus the components X i of the vector X in the initial configuration can be computed by projecting the position vector to the directors By applying Eq. (13) and Eq. (33), the deformation gradient at a point (X 1 , X 2 , s) in the beam can be written as, see [Auricchio et al., 2008], with the rotation tensor The derivatives in Eq. (34) can be written in terms of the beam strain measures, such as ∂d k (s,t) ∂s is related to the strain κ for bending and torsion by and ∂ϕ(s,t) ∂s is related to the strain γ for shear and elongation by Consequently, the deformation gradient F in Eq. (34) can be further formulated in terms of the beam strain measures as The variable a can be formulated in the reference configuration with the beam strain measures Γ and K The determinant of the deformation gradient reads Accordingly, the Green-Lagrange strain and the right Cauchy Green tensor are given by The inverse of the right Cauchy Green tensor can be formulated as The material time derivative of the location of a material point in the current configuration readṡ by which the time derivative of deformation gradient can be evaluated witḣ

Electrical kinematics in the beam
To formulate the electromechanical coupling problem, the electric potential will serve as the extra degree of freedom. In this work, the electric potential at the point (X 1 , X 2 , s) is represented by the electric potential at the beam node plus the increment from the beam node to the point (X 1 , X 2 ) on the cross φ(X 1 , X 2 , s)  13), the electric potential on the cross section is given by with φ o (s) the electric potential at the beam node, α(s) and β(s) the incremental parameters of the electric potential in the directions of d 1 and d 2 , respectively. The Eq. (13) describes the cross section as a plane. Correspondingly, Eq. (45) defines a linear distribution of electric potential on the cross section. If the electric potential is constant within the cross section (i.e. α, β = 0), the electric field only exists in the d 3 direction, which will lead to the uniaxial contraction in beam.
According to the Maxwell equations, the electric field is defined as the gradient of the electric potential φ, see Eq. (9). To compute the gradient of the electric potential for the beam, a similar approach as the deformation gradient in Eq. (34) can be adopted. By applying the electric potential in Eq. (45), the electric field at (X 1 , X 2 , s) in the beam can be computed as For the later computing of the virtual work induced by the electric field, the variation of the electric field is given by Since the electric field expressed in Eq. (47) is not a strain-like variable for beam, we need to formulate the strain-like electrical variable ε, which can be consistently conjugated with the beam electric displacement d e in Eq. (29). For this purpose, we start from the internal virtual work δW int induced by the electric field in the 3D continuum By applying the electric displacement in Eq. (28) and the variation of electric field in Eq. (48) to the internal virtual work, we obtain where the work conjugate of d e and ε can be observed with the strain-like electrical variable ε defined as

Strain energy function for the beam
The strain energy function of the beam model has been derived for the hyperelastic constitutive law, such as the widely used Saint-Venant-Kirchhoff model [Simo and Vu-Quoc, 1986]. However, the strain energy for other material behaviors are not given for the beam model, including the dielectric elasticity. In continuum electromechanics, many material models have been postulated to describe different material behaviors, such as the dielectric elastomers in [Vu et al., 2007]. To apply these material models in a beam formulation, the rewriting of the strain energy function in terms of the beam strains is required, which can be achieved by applying the kinematic relations introduced in Section 4. In this part, we firstly present the derivation of the widely used Saint-Venant-Kirchhoff model for the beam. Then, the strain energy function of the dielectric elastomer is formulated for the beam such that it can be used in the beam constitutive equations directly.

Saint-Venant-Kirchhoff model
The Saint-Venant-Kirchhoff model describes a linear relationship between the Green-Lagrange strain E and the second Piola-Kirchhoff stress S with the strain energy function in which λ and µ are the Lamé parameters. Inserting the linear part of E in Eq. (41) into the previous equation, the reduced strain energy Ω r reads with D = (λ + µ)d 3 (s, 0) ⊗ d 3 (s, 0) + µI.
Integrating Ω r over beam cross section, the reduced strain energy for the beam is obtained as where the material tangents D N and D K are given by with the first moments of area J ii = Σ X i X i dA, the product moments of area J 12 = J 21 = − Σ X 1 X 2 dA and the polar moment of area J p = Σ X 2 1 + X 2 2 dA. It can be observed that the widely used beam constitutive model postulated by [Simo and Vu-Quoc, 1986] and [Simo and Vu-Quoc, 1991] are recovered when λ + 2µ = E. According to the formula of Lamé parameters, this formulation assumes the limit of Poisson's ν → 0, i.e. no lateral contraction. However, the physics is incorrect and it leads to the incorrect material moduli. To overcome this problem, the uniaxial stress assumption can be additionally included within a mix formulation such as the Hellinger-Reissner principle. An alternative solution is to apply a very small Poisson's ratio ν in the simulation, by which a relatively fair comparison between the beam model and the 3D finite element model can be achieved as shown in the numerical examples of this work. Since the reduced beam strain energy Ω r b is derived by only inserting the linear part of the Green-Lagrange strain E, the derived beam force will be just the linear part of the full one. To explain it, we start from the current force n of the beam computed from the first Piola-Kirchhoff stress P It can be observed that force N is composed of two parts, one depending only on the stress and another also depending on strain variable a r . It has been proven that N is work conjugated with the beam strain Γ, see the details in [Auricchio et al., 2008]. When the reduced strain energy in Eq. (54) is applied, the derived force N r equals to the integration of the second Piola-Kirchhoff stress S lin computed with the linear strain E lin which can be seen as the linear part of the full force N in Eq. (57).

Extended Neo-Hookean model for DEA
To model the DEA, the material model of the dielectric elasticity applied in [Schlögl and Leyendecker, 2016b] for the finite element simulation is applied to the beam model in this work, where the strain energy density is given by with ε 0 the vacuum permittivity, c 1 and c 2 the electrical parameters. It can be observed that the strain energy is composed of three parts, the Neo-Hookean part referring to the pure elastic behavior, the polarization part referring to the polarization in the condensed matter and the free space part referring to the effect in vacuum. The last two terms in Eq. (59) characterize the electromechanical coupling. Apart from the dielectric elasticity, the viscoelastic effect in the dielectric material is accounted for by means of the first Piola-Kirchhoff stress P vis , see the Kelvin-Voigt model in [Wriggers, 2008], in which η is the damping parameter. The strain energy function for the beam corresponding to the continuum model in Eq. (59) can be derived with the same procedure as the Saint-Venant-Kirchhoff model. By neglecting the free space term in Eq. (59), the strain energy function for beam is obtained by integrating Ω(C, E e ) over the cross section where the integration can be evaluated with the numerical approach as well as the analytical approach.
As the analytical approach, the beam strain energy function Ω b is explicitly formulated in the Appendix.

Discrete Euler-Lagrange equations
In this work, the electromechanically coupled beam dynamics is approximated within the constrained discrete variational scheme with the null space projection. The Lagrange-d'Alembert principle for the constrained system reads where q is the configuration, L(q,q) is the Lagrangian, g represents holonomic constraints, λ is the Lagrangian multiplier and f ext (t) is the external force. By considering the electrical effect in geometrically exact beam, the electric potential φ o and the incremental variables (α, β) in Eq. (45) are treated as the electrical degrees of freedom φ = φ o α β such that the configuration of the beam model is extended to According to the kinematic assumptions in geometrically exact beams, the directors have to fulfill the orthogonal constraints The continuous Lagrangian contains the difference between the kinetic energy T (q) and the internal potential energy V (q) Since the electrical variables do not contribute to the kinetic energy, the kinetic energy for geometrically exact beams is computed as where A ρ is the mass density per reference arc-length and M i ρ are the principle mass moments of inertia of cross section. In accordance with the configuration defined in Eq. (63), the component of the consistent mass matrix corresponding to the electrical degree of freedom φ will be zero.
For the coupled hyperelastic material in DEA, the internal potential energy is computed by an integration of the beam strain energy density Ω b in Eq. (61) over the beam center line The external force f ext contains all non-conservative forces, such as the viscoelastic effect in this work. Based on the Kelvin-Voigt model in Eq. (60), the non-conservative work contributed by the viscoelastic effect is given by In this case, the external force corresponding to the viscoelastic effect can be formulated as The beam is first spatially discretized with the 1D finite elements. Then the variational integration scheme, see e.g. [Leyendecker et al., 2008], is applied to temporally discretize the action of the dynamic system, by which the good long term energy behavior can be obtained. In the variational integration scheme, the action integral within the time interval (t n , t n+1 ) is approximated with the discrete Lagrangian L d as where the discrete Lagrangian L d is computed by applying the finite difference approximation to the velocityq and the midpoint rule to the configuration q, i.e.
After the temporal discretisation, the discrete Euler-Lagrange equations can be obtained by taking the variation of the discrete action and requiring stationarity. To eliminate the constraint forces λ from the system, see e.g. [Betsch and Leyendecker, 2006], the nodal reparametrisation q n+1 = F d (u n+1 , q n ) and the discrete null space matrix P d are applied to the discrete Euler-Lagrange equations leading to where u n+1 is the generalized configuration acting as the unknown variable, f ext− n and f ext+ n are the discrete generalized external forces evaluated as

Null space matrix and parametrization of rotations
The internal null space matrix P int can be found by expressing the redundant velocityq ∈ R 15 in terms of the generalised velocity t ∈ R 9q where the generalized velocity is composed of the translational velocityφ, the angular velocity ω and the velocity of electric potentialφ, i.e. t = φ ωφ T . The corresponding internal null space matrix at time t n is written as whered i,n denotes the skew-symmetric matrix corresponding to the director vector d i,n at t n and I is the 3 by 3 identity matrix. For a multibody dynamic system composed of flexible beam actuators, rigid bodies, joints and constraints, the null space matrix can be designed by considering the electric potential as extra degree of freedom as well.
To solve the discrete Euler-Lagrange equations efficiently, the system can be reduced further into the minimal possible dimensions by use of the nodal reparametrisation, which can be achieved by introducing a rotation matrix R(θ) parametrised in terms of the rotational variable θ. The rotation matrix can be chosen as the exponential map, see e.g. [Marsden and Ratiu, 2013], The generalized configuration of the electromechanically coupled beam is specified by with u ϕ , θ and v characterizing the incremental displacement, the incremental rotation and the incremental electric potential, respectively. In this case, the nodal configuration for the next time step can be updated as

Tangent matrix
By means of the nodal reparametrisation, the unknowns of the equations system in Eq. (72) is changed from q n+1 to the generalized variables u n+1 . The nonlinear equation system is solved by use of the Newton-Rapson algorithm with the tangent matrix at iteration i in which R L (q n+1 ) is the residual of the discrete Euler-Lagrange equation. When boundary conditions are imposed on some degrees of freedom of a beam node, the corresponding components in null space matrix will be zero. By crossing out the rows involved in boundary conditions, the non-singular tangent matrix can be obtained. Accordingly, the boundary conditions can be imposed by setting specific incremental values in the update of the nodal configuration in Eq. (78). In this work, the residual vector and the tangent matrix are derived by using the automatic differentiation tool CasADi [Andersson et al., 2019].

Legendre transformation for energy evaluation and system initialization
Since the velocity involved in the kinetic energy is approximated with the finite difference scheme, it does not fulfill the hidden constraints and their time derivatives exactly. An alternative energy formulation is the discrete Hamiltonian evaluated with the momentum, which can be obtained via the discrete Legendre transformation, see e.g. [Leyendecker et al., 2008]. Since the consistent mass matrix corresponding to the kinetic energy in Eq. (66) is singular, it can not be applied to compute the Hamiltonian directly where the inverse of the mass matrix is required. Due to the fact that the kinetic energy is independent of the rate of d 3 and the rate of electric potential, a reduced non-singular mass matrixM can be defined at node s asM with I the 3 × 3 identity matrix. Therefore, the discrete Hamiltonian can be evaluated with the reduced mass matrixM and its conjugated reduced momentum To eliminate the Lagrangian multiplier term from the reduced momentump, the projected discrete Legendre transformation can be applied. In this case, the projected reduced momentum can be computed as where Qp− n is the reduced momentum at time t n ,Q is the projection matrix andq n is the vector containing the first 9 elements of q n . Corresponding to the reduced momentum Qp− n and the mass matrixM, the projection matrix is computed in the reduced form with the reduced internal constraint Jacobi matrixḠ 3×9 given bȳ The discrete Hamiltonian H d can be computed as Furthermore, the discrete Legendre transformation is applied to initialize the system at time step t 0 . In the discrete Euler-Lagrange equation (72), the unknown value of q −1 is required to solve q 1 . To avoid the computation of q −1 , the discrete momentum obtained from the Legendre transformation is initialized with an initial momentum p − 0 = p(0), which leads to the equation of motion at t 0 In this case, q 1 can be computed given q 0 = q(0) and p(0).

Numerical examples
The material parameters of the dielectric elastomer applied in this work are shown in Table 1. Instead of using the Lamé parameter µ = 0.233M pa in [Schlögl and Leyendecker, 2016b] for the nearly incompressible case, the Lamé parameter µ is set to be 233M pa in this work for the sake of a smaller Poisson's ratio. According to the discussions in Section 5.1, the smaller Poisson's ratio allows the better comparison of 3D FEM model with the beam model with assumption of no lateral contractions. The geometry of the beam with length l and square cross section b × b is shown in Fig. 3. The beam is fixed at one end. The electrical boundary condition φ i is imposed on node i such that different deformations on beam can be generated, including uniaxial contraction, shear, bending and torsion. The electrical boundary conditions φ i and the beam size for generating different beam deformations in this work are given in Table 2.

Uniaxial contraction
To generate the uniaxial contraction in the beam, the uniform electric potential is applied on the cross section, i.e. α, β = 0. As shown in Table 2, electric potential φ 1 = 0V is applied to the node at the fixed end and electric potential φ 6 = 2 × 10 4 V is applied to the top beam node. The motion of the beam discretized with N e = 5 elements, i.e. 6 nodes, is computed with the time step of 1 × 10 −4 ms. Due to the contractive forces from the dielectric effect, the vibration of the beam is induced as shown in Fig. 4, where the position X z of the beam node at the free end is depicted. In Fig. 4(a), the beam strain energy Ω b derived in the Appendix is compared with that evaluated by the numerical integration over cross section. The derived strain energy shows good agreement with the numerical model in Fig.  4(a) and will be applied in the following parts. In Fig. 4(b), the damping effect with different values of the viscosity parameter η is depicted. Without viscosity i.e. η = 0, the vibration of the beam continues over time. With the increase of viscoelastic effect, the increase of damping effect on the beam can be observed. By setting η = 0.5, the viscoelastic beam is gradually approaching a static state. To evaluate the energy property of the beam, the discrete Hamiltonian H d computed with momentum is plotted in Fig. 5. As shown in Fig. 5(a), the total energy oscillates around a value when there is no viscosity i.e. η = 0. By introducing the viscosity, the total energy decreases at the beginning stage and then preserves a smaller value, where the decrease of total energy is induced by the decrease of kinetic energy. Additionally, with the increase of the damping parameter η, the total energy decreases faster. To evaluate the energy behavior for η = 0 further, the kinetic energy T d and the potential energy V d corresponding to the total energy H d are plotted in Fig. 5(b). It can be observed that the frequency of the kinetic energy is twice as much as the frequency of the displacement in Fig. 4, which is reasonable by considering the velocity changes. However, the potential energy shares the same frequency with the kinetic energy, which is different from the pure elastic case where the elastic potential energy shares the same frequency with displacement. The reason can be attributed to the fact that the potential energy of the dielectric elastomer holding its maximum value in the undeformed state decreases to a minimum value in an intermediate compressed state, and then increases with the further compression, see the first V-shape oscillation of V d in Fig. 5(b).
The magnitude of contractions of the beam depends on the magnitude of electric potential applied on it. With the increase of electric potential, the increase of contraction can be observed as shown in the upper graph of Fig. 6. To validate the beam model, the contractions of the beam model are compared with the results of the 3D FEM model as shown in the lower graph of Fig. 6. In the 3D FEM model, the beam structure is discretized with 5 elements in longitude direction whereas with 2 × 2 elements for the cross section, where the same material model and size as the beam model are applied. To prevent the deformation of cross section, the horizontal movement of the finite element nodes is fixed to zero.
The magnitude of contractions computed by these two approaches is compared in Fig. 7(a). It  can be observed that the result of beam model agrees well with the 3D finite element model. In the beam model, 15 degrees of freedom are attached on each node. In the 3D FEM model, 16 degrees of freedom are required at each beam cross section when the cross section is discretized with one element. However, the number of degree of freedom in the 3D FEM model increases with the increase of number of elements on the cross section as shown in Fig. 7(b). It can be observed that, for the same number of elements in the longitudinal direction, less degrees of freedom are required in the beam model.

Shear
Apart form the uniaxial contraction, the electric potential formulated in this work allows more complex beam deformations by applying a non-uniform electric potential on the cross sections. To generate shear in the beam, the electric potential φ i = iφ o + αX 1 + βX 2 is applied to node i as shown in Table 2. The electric potential on the deformed beam at time t = 0.2ms is shown in Fig. 8(a), where the beam is discretized with 41 nodes and the damping parameter η is set to be 0.5. By applying the same electrical boundary conditions, the deformation is computed by the 3D FEM model as shown in Fig. 8(b). The displacement of the node at the free end of the beam is depicted in Fig. 9, where the 3D FEM model is discretized with 1 × 1 × 40 and 2 × 2 × 40 elements, respectively. It can be observed that the beam model is close to the 3D FEM model. The deviation can be attributed to the small deformation of the cross section in 3D FEM and the approximation in deriving strain energy Ω b , see the Appendix.

Bending
As shown in Table 2, the bending of the beam can be obtained by alternately applying zero and the non-uniform electric potential φ = 2 × 10 2 + 2 × 10 2 X 1 on the cross sections of beam. In this case, the non-uniform contraction of each beam element leads to the overall bending of beam. The deformed state of the beam model at time t = 0.1ms is compared with the 3D FEM model in Fig. 10. The displacement of the top beam node in x-direction is depicted in Fig. 11, where the agreement between the beam model and the 3D FEM model can be observed.

Torsion
The last example is devoted to the torsion of the beam where the torque around the beam axis has to be generated by applying the electric potential. In this case, the spiral distribution of electric potential can be applied with the electric potential at node i given by The torsion of the beam model at time t = 0.1ms is compared with the 3D FEM model in Fig. 12, where the beam is discretized with 80 elements in the longitudinal direction in both models.
The rotation angle of the beam node at the free end is shown in Fig. 13, where the beam is discretized with 20, 40 and 80 elements in the longitudinal direction, respectively. It can be observed that the

Conclusion
In this paper, an electromechanically coupled viscoelastic beam model is developed. Based on the governing equations, the kinematics as well as the strain energy functions in continuum electromechanics, their counterparts in Cosserat beam are formulated consistently. Especially, the proposed formulation of the electric potential allows for all types of the dielectric induced deformations in the beam, such as contraction, shear, bending and torsion. In the uniaxial contractions of the beam, the oscillation of the beam is induced by the contractive electric forces, where the damping effect in the motion of the beam node is observed after introducing the viscoelastic effect. Additionally, the damping behavior of the total energy is observed as well, where the discrete Hamiltonian is evaluated. Another interesting point is that the oscillation frequency of the potential energy in the dielectric elastomer is twice as much as that of the displacement, which can be attributed to the fact that the potential energy is in its maximum value in the undeformed state of the charged dielectric elastomer. The simulation result also shows that the beam model agrees well with the 3D FEM model in contraction, shear, bending and torsion, however less degrees of freedom are required in the beam model. with and where the area moments are given by J 1111 = Σ X 4 1 dA, J 2222 = Σ X 4 2 dA, J 1122 = Σ X 2 1 X 2 2 dA, J 11 = Σ X 2 1 dA, J 22 = Σ X 2 2 dA, other area moments are zero due to the symmetry of the cross section in this work.
In the above formulations, (·) i (i = 1, 2, 3) denote the components of the vector (·), such as e 1 is the first component of the vector e and d 0 12 is the second component of the vector d 1 (s, 0).