Numerical study on the electromechanical behavior of dielectric elastomer with the influence of surrounding medium

ABSTRACT The considerable electric-induced shape change, together with the attributes of lightweight, high efficiency, and inexpensive cost, makes dielectric elastomer, a promising soft active material for the realization of actuators in broad applications. Although, a number of prototype devices have been demonstrated in the past few years, the further development of this technology necessitates adequate analytical and numerical tools. Especially, previous theoretical studies always neglect the influence of surrounding medium. Due to the large deformation and nonlinear equations of states involved in dielectric elastomer, finite element method (FEM) is anticipated; however, the few available formulations employ homemade codes, which are inconvenient to implement. The aim of this work is to present a numerical approach with the commercial FEM package COMSOL to investigate the nonlinear response of dielectric elastomer under electric stimulation. The influence of surrounding free space on the electric field is analyzed and the corresponding electric force is taken into account through an electric surface traction on the circumstances edge. By employing Maxwell stress tensor as actuation pressure, the mechanical and electric governing equations for dielectric elastomer are coupled, and then solved simultaneously with the Gent model of stain energy to derive the electric induced large deformation as well as the electromechanical instability. The finite element implementation presented here may provide a powerful computational tool to help design and optimize the engineering applications of dielectric elastomer.


Introduction
Since its electrically actuated large deformation was identified at the beginning of this century, dielectric elastomer has become one of the intensively researched subjects [1][2][3]. As a class of electroactive polymers (EAP), dielectric elastomer possess many other attractive properties, in the combination of lightweight, low cost, high efficient, and flexibility, showing great potential to realize complicated actuations, especially in soft CONTACT [4,5]. During the past decade, diverse prototype devices have been presented for broad applications, including artificial muscles, adaptive optics, micro vehicles as well as energy harvesting system to recycle energy from ocean waves and human activities [6][7][8][9].
To better understand the response of dielectric elastomer under electric actuation and efficiently control its output deformation, analytic model is necessary to give insight on the nonlinear electromechanical coupling problem, which can be depicted by the following simple positive feedback mechanism. The material is polarized by the electric field established between the compliant electrodes, and expands its area simultaneously under the effect of electric force. The deformation of dielectric elastomer leads to the increasing of electric field and correspondingly strengthens the electric force, whose changes in turn affect the mechanical response until a stable equilibrium status is reached. Pelring et al. [1,2,10] first introduce electroelastic pressure perpendicular to the electrodes to analyze the squeeze and elongation of incompressible elastomer. Due to the one-dimensional model, the equivalent pressure is twice the pressure presented in a rigid parallel capacitor. Later, Kofod [11] combines Maxwell stress and different hyperelastic models, including Ogden, Neo-Hookean, and Mooney-Rivlin stain energy functions in the constitutive model to explain the electric induced large deformation. In general, when the applied electric field can be treated as static, the behavior of dielectric elastomer can be modeled using the theory of nonlinear electroelasticity [12,13]. However, Helmholtz free energy of the dielectric elastomer is employed in most of the previous theoretical studies, only considering the electric field in the solid continuum [14][15][16][17][18][19]. Moreover, it is unable to derive the analytical solutions for this rather complicated problem, effective design and optimization of more complex applications demand powerful computational tools.
In the past few years, several numerical studies on the nonlinear response of dielectric elastomer using FEM have been presented. Wissler et al. [20] use the general purpose finite element program ABAQUS to calculate the time history of desired voltage corresponding to the measured radial stain. In this approach, the problem is treated as uncoupled by imposing the current thickness as kinematic boundary condition to derive the compressive pressure. Then an improved method has been proposed, in which the activation voltage is set as the input by defining the evolution of the electromechanical pressure in a Fortran code [21]. Vu et al. [22,23] present a numerical frame work to characterize the fully coupled electroelasticity, where the finite element method is used to calculate the elastic material domain and the boundary element method is used to calculate the electrostatic domain surrounding the body of interest. However, electromechanical instability, which may limit the achievable deformation, is not considered in their approach. Recently, Park et al. [24] present a three-dimensional FEM calculation to capture the electromechanical instability by including inertial effects and solving the governing equations monolithically. The abovementioned, fully coupled numerical methods are all implemented in homemade codes, which are not easy to master for the majority users. Till now, there is seldom report on the calculation of electro-elastostatics using commercial FEM package.
In this work, we present a fully coupled finite element approach using commercial multiphysics package COMSOL v5.0 to investigate the behavior of a circular dielectric elastomer under electric actuation while considering the electric field in the surrounding space. The basic governing equations of nonlinear electroelasticity, describing the coupled mechanical deformation and electric field, are recalled in the following section. Then the influence of free space surrounding the dielectric elastomer on the electric field is analyzed. By employing Maxwell stress tensor as actuation pressure, the detailed process of solving the Maxwell equation and linear momentum equation simultaneously for the dielectric elastomer with appropriate boundary condition in the FEM package is described. Based on this approach and the adopted Gent model for hyperelasticity stain energy, the electric induced large deformation of dielectric elastomer is calculated together with the electromechanical instability. The approach proposed here may provide a convenient tool for the researchers to investigate the electromechanical properties of more complex applications using dielectric elastomer in engineering.

Governing equations in nonlinear electroelasticity
Our finite element approach is based on the nonlinear electromechanical field theory established in the past decade [5,12,13]. Here, we summarize some essential questions for a dielectric elastomer undergoing large deformation when subjected to a combination of electrical and mechanical load.
In the stress-free reference configuration B 0 , the position vector of each point in the undeformed material is denoted by X. Correspondingly, the material point moves to a place with the spacial position vector x ¼ ϕðXÞ in the deformed configuration B t . By introducing the deformation gradient tensor F ¼ @x=@X, the deformation at every point of the material can be characterized. And the quasistatic linear momentum balance usually can be expressed as the following in reference to the deformed configuration: where σ T is the transpose of Cauchy stress tensor and b 0 denotes the volume force on the body. At the boundary of deformed material @B t , the continuity condition requires the Cauchy stress tensor to satisfy the standard jump condition where n is the normal unit vector pointing outward from the boundary @B t , t m t is the surface traction resulting from mechanical source and ½ ¼ ½ inside À½ outside represents the quantitative jump across a surface.
When the dielectric elastomer is immersed in an electric field, the mechanical deformation field interacts with the static electric field. As the capacitance of dielectric material is relatively tiny, it is reasonable to neglect the electrodynamic effects, which take place in a small tine range. For simplicity, the magnetic field, free currents, and free electric body charges are neither considered, and then the four Maxwell equations [25] covering the electromagnetic effect are reduced to the following two equations where the spatial electric field e is determined by Faraday's law of induction and spatial electric displacement d is determined by Gauss's law of electricity. Correspondingly, at the boundary of @B t , the jump condition for electric field and displacement can be written as where q t is the free surface charge density in reference to the deformed configuration B t . Due to the polarization of electric dipoles distributed inside the homogenous and electrically linear material, the spatial electric displacement in the dielectric elastomer could be calculated by the function of electric field e and spatial electric polarization p: where ε 0 is the permittivity of vacuum and ε r is the relative permittivity of dielectric material. Moreover, the electric field exerts electric force on the polarized body, whose mechanical effect can be described by the body force density method (including Korteweg-Helmholtz and Kelvin force density method) or surface force density method (such as Maxwell stress method) [26]. In our approach, Maxwell stress tensor σ M is directly used as surface stress and added with the elastic component σ e to form the total stress tensorσ,σ Without loss of generality, for an isotropic dielectric, in the absence magnetic field, the Maxwell stress tensor takes the form of [27] where I is the second order identity tensor, α 1 and α 2 are the electrostrictive parameters, expressing the change of permittivity resulting from material deformation parallel and perpendicular to the electric field, respectively. Substitute Equations (6) and (7) into Equations (1) and (2), we could derive the balance equation of linear momentum and the boundary condition with electrostatics considered. In summary, the governing equations of coupled electric field and hyperelastic solid along with the boundary condition in reference to the deformed configuration can be written as It should be notified that the total surface tractiont t and electric surface fluxq t contain the components from the free space surrounding the dielectric elastomer.

Energy in the free space
In most of the previous theoretical studies, the dielectric elastomer together with the applied electric field is taken as a close thermodynamic system. As a result, the total energy is stored inside the space occupied by the material, and a Helmholtz free energy which is a function of electric field and mechanical deformation at every material point can be employed to help construct the thermodynamic equilibrium. However, if a considerable part of energy is stored in the free space around the dielectric elastomer, for more accurate analysis, we should consider not only the solid continuum but also the adjacent air or vacuum domain, making the numerical simulation more cumbersome. Therefore, we firstly give an evaluation for the influence of surrounding space on the electric field by using the electrostatics module embedded in the finite element software COMSOL Multiphysics 5.0 (COMSOL AB, Stockholm, Sweden).

Electric field
As shown in Figure 1, besides the circular dielectric elastomer with 10-mm radius and 0.5-mm thickness, the vacuum domain surrounding the solid continuum is added to the geometric model by a 30-mm radius sphere and an additional spherical shell. The additional layer on the sphere periphery is assigned to an Infinite Element and its outer boundary is kept as zero charge boundary condition corresponding to the fact that there is no stationary charge in the surrounding space sufficiently far away from the body. The employed Infinite Element feature ensures accurate calculation of the electric field by modeling an infinitely large vacuum domain with finite geometry. An electric potential difference is applied on the dielectric elastomer by setting the lower surface as electrical ground whereas the voltage on the upper surface is 5000 V. Mapped meshing with four-node quadrangular elements is performed in the solid domain and infinite element while the finite vacuum domain is meshed with free quadrangular strategy. There are totally 1604 domain elements together with 207 boundary elements, and further increasing the element number has little influence on the derived electric field. The direct stationary solver is employed to calculate the steady state of our numerical model, in which the nonlinearity is automatically detected through analysis of the variables contributing to the residual Jacobian matrix and the constraint Jacobian matrix. With a relative tolerance of 10 -3 , all the calculations converges within 25 iterations. Figure 2(a) shows the distribution of the electric potential in the dielectric elastomer and the surrounding vacuum. The direction of electric field vector is also illustrated with black arrows. Interestingly, the electric potential in the space far away from the dielectric elastomer almost equals to 2500 V, half of the applied voltage. Nevertheless, the electric potential difference remains zero, coinciding with our daily life experience that it is safe for human body to keep a certain distance from the equipment operated with high voltage. And the electric potential variance along the z-axis is plotted in Figure 2(b), it is clear that the electric potential decreases along with the increasing of vertical deviation distance just above the upper surface, while this tendency reverses in the lower space. The inserted figure magnifies the linear distribution of electric potential in the dielectrics along thickness direction.
The z-axis and radial component of the electric field intensity in the dielectric elastomer and its vicinity is illustrated in Figure 3(a) and 3(b), respectively, along with the arrows for the direction of total field at each location. Inside the dielectric elastomer, the electric field is almost uniform with the intensity of 10 7 V=m, pointing from the positive electrodes to the ground. In most of the surrounding free space, the direction of electric field vector is from the location with higher electric potential to the lower potential places. Correspondingly, the electric field intensity variance along the z-axis is plotted in Figure 3(c). The field intensity in the vacuum adjacent to the dielectrics is about 2 Â 10 5 V=m, two orders smaller than the field inside the solid continuum and this value further decreases along with the increase of vertical deviation distance. However, in the space near the circumference of dielectric elastomer, the distribution of electric field is quite complex, which can be seen in Figure 4 (a), a magnification of this area. The maximum upward electric field with 4:27 Â 10 6 V=m intensity appears at the outside corner in the vacuum. It should be notified that the FEM solution runs into a geometrically singular result in this area. Though electric potential itself remains smooth and well-defined at the sharp corner, its gradient, the electric field is in theory infinite. In numerical calculation, this value will tend toward infinity as you refine the mesh. Moreover, unlike the sudden jump of electric field intensity at the upper and lower interface shown in Figure 3(a), the electric field is continuous at the circumference and the field intensity in vacuum decreases to an average level 10 5 V=m at a distance about 1-mm away from the outer edge. Figure 4(b) and 4(c) show the z-axis component of the electric field along different lines marked in Figure 4(a) by (a)-(d). The field intensity in this area is much larger than the value in the region near the top and bottom surface. From above analysis, we can conclude that the energy distributed in the free space is rather weak except that in the small area adjacent to the circumferential edge of dielectric elastomer.

Distribution of Maxwell stress tensor
Maxwell stress tensor can be derived directly in the electrostatics module, which is the outstanding advantage of this multiphysics FEM package. Figure 5  calculated radial component of Maxwell stress tensor is smaller than the theoretical prediction.
In both of these two figures, the numerical values at the sharp corners severely deviate from the expectations. Indeed, the algorithm for Maxwell stress tensor employed in COMSOL is slightly tricky, which is based on a few simplifications and may make use of the steep gradients, such as electric potential. In the former section, we have pointed out that the calculated electric field in the sharp corner is not convincible because of the geometric singularity. So as the Maxwell stress tensor, a function of the electric field. Moreover, according to the user guide of the Maxwell stress tensor in COMSOL, dense and symmetric mesh is expected at the target to get more reasonable results. From the inserted magnification in Figure 2, we may find the symmetry of employed mesh in the sharp corners is not satisfying, which may also account for the deviation. If we want to use Maxwell stress tensor to simulate the response of dielectric elastomer under electric

Numerical implementation
In the following, we will illustrate how to perform the formerly summarized formulation of coupled electrostatics and dielectric hyperelastic material in the finite element software COMSOL. For any dielectric elastomer actuator with arbitrary geometry, the governing Equation (8) can be implemented in the weak form PDE module of this FEM package by using its corresponding weak form governing equations: where η is the test function and σ e can be derived with the selected energy density function of the solid continuum. After meshing the geometric model, this fully coupled nonlinear system can be solved with a variety of direct and iterative solvers embedded in the software. As the numerical implementation of nonlinear electroelastic problems in weak form is still not user friendly, here, we introduce an alternative approach only using the pre-established physical modules solid mechanics and electrostatics.

Geometric and material parameters
Let us consider a circular dielectric elastomer actuator made from VHB4905 (3 M corporation) as shown in Figure 6, which has been comprehensively investigated by previous theoretical analysis and experiments. In the absence of applied load, the circular membrane with 0.5-mm thickness and 10-mm radius is sandwiched by two compliant electrodes. Then the thin membrane is biaxially stretched by a uniform distributed constant force P along its boundary. When a voltage is applied between the compliant electrodes, the actuator squeezes in thickness and expands in the radial and latitudinal directions. Two-dimensional axisymmetric model is used for the dielectric elastomer, in which the mechanical effects of compliant electrodes and viscoelasticity of VHB4905 have not been considered. The origin of the coordinate system is set at the center of the circular membrane. According to the reported experimental data [28], the value of relative permittivity ε r for VHB4905 is 4.7, independent of the deformation. Therefore, electrostrictive parameters α 1 and α 2 are assumed to be zero in out numerical simulation. And the remaining mechanical properties of the hyperelastic elastomer are given through the stain energy function defined by Gent model [29] W ¼ À μ 2 log 1 À where J ¼ det F, I 1 is the first principle invariant of the right Cauchy-Green deformation tensor, the macroscopic shear modulus μ is chosen to be 45kPa, j m is 80, and the initial bulk modulus K is 500MPa corresponding to the nearly incompressibility.

Modified approach using actuation pressure
Once the vacuum domain surrounding the dielectrics was considered in the numerical model, the simulation becomes rather cumbersome and time-wasting. Besides the additional elements employed for the free space, when dealing with the large deformation of hyperelastic material, the major challenge of FEM lies not only in the amount of effort spent on computing the mechanics and electric fields inside the solid continuum, but also in the effort Figure 6. Schematic of a circular dielectric elastomer with undeformed 10-mm radius and 0.5-mm thickness. Electrodes coated on the either side of the dielectric elastomer membrane are sufficiently complaint such that they do not constrain the material deformation. A uniform distributed constant force P is applied along the circular boundary and biaxially stretches the membrane before the electronic actuation.
to remesh the vacuum domain after several iterations [22,23]. The mesh will be distorted along with the expanding of dielectric elastomer, however the critical status, at which the current element mesh is not acceptable, has to be chosen manually in the calculation. As the thickness of dielectric elastomer is much smaller than its radius, considering the given analysis that the electric field intensity in the vacuum above the upper surface and below the bottom surface is tiny, the surrounding free space will not be included in our model, whereas the mechanical effects resulting from the vacuum adjacent to the circumferential edge of the dielectric elastomer should be taken into account.
Across the circumferential electric-air boundary, the electric field is continuous, so the radial component of Maxwell stress in the air does not vanish and takes the following form Correspondingly, surface traction on this interface contains the electrical component given rise by this source. Thus, the mechanical boundary conditions for the electric-air interfaces can be written as, where the t e is the surface traction resulting from the mechanical load P and σ M die is the Maxwell stress tensor in the dielectrics. Taking advantages of the nearly uniform distributed electric field in the dielectrics, we transfer σ M die from the left side of Equations (9) and (10) to the right side as part of the surface traction, then the material deformation could be simulated by this imaginary actuation pressure approach, in which the following steps are included.
(i) The circular dielectric elastomer is biaxially stretched by imposing a surface load in mechanics module. (ii) The electrostatics module is added, where the electric potential difference is applied on the deformed configuration. And the actuation pressures, functions of the calculated electric field from electrostatics module is coupled in the mechanics module as additional surface load. Figure 7(a)-(d) demonstrate the finite element results for the equilibrium deformation of a circular dielectric elastomer with an identical 4000 V actuation voltage and various prestretch conditions. For larger pre-stretch values, the electric induced area expansion is larger, which agrees well with the tendency derived from previous analytical studies. As the slightly nonuniform electric field has been neglected in our numerical model, the elastic stress and strain are generally homogeneous over the membrane. However, for the huge deformation case as shown in Figure 7(d), the radial strain in the vicinity of the outer circumference becomes imperceptible nonuniform. Figure 8 shows the relations between the applied voltage and total stretch. For directive comparison, the curve without including the influence of surrounding space is also illustrated, in which the Maxwell stress in the vacuum adjacent to the circumference boundary of dielectric elastomer is not considered. Table 1 illustrates the overestimated relative deformation without the influence of surrounding space. It can be observed that under the same actuation voltage, the radial elongation is smaller when the surrounding space has been taken account, especially at large deformation conditions. As it is wellknown in dielectric elastomer theory, the extent of electric induced inplane strain is limited by electromechanical instability and larger pre-stretch level will improve the stability, a large deformation, more than 4.5 times of its initial radials is only achieved for the case with a pre-stretch value of 2.1. For other pre-stretch conditions, the finite element calculation fails due to nonconvergence at a relative small deformation status.

Electromechanical behavior with the influence of surrounding space
The electromechanical instability of dielectric elastomer can be analyzed systematically by using the Hessian, which must be positive definite at the equilibrium state of a large deformation. According to Zhao and Suo's study [30], the Hessian can be analyzed from the curves of normalized voltageẼ ffiffiffiffiffiffiffi ε=μ p versus the normalized chargẽ D= ffiffiffiffiffiffiffi ε=μ p as plotted in Figure 9. The solid curves are derived by modifying our numerical approach with the final equilibrium status as an input Global Equation restriction. The  left-hand side of each curve corresponds to a positive definitive Hessian, while the right hand-side represents the non-positive condition, and the peak accords with the condition of detðHassianÞ. It is clear that the peak falls and finally vanishes by increasing the pre-stretch level. However, for the smaller pre-stretch case, the former calculation method with voltage as input collapses at the vicinity of the peak, the beginning of electromechanical instability.

Conclusion
In this work, an FEM computation approach for the electromechanical response of a prestretched circular VHB4905 dielectric elastomer is presented within the commercial multiphysics package COMSOL. Maxwell stress tensor as an effective actuation pressure in our implementation. Electrostatics and mechanics modules, being coupled through the applied actuation pressure as surface load, are used to solve the Maxwell's equation and the linear momentum equation for the dielectrics. As the analysis of electric field distribution in the free space surrounding the solid continuum shows that the special field intensity is rather tiny except that in the region adjacent to the circumferential edge, the vacuum domain is not included in our numerical model. However, its mechanical influence is take into account by applying the electric surface traction on the outer circumference, which equals to the Maxwell stress tensor in vacuum. The numerical calculation shows that the electric-induced deformation will be overestimated if the influence of surrounding medium is not taken into account and the discrepancy is considerable for the large deformation case. Despite that only static case is considered in this work, it can be easily extended to the simulation of time dependent behavior by including the inertia terms in the governing equations. As the commercial package COMSOL is user friendly and relatively easy to master, the proposed approach may help improve the design of more complex dielectric elastomer actuators.