Energy-based nonlinear dynamical modeling of dielectric elastomer transducer systems suspended by elastic structures

This paper presents the Hamilton principle approach to model, design and control mechatronic systems using dielectric elastomer transducers (DET) suspended with elastic structures. An overall dynamical modeling approach for dielectric elastomer-based actuators is presented, taking into account the dynamical effects, e.g., electrical input quantities, inertia, viscous effects, and the nonlinear behavior of DETs and elastic structures. Energy-based techniques are used to obtain a coherent modeling of the electrical and mechanical domains. Based on the variational principle and using the Rayleigh–Ritz method to approximate the field variable, a nonlinear state space model is derived considering various geometric deformations and boundary conditions. The presented approach leads to a set of ordinary differential equations that can be used for control and engineering applications. The proposed method is finally applied to a multilayer DET coupled with a nonlinear buckled beam structure and analyzed based on analytical considerations and numerical simulations.


Introduction
The emerging field of soft robotics offers the prospect of applying soft actuators as artificial muscles in robots, replacing traditional actuators based on hard materials such as electric motors and piezoceramic actuators. Dielectric elastomer transducers (DET) are one class of soft actuators that can deform in response to voltage and substitute biological muscles in the aspects of large deformation, high energy density, and fast response [1][2][3]. All that is required is an elastic and soft dielectric surrounded by two flexible electrodes with layer thicknesses in the lower micrometer range and electrical field strengths up to E = 60 V μm , as shown in Fig. 1. This enables low weight and, due to the small dimensions, high energy density. The DET's operation can be described as a deformable capacitor, and its principle of operation is based on electrostatic pressure, being able to offer a large amount of deformation by applying sufficient voltage [4]. Due to their deformation under the influence of an electric field, they are suitable for certain actuator applications, e. g. proportionally controlled valves, haptic displays, and soft robotics applications [5,6].
Many research groups have recently focused on the modeling and design of soft robotics actuators based on dielectric elastomer transducers (DET) and are promising in the emerging field of soft robotics due to the advantages of light weight, high deformation, high energy density, inherent soft nature, shape control and selfsensing capabilities [7,8]. For example, humanoid robots have been developed and consider the antagonistic arrangement of DETs to mimic the function of human muscles [9,10]. DET-driven walking robots have been designed using stacked DETs [11]. A soft gripper based on bistable DET has been designed that is capable of grasping various objects with various sizes and shapes [12]. Thus, DETs have been coupled with various structures and mechanisms to achieve certain mechanical characteristics such as bistability, amplified stroke, minimum energy structures, and low energy consumption.
There are various methods and approaches for modeling DET systems combined with different structures. Models have already been developed and published describing the static [13,14] as well as transient behavior [15,16] of either the DE material in general or a certain transducer design. Advanced dynamic nonlinear continuum models have been developed taking into account the electromechanical coupling of the DET, which leads to a set of partial differential equations, and therefore, they are solved numerically using the finite element method [17,18]. This novel modeling approach is applicable for a detailed analysis and design of the transducer's coupled behavior, although it is difficult to use in control applications. For various types of DE transducers, different control-orientated modeling approaches have been previously presented by describing the DET as a lumped parameter system or using energy methods. A positioning control system based on a DE membrane [19] and a self-sensing control for a DE stack transducer [20] have been developed describing the DE actuator by a lumped parameter system. A soft robotic module actuated by rolled DET and a multi degree-of-freedom soft robot based on a double cone DET have been modeled using the Lagrange equations for a rigid body considering the DET as an external force acting on the system [21]. In addition, a thermodynamic model has been developed to describe the viscoelastic, electromechanical and dynamic response of a DET [22]. In this context, the dynamic modeling of DETs coupled with elastic structures has rarely been studied and mostly replaced by concentrated elements to model the dynamic behavior of DET systems, which have ultimately been used in design and control applications.
Within this contribution, we present an energy-based nonlinear modeling approach for DET systems combined with a flexible structure using the Hamilton principle for electromechanical systems [23][24][25]. The interactions between DE transducers and flexible structures connected with mechanical interfaces are investigated, taking into account the nonlinear and dynamic behavior of both the deformation of the elastic structure and the DETs. The approach comprises the electrical domain of the DET and is combined with the mechanical model by considering a lumped electrical circuit. In addition, multiple DETs are coupled with the structure, and different measurement points are concluded in the model, resulting in a MIMO System. The deformation of the elastic structures generally depends on the coordinate system and leads to a set of partial differential equations. In order to exclude the coordinate dependencies, the field variables are discretized with assumed elastic field shapes using the Rayleigh-Ritz formulation, where the geometric deformation of the soft actuator is guessed by polynomial functions. Thus, we obtain a state space model for the coupled DET system that reflects the actuator dynamics. As an predestinated example to demonstrate the modeling approach, a state space is developed for a buckled beam actuated by a multilayer DET.
Therefore, the remainder of this paper is structured as follows. In Sect. 2, we will present the Lagrange function, define the conservative energies and the corresponding virtual works of the mechanical and electrical domains for the DET system. The energy and work expression are approximated using the Rayleigh-Ritz method. After defining the generalized coordinates of the electromechanical system, the Hamilton principle is applied to derive the governing Lagrange equations, which can be easily transformed into a nonlinear state space model. The method is finally applied to a buckled beam structure actuated by a multilayer DET in Sect. 3. The novel concept uses elastic and dynamic instabilities of the buckled beam, while the DET provides the actuation. The actuation characteristics are investigated in Sect. 4 based on analytical solutions and numerical simulations, before finally Sect. 5 summarizes the findings.

Energy-based modeling for DET systems
A sophisticated approach for modeling the DET system is offered by the Lagrange formalism, also called the Hamilton principle, based on the calculus of variations. By using the generalized energy and power variables, the dynamic equations for mechatronic systems can be obtained. One important modeling decision is the choice of generalized physical coordinates. In the case of mechanical systems, the displacements q m are commonly chosen as the generalized coordinates (for the continuous system, the coefficients of the Ritz approach are used), and for the electrical system we choose the charge q p of the dielectric polymer [23]. Applying the Hamilton principle and variational calculus, the resulting system of equations leads to a set of nonlinear differential equations in terms of the generalized minimal coordinates q m and q p . Table1 provides the definitions of variables and quantities used in this Section.
The following steps are carried out to derive the dynamic model of a DET system: 1. Definition of the Lagrangian and corresponding conservative energy terms 2. Definition of the nonconservative virtual works of the mechanical and electrical domain 3. Field variables approximation using Rayleigh-Ritz method and definition of generalized coordinates for the electromechanical system 4. Applying Hamilton's principle and formulation of the nonlinear Lagrange equations

Lagrange functional
In the following, we will derive the electromechanical model of DETs coupled with an elastic flexible structure considering the system presented in Fig. 2. The DET system comprises an elastic structure Ω, DE materials Ω DE , and concentrated mass, spring and damper elements. External forces and electrical voltages through the DEs can act on the system. Furthermore, multiple DETs are coupled with elastic structures with voltage inputs u, and different deflection points y are concluded in the model, resulting in a MIMO System. Finally, boundary conditions are also defined, which will be later considered by the Rayleigh-Ritz approximation.
The Lagrangian L for a DET system is given by: where T , U s , and U DE are the kinetic energy, the elastic structure energy (in addition nonlinear spring elements), and free energy of the DET, respectively. The kinetic energy comprises the energies of the concentrated mass elements, the elastic structures and DETs: Next, the energy of an elastic structure is given by with the elasticity tensor C i jkl and the strain tensor ε i j depending on the displacement field u 1 , u 2 , and u 3 , considering geometric nonlinearities of the elastic body. The Einstein summation convention is used to calculate the first term in U s and the strain tensor. The components of the elastic strain energy can be simplified and adapted to the corresponding model depending on the actuator application. Therefore, Eqs. (3) and (4) can be simplified for well-known structures such as bar, beam, plate, and shell elements as applied in Sect. 3 for a buckled beam model. The second term in Eq. (3) corresponds to the external nonlinear spring elements as illustrated in Fig. 2. Finally, the free energy of the DET is introduced: which comprises a hyperelastic material model described by the Helmholtz free energy function ψ and the electrical field energy conserved in the DET per unit reference volume. The strain energy function is expressed in terms of the principle invariants of the right Cauchy-Green tensor C [26] by and it is valid for isotropic hyperelastic materials, where the eigenvalues of C are the squares of the principle stretches λ 2 x , λ 2 y , λ 2 z . Thus, we may consider ψ(λ x , λ y , λ z ) as a function of the principle stretches. Generally used hyperelastic material models are Neo-Hooke, Mooney-Rivlin, Yeoh, Gent and Ogden [27][28][29]. A uniform distribution is considered for the charges q p,i over the area of the DET so that the true electric displacement field distribution can be written as where A i and A 0,i are the area of the i-th DET in the deformed and undeformed state, respectively Altogether, Eq. (5) can be rewritten as In this paper, we will assume homogenous constant stretches, so that the integral in Eq. (8) can be replaced by the DET volume V DE . The well-known electrostatic pressure for an incompressible material σ e = ε 0 ε r E 2 can be recovered considering the energy of the DET shown in Fig. 1. The stored electric field energy is given by assuming an incompressible material λ x λ y λ z = 1. The electrostatic force in the thickness direction F e can be calculated as Thus, the electrostatic pressure can be finally obtained using the electric field relation D = ε 0 ε r E:

Nonconservative virtual works
The mechanical virtual work acting on the DET system is given by which comprises the external forces, the damping elements, the inner dissipation of the structure, and viscose losses of DETs as depicted in Fig. 3.
The Einstein summation convention is also used to calculate the dissipation of the structure Ω and can be simplified for the mentioned well-known elements. Within this contribution, the virtual work density δw v,DE of a DET is described by the viscosity η v and the Maxwell element in parallel comprising a spring E M and dashpot η M (creeping material) as experimentally proven in [15]: where ε E , ε η , and ε z = λ z − 1 are the strains of the Maxwell spring element, the dashpot, and the strain in thickness direction. Due to the connection in series of E M and η M , the tension is the same for both elements, and the sum of the strains is described by the strain of the DET in z-direction, Inserting Eqs. (14) and (15) into Eq. (13), the virtual work density δw v,DE depends only on ε E and ε z , The Maxwell strain ε E will be defined as an additional generalized coordinate modeling the creep behavior of the DET.
Next, we will describe the electrical virtual work for the lumped circuit as illustrated in Fig. 3. According to the simplified approach in [30], the electrical model comprises a series resistance R s representing losses of the electrode, a parallel resistance R p of the non-ideally insulating polymer, and its capacitance C p , which is already considered in the energy function U DE . Beside material and design parameters of the DET, these electrical parameters generally depend on the current deformation state. The resistances R s and R p are assumed to be constant and should not depend on the deformation. The electrical virtual work for n p DET elements under consideration of the input voltage v DE is described by According to Kirchhoff's current law, we obtain: The current relation (18) is inserted into Eq. (17), resulting in with the generalized voltages

Rayleigh-Ritz method
The Rayleigh-Ritz method is an approximation enabling a transformation of a partial differential equation (PDE) into a set of ordinary differential equations (ODEs). It allows us to represent a continuous system by a discrete approximation, which is expected to approximate the behavior of the continuous system. To achieve this, it is assumed that all distributed field variables (u 1 , u 2 , u 3 ) can be approximated by a general function where φ j is a set of assumed modes, which are continuous. The coefficients a j (t) are reduced by applying the geometric boundary conditions, which leads to a reduced function approximation, The n r functions of time q m (t) are defined as the generalized coordinates of the continuous mechanical system, and ψ m are the reduced modes. If the set of assumed modes is complete, the approximation converges toward the exact solution as their number n r increases. In this paper, we will approximate the modes ψ m (x, y, z) with power series transforming the problem into a finite number n r degrees of freedom. The geometric deformation and kinematic assumptions of actuators are usually known for many applications so that we can approximate the field variables by a few sets of generalized coordinates. By the volume integral in the Lagrangian formulation, the field dependency of the distributed variables can be eliminated. This yields an ordinary differential equation in time describing the dynamical model by a few time-dependent state variables similar to those of discrete systems. There are different discretization methods (e.g. finite elements, finite difference, etc.) to derive an ODE system from a PDE and apply this for control systems. However, this leads to a very large ODE order and considers higher shape components in the system, which are mostly not excited by the actuation of the DE.
Using the Ritz approach, the number of states can be significantly reduced by selecting only the most significant excitable modes of the solution and at the same time including the boundary conditions. A requirement for this is obviously the prior knowledge of the solution form or rather the geometric deformation of the structure, which is usually known for actuator applications. Thus, the electromechanical system can be fully characterized by n r mechanical variables q m (m = 1, . . . , n r ), n p Maxwell strains ε E, j , and n p electric charges q p, j j = 1, . . . , n p introducing the state vector of generalized coordinates of the mechanical domain, and the state vector of generalized charges for the electrical domain Considering Eqs. (23) and (24), the variables u 1 , λ y , λ z and their time derivatives that appear in the energy and work expressions can be described as functions of q m andq m using the kinematic assumptions and constraints of the overall system. Hence, the Lagrangian functional in Eq. (1) is formulated as a function of q m ,q m , and q p , Furthermore, the virtual works of the mechanical and electrical part are also described by the generalized coordinates and charges. Therefore, the virtual position vectors δx F,k , coordinates δx d,l , strain δε z , and strain tensor δε kl are described by the virtual generalized coordinates δq m ,  (12), with the generalized forces where the latter two expressions (28.4) and (28.5) are determined using Eq. (16). Finally, the virtual work in Eq. (27) can be rewritten in vector form as follows: with the general virtual force vectors Q F = Q F,1 , . . . , Q F,n r , The electrical virtual work in Eq. (19) is independent of q m and is also presented in vector form, with the generalized voltage vectors It should be noted that the structure of the equations remains the same if one uses the Ritz approach since the shape function and its derivatives appear as integrals over the volume and generate constant values. Therefore, the structure of the equations only changes by the corresponding order of the Ritz approach or number of independent generalized coordinates. From the perspective of systems theory, unknown system parameters can be determined by an identification procedure without selecting a specific mode or shape of the solution and hence can be determined by experimental data.

Hamilton principle
The derivation of the electromechanical system can be obtained from a variational principle applied to the action of the Lagrange energy function considering the electrical and mechanical virtual works and the Rayleigh-Ritz approximation of the field quantities [31]. Thus, the Hamilton principle for the coupled electromechanical system is given by Inserting Eq. (25), (29), and (31) into Eq. (33), we obtain Performing variation calculus to Eq.(34), we finally obtain the equation of the coupled electro mechanical DET system (Lagrange equations), with n = 2 n q + n p states, p = n p inputs and q outputs, this can be easily transformed into a state space model of the formẋ = f (x, u) , y = g (x, u) .
Equation (37) can be finally used to design and control a complex DET actuator system suspended with flexible elastic structures. The proposed method includes various nonlinear and dynamical effects for both a DET and elastic structure, which leads to a complex high order nonlinear system. The model can be reduced, if for example we assume that the dynamic behavior of the elastic structure dominates the dynamic behavior of DETs. Thus, the creep losses of the DET can be neglected, and therefore, the Maxwell strain vector q ε can be eliminated. In this context, when modeling DET systems, not all energy or work terms have to be considered. Depending on the application, some physical effects can be neglected in order to generate simplified reasonable models that are suitable for further investigations avoiding complex incomprehensible systems.

Example: modeling DET coupled with a flexible buckled beam structure
In this Section, the general method derived in the previous Section is applied to a DET, coupled with a buckled beam as depicted in Fig. 4. The electromechanical system comprises a DET with n layers, two symmetric bending beams, and a mechanical interface connecting the elastic structure with the DET. Material and geometric parameters are listed in Table 2. The DET is stretched in such a way that the beam reaches a buckled state with a maximum deflection q. By applying a voltage to the DET, the actuator moves in vertical direction. Thus, the voltage v DE is defined as the input and the vertical displacement y as the output, resulting in a SISO System. Furthermore, two different mechanical interface configurations are examined and compared with each other, which can be easily considered by various boundary conditions. We will consider a simply supported beam and a fixed beam evaluated for the DE actuator system. In the first step, before starting to derive the model equations, we will introduce the kinematic assumptions applied to the system. Energy and virtual work expressions are formulated for the DET material and beam elements considering the mechanical and electrical behavior of the model. The Rayleigh-Ritz approximation is performed, which accounts for different boundary conditions and in addition defines the generalized coordinate.  Finally, the Hamilton principle is applied to the electromechanical system, resulting in a nonlinear ODE system that can be easily transformed into a state space model and utilized for further investigation. Geometric nonlinearities are considered to accurately describe the instable buckling behavior of the beam elements and therefore induce multiple equilibrium solutions [32], which appear in our actuator concept.

Beam element assumption
The Euler-Bernoulli beam theory is considered, which provides a simple means of calculating the load and deflection characteristics of beams by introducing the kinematic assumption that a cross-sectional plane does not deform and thus only undergoes rigid displacements and rotation. This assumption is equivalent to neglecting the shear deformation of the material. A beam of length l, mass density ρ, and a cross-section area A is considered. The transverse displacement on it is uniform and limited to the transverse displacement in the xz-plane, which depends on the coordinate x, and time t. The axial displacement component results from the rotation of the cross section such that it remains orthogonal to the neutral axis, and depends on the coordinate x, z, and time t. With Eq. (39), large rotations can be considered, which is necessary for an accurate characterization of the buckled beam structure [33]. The linearized time derivativė is considered when the rotatory inertia of the beam cross section is included. Based on the small strain assumption, the nonlinear strain terms in Eq. (4) are neglected, Considering the assumption of the Euler-Bernoulli beam, the strain tensor is reduced to the strain in x-direction: κ(x, t) is the physical curvature of the deformed beam, and it is a function of w only. It is valid for small strains and does not contain any restriction on the magnitude of rotations, which is necessary for an accurate consideration of the buckling behavior [33]. A linearized form of the curvature κ is also necessary, which will be used for the small dissipation resulting in the beam elements. Finally, for isotropic material behavior with a uniaxial stress state, the elastic strain energy simplifies to DET assumption An n-layer DET with length l 0 , area A 0 , and thickness d 0 is considered. Within this contribution, the Neo-Hooke model is used for modeling the static mechanical behavior of the DET. Due to the incompressibility of the elastomers considered, on the one hand the volume is constant. On the other hand, a planar deformation in the x-direction occurs due to the design of the DE transducer considered, i.e., the stretch λ y is equal to one, and λ x directly depends on the stretch in the z-direction λ z : The stretches are assumed to be constant and do not depend on the coordinates x DE , y DE , z DE . Finally, the kinetic energy and viscose losses are assumed to be small compared to those of the beam elements, and therefore, they are neglected. Nonetheless, the viscose losses of the DE material are considered in the calculations, which induces an additional internal state variable q ε = ε E .

Energy functions and virtual works of the coupled DET beam structure
We will start by calculating the Lagrangian of the DET system, which comprises the kinetic energy, the strain energy of the symmetric buckled beam, and the free energy of the n-layer DET.

Kinetic energy
The kinetic energy T is determined for the mechanical interface as a concentrated mass element m, and the kinetic energy for the two beam elements B 1 = B 2 = B is simplified considering Eq. (2) and the assumptions of the Euler-Bernoulli beam theory: Inserting Eqs.(40) into (48), the kinetic energy of the beam can be written depending on w: where I = A z 2 d A is the geometric moment of inertia of the cross-section area A. The first term in the integral describes the kinetic energy for the rotation of the cross sections, while the second is the translational kinetic energy.

Strain energy
Given Eq.(42), the strain energy U s of the system is written in terms of w and its derivative w , w : The expression 1 1−w 2 is typically approximated by a second-order Taylor series to simplify the calculation of a beam structure [34], resulting in where E I represents the bending stiffness of the cross-section.

Free energy of DET
Finally, the free energy U DE of the n-layer DET is given, taking into account the Neo-Hooke model and a uniform distributed charge over the DE area A 0 : Substituting the assumptions (46) and (47) into Eq. (52), the free energy is simplified: (53)

Mechanical virtual work
We will consider the dissipation of the beam and the viscose losses of the DE material as a mechanical virtual work acting on the DET system (see Eq. (12)), which can be expressed by: where η, η v , η M represent the damping coefficients of the beam and DE material (viscose losses, Maxwell element) and E M the elasticity modulus of the Maxwell element. By substituting the linearized expression (43) into the first expression of Eq. (54), the mechanical work is given depending on w(x, t), ε z = λ z − 1 and ε E :

Electrical virtual work
The electrical virtual work is simply rewritten for the single electrical circuit with one mesh, having the generalized voltages

Rayleigh-Ritz approximation and generalized coordinate
In the following, the Rayleigh-Ritz approximation is applied to the suspended DET beam system which allows us to represent the continuous system by a discrete approximation. The single field variable appearing in the equations above is w(x, t), which is approximated by the mentioned approach. Thus, the one-dimensional displacement field w(x, t) can be written as a polynomial function, As previously mentioned, we will consider two different mechanical interfaces in Fig. 4, which leads to a different approximation of w(x, t) and therefore two different sets of equations with the same structure. In order to describe the beam deflection by only one general coordinate, the simply supported beam model will be approximated by a third-order polynomial, while the fixed beam model is described by a fourth-order polynomial, Both defined functions (59) and (60) result in one degree of freedom models by satisfying the following geometric boundary conditions for both cases as follows: • simply supported beam model of Fig. 4: • fixed beam model of Fig. 4: Inserting the geometric boundary conditions (61) and (62) into the corresponding approximations (59) and (60), the field variable can be simplified by only one generalized coordinate q(t) regarding the simply supported beam, and the fixed beam, Thus, all kinematic variables appearing in the energy and work terms can be described by just one generalized mechanical coordinate q (n r = 1), and the corresponding shape function S(x) (S(x) represents S s (x) or S f (x)), leading to a one degree of freedom motion as aimed before. The vertical deflection y that appears in the kinetic energy of the mass m, can be written depending on w and consequently as a function of q using the Rayleigh-Ritz approximation [35]. As we will see later, Eq. (65) represents the measured output of the system, and it is a nonlinear (quadratic) function of q. The assumed constant stretch λ x is also introduced as a function of q describing the ratio between the actual deformed length l − y and the initial length of the DET l 0 by: Finally, the virtual mechanical work in Eq. (55) is rewritten depending on δq and δε E according to: where Q s and Q v comprises the generalized damping force of the beam and DE material, respectively. The generalized force Q ε represents the viscose losses caused by the Maxwell element, and it adds a new internal generalized coordinate ε E to the system. In summary, the Lagrange function and the virtual works of the holistic electromechanical model are represented by the generalized coordinate q, the internal coordinate ε E , and the generalized charge q p . The final step is to derive the dynamical model of the electromechanical system using the Hamilton principle, which leads to a set of nonlinear differential equations (Lagrange equations).

Hamilton principle and Lagrange equations
The Hamilton principle is applied to the DET beam model, resulting in the following variational equation (see Eq. (34)): where Q s , Q v , V p , and V e are determined from Eqs. (67) and (57). The Lagrange function of the electromechanical system is summarized by the following expression: using Eq. (49), (51), and (53) with It should be noted that the structure of the equation is not changed by the choice of the shape function S(x), except the constant terms in the above equation. Using Eqs. (63)-(66) in Eq.(70), the Lagrange function L q,q, q p is represented by the generalized coordinates q,q, and q p . According to Eq. (35), we finally attain the Lagrange equations of the DET beam model: where the fourth equation of (71) is substituted into the third equation to eliminate the variableq DE . The second equation represents the viscose losses of the DE due to the Maxwell element, and it is specified by the inner generalized coordinate ε E . The calculations were performed and lead to the following set of equations: of the DET for the two different boundary conditions. Considering the blue curves (v DE = 0 V), the initial vertical state of the beam equilibrium is stable when ε b < ε b,min , and it remains vertical up to ε b = ε b,min , which corresponds to point A on the equilibrium paths. When increasing this value, the initial vertical equilibrium state becomes unstable, and any arbitrarily small disturbance definitely brings both beams out of the vertical state.
Since we have a stable deflected equilibrium state in the vicinity of point A, the beam will switch to a new stable state to either the right-hand q + (convex shape) or left-hand q − (concave shape), depending on the disturbance. At this bifurcation point, the solution splits into several, here two, paths. For ε b > ε b,min , the DET beam actuator has two stable equilibrium points B and B beside the unstable point at q = 0. The stability can be analyzed considering the total potential energy of the system, which comprises the potential energies of the beams U s , the DET's energy U DE , the constant initial pre-strain x,pre − 2 with λ x,pre = l l 0 , and the driving voltage source v p q p . The potential energy depends on the coordinate q and charge q p stored as electrical energy in U DE,e (q, q p ). The energy expression in Eq. (78) is appropriate if the charge q p is controlled in the system. In our case, we adjust voltages to deform the actuator, and therefore U pot (q, q p ) is reformulated depending on the polymer voltage v p using the Legendre The energy U DE,e (q, v p ) describes the electric co-energy of the DET, and it is calculated for a homogenous field by where λ x is determined using Eq. (66) and C 0 = nε 0 ε r A 0 d 0 is the capacitance of the DET in the initial state. At steady state, the polymer voltage v p is related to the input voltage v DE considering the electric circuit in Fig. 3 for this case: The normalized energy function is illustrated in Fig. 5c and d for the three chosen pre-strains ε b . The minimum points B, B indicate the stable equilibrium points, and the maximum point at q = 0 corresponds to the unstable one A. By increasing the voltage up to v DE = 3 kV, the stable states move from point B to C and from point B to C, respectively, as presented in Fig. 5a-d. The dashed lines between C and C in Fig. 5 describe the unstable regions and therefore cannot be actuated within the voltage range 0 ≤ v DE ≤ v max . However, it is possible to jump between the two stable states if the actuation force momentum is sufficiently large, enabling a bistable operation of the DET beam system between the convex and concave shapes, and vice versa. Figure 5e and f demonstrates the vertical displacement of the DET beam actuator, which is defined as the difference between the initial vertical displacement due to the pre-strain y 0 at v DE = 0 and the actual displacement y = g(x) according to Eq. (74) due to the actuation v DE = 0 by: Hence, the DET beam actuator has two operation modes, namely in the vertical direction with the variable Δy and the horizontal direction q. Comparing different mechanical interfaces, the horizontal deflection q as well as the vertical displacement Δy of the fixed-end beam (right half of Fig. 5) are larger than those of the simply supported beam shown on the left half in Fig. 5. Furthermore, the bifurcation point is reached at a higher pre-strain because the motion is hindered by the fixed clamping. While both actuator designs strongly differ considering the horizontal displacement q at different pre-strains, the vertical deflections Δy hardly differ from each other, as shown in Fig. 5e and f.
The dynamic behavior of the actuator is studied using the step responses illustrated in Fig. 6. The simulations are performed for the aforementioned pre-strains by different beam lengths. Figure 6a and b shows the evolution of the beam deflection q with respect to time t.
The transient behavior for t < 1.5 s without input voltage appears due to the pre-strain of the DET. Oscillations increase with increased pre-strain. After a short settling time, the steady states are reached, where the actuator can be adjusted with input voltages. At time t = 1.5 s, the actuator is excited with the maximum voltage v DE = v max = 3 kV as depicted in Fig. 6a and b for the simply supported and fixed-end beam. Due to the momentum of the beam, the actuator switches from the state q + (B) to q − (C) for ε b = 0.15, 0.2 and moves toward (B) within 3 < t < 4.5 s when the input voltage is switched off. Thus, the actuator can switch between the two equilibrium states depending on the actuation, which demonstrates the bistable characteristic of the system. However, considering a large pre-strain (ε b = 0.25), the instability range between C and C is too large to cross the maximum energy at q = 0 with v DE = 3 kV, and therefore, it remains in the same path q + or needs much higher voltage momentum to jump into the concave state at C.
In order to investigate the vertical deflection of the actuator, the displacement Δy = y 0 − y is investigated, which leads to the graphs shown in Fig. 6c and d. In addition, the actuation strain is defined as where l y = l − y 0 is the length of the actuator in the initial state without actuation, as shown in Fig. 6e and f. The vertical deflection of the simply supported beam actuator is about ε y ≈ 6.5 %, while the vertical deflection of the fixed-end beam is approximately ε y ≈ 8%. It can be observed that the fixed-end beam performs much better than the simply supported beam. Nevertheless, the vertical motion is independent from the shape of the beams. Both the convex (q + ) and the concave (q − ) beam shapes result in the same values for Δy and ε y . However, compared to actuators with poor DE films the deflection of DET beam structures increases the strain ε y . Within this Section, we have demonstrated the functionality of a novel actuator concept that uses the instability of flexible structures to improve the actuator performance by significantly amplifying the deformations of DET systems. Besides the application in soft robotics, the actuator can be used for example as a valve in fluidic applications, which ensures a safe operation due to the supporting beam elements preventing fatigue of the DET. Furthermore, an energy efficient operation is possible due to the bistable characteristics of the DET beam actuator. Another example is to utilize the DET beam system as a sensor to measure pressures in fluids or perhaps forces within its stable operation region.
Further configurations of a beam-supporting DET structure are feasible offering more degrees of freedom in motion. Here, two beam elements are buckled by one DET in the concave and convex shape, respectively. For example, when two DETs are connected outside to a T-shape beam, a rotation of the structure is obtained in addition to the longitudinal deflection, also making use of buckling to amplify its deflections, as illustrated by certain configurations in Fig. 7. In general, the introduced buckling concepts can be multifariously combined to generate more complex kinematic motions, and in particular they are suited for soft robotic applications. The proposed modeling approach is predestined to analyze and optimize such MIMO-DET configurations in future work.

Conclusions
In this paper, an energy-based modeling approach has been presented to model the nonlinear dynamic behavior of soft material structures comprising dielectric elastomer transducers (DET) suspended with flexible elastic structures. For this purpose, the Hamilton principle was applied considering Lagrange energies and virtual works of both DETs and elastic structures as well as electrical and mechanical inputs. Thus, the overall electromechanical dynamics of the coupled DET structure were derived without the need to consider the inner dynamics and forces of the system. In order to obtain a set of nonlinear differential equations for the soft suspended DET structure, the Rayleigh-Ritz formulation was used to approximate the continuous field variables appearing in the system by generalized coordinates under consideration of defined boundary condition. The modeling approach is predestined for multiple-input-multiple-output (MIMO) systems. Nonlinear ODEs obtained by this approach can be easily utilized to design and control the soft transducer systems based on DETs. The method was successfully applied to a buckled beam structure actuated by a multilayer DET, resulting in complex nonlinear equations. As analyzed in detail based on analytical and numerical investigations, the soft actuator system shows structural instabilities and nonlinear characteristics of the DET. The behavior obtained significantly improves the actuated deflection of the soft actuator in two directions of motion and offers an energy efficient operation due to the bistability of the buckling beams. In the future, we will perform experimental results to prove the actuator's performance and will further extend the work by designing and optimizing more complex DET beam systems with multi-degree-of-freedom motions, e.g., for soft robotic applications.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.