A detailed solution framework for the out‐of‐plane displacement of circular dielectric elastomer actuators

Abstract

The out‐of‐plane displacement (OPD) of a circular dielectric elastomer actuator (DEA) membrane has been explored in recent years for a variety of important applications. Circular DEAs consist of an elastomer membrane that is prestretched to a rigid frame and coated with compliant electrodes. Such a simple configuration has found many interesting applications such as in pumps, pulse tracking, hopping search and rescue robots, dielectric elastomer generators for renewable energy harvesting, linear actuators, and many others. Here, we present an easy‐to‐implement solution framework of the OPD for circular DEAs via a detailed numerical solution using the shooting method (SM) on a system of differential algebraic equations (DAE). The material model incorporated follows the robust Gent hyperelastic representation, and a comparison (of the OPD of a single‐cone DEA) against a neo‐Hookean based solution is included. The SM‐based numerical solution to the resulting boundary value DAE problem is applied to a practical range of values based on the reported literature. Analysis of the results and comparison against other studies are provided. The current work provides a go‐to framework for implementation in further research and development that can be useful to a broader audience including nonexperts in the field of soft robotics.

INTRODUCTION

A circular dielectric elastomer actuator (DEA) can consist of an elastomeric membrane that is prestretched over a rigid frame and is coated with compliant electrodes. Depending on the desired actuation, numerous configurations of electrodes can be used and applied on opposite areas of the membrane. The membrane is then displaced out of plane using a variety of means such as a spring, biasing mass, or magnet.¹ This “cone” is then utilized as a means of inhomogeneous actuation, sometimes in conjunction with other mechanisms (such as attached to a multi-linkage system, as in Reference 2). Actuators made in such a way can be used for multiple applications as they can exhibit motion with large ranges, high prescribability, smoothness, reasonable bandwidth of $\mathrm{𝒪}(10^0)$ to $\mathrm{𝒪}(10^3)$ Hz, and capability of resonance. Due to their broad potential, many applications are pursued by nonexperts in the field of soft robotics (particularly DEAs). However, no comprehensive and robust solution framework for the combined analytical and numerical solution has been developed for the out-of-plane displacement (OPD) of circular DEAs. Such a framework can be used as a go-to resource by a broader group of both developers and researchers.

Recognizing its various applications, He et al.³ pioneered a derivation for the equations of state for a DEA circular membrane undergoing axisymmetric OPD. The obtained system of differential algebraic equations was based on the neo-Hookean (NH) hyperelastic model and constant-permittivity assumptions. While constant permittivity is a reasonable assumption for a wide range of stretches, it is known that the NH model breaks down even at relatively low stretches. Furthermore, for materials that possess low limiting stretch invariants, the NH model can inaccurately represent the material's mechanical behavior. More robust models such as the Arruba–Boyce (AB)⁴ and the Gent model (GM)⁵ can predict not only strain softening behavior at low stretches but also strain stiffening behavior at large stretches, which are typical behaviors of many elastomers. In a comparative study, Boyce⁶ compared their AB model against the GM and found that they led to almost identical results. However, the GM model is simpler to implement than the AB modeling scheme. Therefore, herein a detailed solution framework for the OPD of a circular DEA is developed using the robust Gent hyperelastic model, which allows for more accurate modeling at higher stretches. In addition, the effects of varying force, voltage, and geometry are explored; and a comparison between GM-based versus NH-based OPD is provided. The resulting work provides an easy-to-follow framework for implementation of the governing equations that can be useful to a broader audience of both researchers and developers in various fields. Examples of how this solution framework can be implemented are presented using practical values of the parameters involved.

To emphasize the relevance of the proposed work, in the sequel some applications of devices that rely on the OPD of circular DEAs are briefly surveyed. Cao et al.⁷ proposed a novel pneumatic diaphragm pump where two circular membranes are held apart via permanent magnets and deformed out of plane by applying a voltage across the compliant electrodes. Linnebach et al.⁸ proposed a different variation of a pneumatic diaphragm pump. Their pump consists of a circular DEA that is deformed out of plane using a spring. The DEA is used in conjunction with a diaphragm that is controlled by actuating the DEA. Ye et al.⁹ used a circular DEA deforming out of plane to validate the pulse-tracking ability of a conical DEA. Using an $H_{\infty }$ control, the device was able to track the pulse signal albeit with better results for the low-frequency components of the signal. Dubowsky et al.¹⁰ investigated the use of a circular DEA deforming out of plane to power a hopping search and rescue robot. This 46 g robot was comprised of a spherical housing with a power source providing 8.8 kV and a single DEA that would be actuating 35 times to store energy in a spring via a ratcheting mechanism prior to launch. Bortot and Gei¹¹ explored the feasibility of using a circular DEA deforming out of plane as a dielectric elastomer generator to harvest renewable energy. This was accomplished through a four-phase cycle of mechanical loading, electrical charging, mechanical discharging, and electrical discharging. An efficiency higher than 26% was achieved at a max load of 80 N. Berselli et al.¹² used a singular, circular DEA that was deforming out of plane to generate linear motion. A compliant frame was used to create either unidirectional or bidirectional devices that could reset themselves to the starting position. Cao and Conn¹ optimized the stroke and work output of a circular DEA that was deformed out of plane through multiple means: spring, biasing mass, and rigid support in a double cone configuration.

The current manuscript is organized in the following manner: first, the OPD of a circular DEA is rederived in detail following the work of Reference 3. Then the hyperelastic Gent material model is incorporated and the corresponding state equations are derived. Using a dimensional analysis, dimensionless quantities are obtained and incorporated into a system of algebraic differential equations. Next, a numerical method that relies on the shooting method is implemented on the resulting boundary value problem. Finally, the results and analysis are provided to exemplify the use of the presented framework solution.

SYSTEM MODELING Out-of-plane deformation

The deformed shape of a circular membrane is an axisymmetric cone if the out-of-plane deformation occurs at the center of its circular boundary. When this deformation is applied to dielectric elastomeric membranes, those membranes are referred to as a single cone dielectric elastomer actuator(SCDEA)*

Obviously, this is true for any deformable membrane but our focus is on DE membranes. In addition, an extension of SCDEA, the double conical dielectric elastomer actuator (DCDEA)¹³ involves other complexities, due to additional degrees of freedom,¹⁴ not considered here. For DCDEA applications see References 15-17.

as is shown in Figure 1. To avoid very large stresses, a distributed loading is applied over some finite area in the neighborhood of the membrane's center.¹⁸ In practice, a SCDEA is created as follows. A circular disk is cut out of an elastomeric membrane of thickness T. This membrane is then prestretched over a rigid ring of radius b after which it is coated on both sides with a compliant electrode. Subsequently, the membrane is deformed out of plane by means of a force F applied to the rigid disk of radius a. There are multiple ways of applying this force seen in Figure 1C.^1,19 In this case, the force is assumed to be constant. The origin of the system O is taken at the intersection of the axis of symmetry and the top face of the rigid ring $\mathrm{ð’\circledR }$. Points Aand B on the membrane are located at the boundaries of the rigid disk and rigid ring, respectively, as seen in a cross-sectional view (see Figure 1C). An arbitrary point of interest is denoted as R. Point R starts on the undeformed membrane and moves to a new position on the deformed membrane that is denoted by the coordinates $r(R)$ and $z(R)$.

FIGURE 1. Diagram of a single cone dielectric elastomer actuator (SCDEA). (A) Side, perspective, and top views of unstretched and undeformed dielectric elastomer membrane. (B) Perspective view of prestretched and deformed DE membrane due to the action of a distributed loading acting downward (not shown in the figure) via a rigid center disk, forming a conical shape. (C) Cross-section of stretched and deformed DE membrane showing a concentrated force, F, and direction of stretches. (D) SCDEA showing a unit volume ring with infinitesimal thickness, δz, which is emphasized in (E). (F) Stretch state of unit volume in the ring material showing main stretches, λ1, λ2, and λ3. Compare these stretches with those shown at point (r(R),z(R)) in subfigure (C)

In order to properly examine the kinematics of point R, the stretches present in the membrane must first be defined. This elastomer is modeled using three principal stretches that characterize the membrane in the deformed state: ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$. The longitudinal stretch ${\lambda }_1$ describes the deformation in the direction from the rigid ring to the center disk. The latitudinal or circumferential stretch ${\lambda }_2$ points in the circumferential direction around the cone and at a right angle to ${\lambda }_1$. Lastly, ${\lambda }_3$ is in the thickness direction, and orthogonal to both ${\lambda }_1$ and ${\lambda }_2$. The stretches are illustrated in Figure 1F.

Consider a point R on the undeformed membrane and a point $R+dR$ that is infinitesimally adjacent to it. As the membrane is stretched and deformed out of plane the two points move away from each other so that in the final state they are separated by a distance $dl$. Assuming that this infinitesimal distance is a straight line yields the right triangle shown in the cutout in Figure 1C, where $\theta$ is defined as the angle between the horizontal axis and the tangent to the curve at point R in the deformed state. This causes the two legs of the triangle to be $dz=-\sin \theta dl$ and $dr=\cos \theta dl$. Therefore, $dl$ can be written in terms of $dr$ and $dz$: $dl=\sqrt{{(dr)}^2+{(dz)}^2}$. Dividing through by $dR$ (the distance between the points R and $R+dR$ in the undeformed state) provides the longitudinal stretch ${\lambda }_1$: [Image Omitted. See PDF]

The latitudinal stretch, ${\lambda }_2$, is defined as the ratio of the circumference at the point R in the deformed state versus the undeformed state due to rotational symmetry of the cone about the z-axis. This results in ${\lambda }_2$ being defined as $2\pi r(R)/2\pi R$. By canceling the $2\pi$ out of the numerator and denominator, we arrive at [Image Omitted. See PDF]

The third stretch, ${\lambda }_3$, is a ratio of the thickness of the elastomer in the deformed state versus the undeformed state. Using the assumption of incompressibility, the following equation holds true:²⁰ $1={\lambda }_1{\lambda }_2{\lambda }_3$. Therefore, ${\lambda }_3$ can be written in terms of ${\lambda }_1$ and ${\lambda }_2$ as [Image Omitted. See PDF]

The compliant electrodes on each side of the elastomer act like soft plates in a compliant parallel-plate capacitor, and the membrane is the medium. In this idealized model, it is assumed that the electrodes have negligible stiffness. As a voltage, V, is applied across these electrodes, a charge accumulates. Here, Q is the total electric charge on the conical elastomer. As the charge builds up, the two opposite electrodes attract each other. However, since the material is considered incompressible, ${\lambda }_1$ and ${\lambda }_2$ increase as ${\lambda }_3$ decreases, resulting in the motion of the DEA.

The charge on the electrodes is found using the nominal electric displacement $\tilde{D}$ which is the charge on a deformed element divided by the area of an undeformed element. This definition lends itself to easy calculation of the total charge since all calculations can be done with respect to the undeformed elastomer which has a much simpler geometry. Since $\tilde{D}$ varies at each R position, the total charge on the deformed elastomer has to be calculated by integrating across the surface of the undeformed elastomer using infinitesimal rings of width $dR$. The lower and upper bounds of the integral would be A and B respectively since they bound the active area. [Image Omitted. See PDF]

This system is then analyzed using the Helmholtz free energy, $\mathscr{H}$, under the assumption that the DEA is an isothermal and isochoric system.²⁰ By letting W be the Helmholtz free energy density of an element in the deformed state with respect to the volume of the undeformed elastomer, one can calculate the total Helmholtz free energy of the system: [Image Omitted. See PDF] where $\phi$ is the azimuthal angle that sweeps around the z-axis. In our case this equation could be simplified to a single integral since there is symmetry about the z-axis.

Using the Helmholtz free energy approach, the energy can be analyzed via two pairs of work-conjugates: the force and displacement and the voltage and charge. This approach combines the work done through mechanical and electrical means to achieve equilibrium. A small change in the Helmholtz free energy can be represented by a small variations of the vertical position or charge: [Image Omitted. See PDF]

By examining a small delta in W, we get the following equation due to the mechanical work done in the ${\lambda }_1$ and ${\lambda }_2$ directions and the electrical work: [Image Omitted. See PDF] where $s_1$ and $s_2$ are the nominal stresses in the ${\lambda }_1$ and ${\lambda }_2$ directions, respectively. Also, $\tilde{E}$ is the nominal electric field defined as $V/T$. Suo's work on dielectric elastomers²⁰ provides a detailed derivation of these fundamental equations.

In order to further understand Equation (7), consider a small variation in ${\lambda }_1$ that is caused by small changes in r and z. The small change in ${\lambda }_1$ can be obtained by manipulating Equation (1). Using the cutaway in Figure 1C and the Pythagorean theorem, the equation $d{l}^2=d{r}^2+d{z}^2$ follows naturally. By writing one of the $dr$ terms and one of the $dz$ terms in terms of $\theta$ and $dl$, $d{l}^2=dl\cos \theta dr-dl\sin \theta dz$. Dividing through by $dl$ reduces the equation to $dl=\cos \theta dr-\sin \theta dz$. Dividing through by $dR$ results in an alternative equation for ${\lambda }_1$: ${\lambda }_1=\cos \theta \frac{dr}{dR}-\sin \theta \frac{dz}{dR}$. A small variation in ${\lambda }_1$ is achieved by varying r and z by small amounts. This results in ${\lambda }_1$ becoming $\delta {\lambda }_1$, r becoming $\delta r$, and zbecoming $\delta z$ which leads to Equation (8): [Image Omitted. See PDF]

By the same reasoning, Equation (2) can be rewritten as [Image Omitted. See PDF]

Substituting Equation (7) into Equation (5) and setting the result equal to Equation (6) leads to [Image Omitted. See PDF]

By substituting Equations (8) and (9) into Equation (7) and the result into Equation (10), the following expression emerges for the left-hand side of the aforementioned equation: [Image Omitted. See PDF]

Splitting the integral into individual terms leads to [Image Omitted. See PDF]

The last two terms of this expression are already completely reduced. However, the first two terms can be further simplified. Integration by parts will be used for this ($\int udv=uv-\int vdu$). For the first term, choose $u=R{s}_1\cos \theta$ and $dv=\frac{d(\delta r)}{dR}dR$. Integrating by parts results in $(R{s}_1\cos \theta \delta r){|}_A^B-{\int }_A^B\delta r\frac{d(R{s}_1\cos \theta }{dR}dR$. The second term in solved in a similar manner with $u=R{s}_1\sin \theta$ and $dv=\frac{d(\delta z)}{dR}dR$ resulting in $(-R{s}_1\sin \theta \delta z){|}_A^B+{\int }_A^B\delta z\frac{d(R{s}_1\sin \theta )}{dR}dR$. Grouping like terms together from these and the other two terms in Equation (12) results in Equation (13). [Image Omitted. See PDF]

Since $\delta r(R)$, $\delta z(R)$, and $\delta \tilde{D}(R)$ vary independently of each other, their respective terms can be isolated and examined independently of the others. This is accomplished by setting the two variables that are not being considered to zero.

Set $\delta r(R)=0$ and $\delta \tilde{D}(R)=0$ to look at the $\delta z(R)$ terms while using Equation (12) as the left-hand side for Equation (10).

[IMAGE OMITTED. SEE PDF.]

[IMAGE OMITTED. SEE PDF.]

Eliminating $dR$ from the left-hand size of Equation (15) and noticing that the linear distributed loading $2\pi s_{1}R\sin (\theta )$ is constant†

This is based on the assumptions that the material is homogeneous and that after a force, F, is applied, equilibrium is allowed to be reached.

throughout the membrane (i.e., from A to B), then, in particular, it must independent of $\delta z$. Therefore, one can rewrite Equation (15) as: [Image Omitted. See PDF]

The small vertical displacement $\delta u$ of the center disk must be equal to the total contribution of each small $\delta z$ at every unit volume of the material (see Figure 1), that is, [Image Omitted. See PDF] where the negative sign is due to the convention originally established (see Figure 1). Therefore: [Image Omitted. See PDF]

In order to look at $\delta r(R)$ terms independently, set $\delta z(R)=0$ and $\delta \tilde{D}(R)=0$.

[IMAGE OMITTED. SEE PDF.]

[IMAGE OMITTED. SEE PDF.]

Let us consider Equation (20) briefly. The term on the right-hand side of Equation (20) must vanish since all small contributions of $\delta z$ have been set to zero; therefore, there is no external work in the vertical direction. Now, let us examine the first term in Equation (20). For every initial particle of material R, the component principal internal force $2\pi RT{s}_1\cos \theta$ acts perpendicular to the direction of the external force F in Figure 1. Therefore the small internal work produced by such a horizontal force due to small horizontal displacement $\delta r$ is given by $2\pi RT{s}_1\cos \theta \delta r$. By problem definition, the work done at the fixed boundary particles A or B must vanish and thus their difference $2\pi RT{s}_1\cos \theta \delta r{|}_A^B$. Therefore: [Image Omitted. See PDF] or, since $2\pi T$ is nonzero, [Image Omitted. See PDF] which corresponds to the first term in Equation (20). In turn, it can be rewritten as: [Image Omitted. See PDF] which leads to a function for the nominal stress in the ${\lambda }_2$-direction as: [Image Omitted. See PDF]

Finally, by setting the two mechanical variables to zero, $\delta r(R)=0$ and $\delta z(R)=0$, the electrical equation can be acquired. For this case the equation reduces to

[IMAGE OMITTED. SEE PDF.]

We can then turn the charge equation, Equation (4), into a small charge, $\delta Q$, by changing $\tilde{D}$ to a small change, $\delta \tilde{D}$. The area for the integration remains the same; and by substituting the resulting $\delta Q$ expression into Equation (25), we get:

[IMAGE OMITTED. SEE PDF.]

However, the nominal electric field, $\tilde{E}$, is not dependent on R since it is by definition only dependent on the voltage (which is invariant to deformation) and the original thickness, T. This allows $\tilde{E}$ to be pulled out of the left-hand side of the integral and the remaining integrands to be canceled. The resulting equation verifies the nominal electric field, $V=T\tilde{E}$.

Material model

Hyperelastic material models are used to describe the mechanical behavior of the elastomer, providing the equations for $s_1$ and $s_2$. Numerous models are used for DEAs, and descriptions of some of the most common models are described elsewhere.^21-24 Models such as the Gent and Ogden model are widely used. Generally, these models are based on strain energy density. However, since the DEA acts as an electromechanical system, the complete representation must include a model for the dielectric energy density. The material model will be denoted a $W_M$. The dielectric energy density is $W_{DE}$. Combining these two results in: [Image Omitted. See PDF]

In order to simplify the calculations, certain assumptions are made that result in an ideal dielectric elastomer. According to Reference 20, the assumptions for ideal dielectric elastomer are as follows: (i) constant permittivity ‡

It is known, however, that permittivity is dependent on stretch,^25,26 but this assumption is reasonable for a wide range of stretches.

($E=D/\epsilon$) and (ii) incompressibility (${\lambda }_1{\lambda }_2{\lambda }_3=1$). This allows the dielectric energy density to be expressed as: [Image Omitted. See PDF] using the relationship $E=D/\epsilon =\tilde{D}/(\epsilon {\lambda }_1{\lambda }_2)$ where E is the true electric field, D is the true electric displacement, and $\epsilon$ is the absolute permittivity of the material.

The internal stress state and electric displacement are crucial in analyzing the material behavior. The in-plane nominal stresses, $s_1$ and $s_2$, and nominal electric field, $\tilde{E}$, are defined as follows:²⁰ [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

One of the most common material models used in modeling DEAs is the Gent hyperelastic material model.^27-30 The Gent model more accurately describes the behavior of hyperelastic materials for a wider range of stretches than the NH model because of the incorporation of a limiting stretch that models well strain stiffening behavior. This contrasts the current work with that reported in Reference 3, since they used the NH model. The material model part of Equation (27) is displayed below: [Image Omitted. See PDF]

Here $\mu$ is the small strain shear modulus, and $J_m$ is the limiting stretch. Combining this with the dielectric energy density results in: [Image Omitted. See PDF]

Using Equations (29a)–(29c) the following three equations can be derived: [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

NUMERICAL SOLUTION System of equations

The following system of four equations completely describes the motion of the system due to electrical and mechanical inputs. Equations (35a) and (35b) naturally precipitate from Equation (1). Equation (35c) comes from solving Equation (18) for the product $R{s}_1$ and plugging the result into (24). Then the differentiation is carried out, and the result is substituted back into (18). Equation (35d) is just Equation (18) rewritten, and added here for completeness. [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF] [Image Omitted. See PDF]

Besides the aforementioned, several other relationships in this article were used to simplify the system (35) into terms of only four variables, namely, r, $\theta$, z, and ${\lambda }_1$. The final algebraic constraint equation becomes: [Image Omitted. See PDF]

The following expressions were used to nondimensionalize the system of equations for easier comparison: $u/a$, $Q/(2\pi a^2\sqrt{\epsilon \mu })$, $F_n=F/(2\pi aT\mu )$, and $V_n=V/(T\sqrt{\mu /\epsilon })$, as suggested in Reference 3. Applying the nondimensionalization and $V/T=\tilde{D}/\epsilon {\lambda }_1^2{\lambda }_2^2$ from Equation (34) and $V=T\tilde{E}$ results in the following algebraic constraint equation: [Image Omitted. See PDF]

Additionally, the system is subjected to the following boundary conditions (BCs): [Image Omitted. See PDF] which result from geometrical constraints. The first two BCs are defined by the radius of the inner and outer rigid supports, respectively. The third BC is due to the definition of the origin of the r–z coordinate system.

The resulting system of equations is what is known as a differential-algebraic equation (DAE) system. As the name implies these DAE systems contain both differential and algebraic equations, where the algebraic equations are often constraint equations. We will focus on semiexplicit DAE systems such as ours that take the form:³¹ [Image Omitted. See PDF] [Image Omitted. See PDF]

Here the $\mathbf{\text{x}}$ variables are those that occur both in derivative form and algebraic form. Whereas the $\mathbf{\text{y}}$ variables only occur in the algebraic form even though they may appear in both the $\mathbf{\text{f}}$ and $\mathbf{\text{g}}$ function, and t is the independent variable. This is what separates an ODE system from a DAE system.

Implementation of numerical method

A boundary value problem (BVP) DAE system as given by Equations (35)–(38) can be solved numerically. Here, a more comprehensive (as compared with References 3 and 1) numerical algorithm is provided for solving this system, with the goal of providing a framework for easy implementation. Information on how to solve BVP DAE systems is scarcely found. Ascher and Spiteri³² created a DAE solver (COLDAE) in FORTRAN capable of solving semiexplicit BVP DAEs up to an index of 2. However, as cautioned by the authors in Reference 32, their algorithm is not successful for all applications. Since the authors did not find a suitable package that can directly handle these BVP DAE systems, a customized solution was devised. One way of solving a BVP is to transform it into an initial value problem (IVP). This is accomplished through the shooting method, which is summarized as follows (see flowchart shown in Figure 2). An initial value is guessed; the solution is calculated and the result is compared with the known BC; then an adjustment is made to the initial guess; this updated guess is used to repeat the process until the solution converges. An existing ODE solver was used inside of our implementation of the shooting method algorithm, which is shown in Figure 2. The algebraic equation was used to generate consistent initial conditions.

FIGURE 2. Flowchart of algorithm for solving the DAE system based on a shooting method. Note that an extra step is needed in order to ensure that consistent initial conditions are computed based on the nondimensionalized algebraic obtained from Equation (18), and given by Equation (37)

The developed algorithm follows the logic of the flow chart in Figure 2. $MATLA{B}^{\circledR }$ was chosen to solve this system of equations due to its perceived familiarity and availability. The algorithm begins with the definition of the global parameters and an initial guess for ${\lambda }_1$. The input parameters defined were $F_n$, $V_n$, a, b, and ${\lambda }_p$. The ${\lambda }_1$ guess is based on the prestretch for the particular application. Next, Equation (37) is solved for $\theta$ to directly calculate a value for $\theta$ based on initial conditions: r, R, and the guessed ${\lambda }_1$. However, the ${\lambda }_1$ value is kept as a variable, and only the r and R values are defined using the BCs at B. Now that a set of consistent initial conditions is defined, the system of equations, Equations (35a)–(35c) and (37) can be written in matrix form according to Equations (41) and (42), where M is the mass matrix: [Image Omitted. See PDF]

Where [Image Omitted. See PDF]

This system is then solved over the span $[B,A]$ using the two conditions at B defined above, the calculated $\theta$ value, and the variable for the unknown ${\lambda }_1$. The first-order DAE system is then solved using the ode15s function in ${\text{MATLAB}}^{\circledR }$. Then, the remaining BC is applied ($r(A)=a$). This means that the calculated value of $r(A)-a$ should equal zero if the BC is satisfied. The fsolve function is used to drive this quantity to zero by adjusting the value of ${\lambda }_1$. This iteration is the shooting method portion of the code. Once convergence is achieved, the ${\lambda }_1$ guess is updated with the current ${\lambda }_1$ value. This will be used as the next guess for the next voltage term. This becomes important as the voltage terms become larger and the difference between ${\lambda }_1$ values becomes larger. Increasingly smaller steps might need to be taken to facilitate convergence. Once this occurs, the program ends the loop and reports the following results: R, r, z, $\theta$, ${\lambda }_1$, and ${\lambda }_2$.

RESULTS AND DISCUSSION: IMPLEMENTATION USING MATERIAL MODEL Comparison of OPD of SCDEA using neo-Hookean versus Gent model

The work reported by He et al.³ relied on the NH hyperelastic model. It is well known that although the NH can predict strain softening behavior at low stretches, it does not model strain stiffening behavior (at larger stretches), which is typical in many elastomeric materials. Therefore, it is expected that the stretch and stress states of the cone configuration at locations of relatively high stretches would be misrepresented when using the NH model. This state of affairs can be visualized from the results shown in Figures 3 and 4, which were produced for identical geometries and for a material having a $J_m$ of 100 under a 120% prestretch. Both longitudinal stretch and stress at the location closest to the inner rigid disk were plotted, versus dimensionless force $F_n$ in Figure 3, and versus dimensionless voltage $V_n$ in Figure 4. As it might be noted, for this particular material:

• The stretch–force behavior is linear if NH model is used; however, the more realistic GM exhibits a nonlinear stretch–force relationship.
• The larger the force, the larger the discrepancies between the NH stretch and the GM stretch.
• The longitudinal stress predicted by the NH model is overestimated for $F_n\lesssim 2$. Furthermore, this overestimation becomes more pronounced as $F_n$ increases. In some practical cases, the value of $F_n$ can reach a value of about 5, in which case the NH model would predict a ${\sigma }_1/\mu$ ratio 35% higher than the more realistic value predicted by the GM model.
• Although both the stretch–voltage and stress–voltage behaviors can exhibit similar trends for both NH and GM, yet the NH values are consistently larger than the GM values. It might be worth mentioning that for larger $F_n$ values, the stretch–voltage and stress–voltage behaviors can be different (see Figure 5).

FIGURE 3. Comparison of the neo-Hookean and Gent models for a range of maximum λ1 stretches as a function of force, Vn=0.2, Jm=100, b/a=4, λp=1.2

FIGURE 4. Comparison of the neo-Hookean and Gent models for a range of maximum λ1 stretches as a function of voltage, Fn=1, Jm=100, b/a=4, λp=1.2

FIGURE 5. Comparison of the neo-Hookean and Gent models for a range of maximum λ1 stretches as a function of voltage, Fn=3, Jm=100, b/a=4, λp=1.2

In addition, it is a known that GM approaches the NH model as the limiting stretch invariant, $J_m$, goes to infinity. Therefore the smaller the value of $J_m$ the less accurate the NH model becomes and the larger the differences between values predicted with both models. For example, values of $J_m$ for commonly available DE materials can varied from about 10²⁸ to about 200.³³ The impact of such variation of values can be observed in the plots shown in Figure 6, which were developed using the GM hyperelastic model. Note the remarkable difference between, for example, the plots shown in Figure 6A,C. An identical plot to that shown in Figure 6C (for $J_m=1000$) would be obtained if the NH model were to be used.

FIGURE 6. Effect of increasing limiting stretch on λ2 for varying voltage, Fn=2, b/a=4, λp=1.2. Note that the dimensionless voltage, Vn, is varied from 0 to 0.3

General behavior

In order to gain more understanding into the parameter-dependent behavior of the SCDEA, the force, voltage, and geometry were varied. In addition, a few example cases are evaluated using a typical range of reported values in the literature, which provides the reader with a sense of practical OPDs, stretches, and stresses. Figure 7 is a cutout from Figure 1 and is included here for convenience as a reference to visualize the stretches. Unless otherwise specified, the following parameters were used throughout this section: a rigid disk radius of $a=0.008$ m, a rigid ring radius of $b=4a$, ratios relating the prestretch $A=a/1.2$ and $B=b/1.2$, and $J_m=100$. However, in the following subsection, some other examples were considered that used data directly from literature that represents a wider use case. It might be worth mentioning that a prestretch of 1.2 is on the lower bound and more typically used for silicone rubbers. For some acrylic-based DEA, such as VHB 4910/05, higher values are typically used. (For example, in other experimental works,^16,24 the authors have typically used a prestretch value closer to 3 for VHB 4910.)

FIGURE 7. Clarifying cutout for stretches

The first scenario explored was that of a constant voltage with varying applied forces. Figures 8–11 show how $F_n$ ranging from 1 to 3.5 affects the shape of the cone as well as the stretches within the cone at certain fixed voltages.

FIGURE 8. Dimensionless deformation and stretch graphs for various applied dimensionless forces at Vn=0. Note that the dimensionless force, Fn, is varied from 1 to 3.5

FIGURE 9. Dimensionless deformation and stretch graphs for various applied dimensionless forces at Vn=0.1. Note that the dimensionless force, Fn, is varied from 1 to 3.5

FIGURE 10. Dimensionless deformation and stretch graphs for various applied dimensionless forces at Vn=0.2. Note that the dimensionless force, Fn, is varied from 1 to 3.5

FIGURE 11. Dimensionless deformation and stretch graphs for various applied dimensionless forces at Vn=0.25. Note that the dimensionless force, Fn, is varied from 1 to 3.5

The purely mechanical case (i.e., when the voltage is zero) is first explored and results are shown in Figure 8. Here the general behavior can be observed. As would be expected and can be seen in Figure 8A, the cone elongates in the z-direction under increasing load from approximately $-$1.6 to $-$4.7. Figure 8B shows how the longitudinal stretch increases as more force is applied. Note that there is a much more drastic change in longitudinal stretch closer to the rigid disk. This is due to the fact that the tangent angle farther down the cone (closer to the tip) is larger; therefore, the force is directed more parallel to the ${\lambda }_1$-direction near the tip of the cone while less near the base. Figure 8C§

For direct comparison, the axes were maintained the same as in Figures 9–11.

shows the circumferential stretch (${\lambda }_2$) in the cone. This stretch is restricting the cone from expanding in diameter. The range of ${\lambda }_2$ is smaller than that of ${\lambda }_1$. Also, notice that for $V_n=0$, ${\lambda }_2$ does not reach the initial prestretch as the force exerted on the cone is increased, except at the boundary values. Since ${\lambda }_2$ is less than the initial prestretch, this means that there is a small degree of relaxation in the circumferential direction for all applied forces. However, the stretch remains above unity, indicating that there is still some restrictive force on the cone; otherwise, electromechanical instability could occur.

Similar results were produced for $V_n=0.1$, $V_n=0.2$, and $V_n=0.25$, and shown in Figures 9,10, and 11, respectively. Note that in those figures, the overall trend for subfigures (A) is an increase in the $z/a$ deformation, as should be expected. The $z/a$ deformation at the highest force was around 4.7 for the zero voltage case compared with about 5.1 for $V_n=0.25$. This difference, while relatively small, is due to the actuation and thinning of the elastomer. As the elastomer is compressed in thickness the other two stretches must compensate by expanding the overall surface area. This extra area is invested toward the deformation in the height of the cone. Subfigures (B), in general, maintain a typical shape with only minor change in amplitude more on the $r/a=1$ side (note that as force increases, these changes become more noticeable). Overall, the subfigures (C) undergo the most noticeable change. Even with slight actuation, a large spike in the circumferential stretch can be observed. This is caused when the cone is deforming in the z-direction and tries to stretch the cone into a more straight-line case where the radius, r, would increase. The circumferential stretch is the counteracting stretch in this case. However, points that were once at a tighter radius inevitably get repositioned to larger radius. A large difference can be seen here due to actuation.

The second scenario explored was that of a constant force with varying applied voltages. Figure 12 shows how dimensionless voltages ranging from 0 to 0.25 affect the circumferential stretch within the cone at different discrete dimensionless force values ranging between 1 and 4. Only ${\lambda }_2$ is shown here because the various $z/a$ and ${\lambda }_1$ graphs did not exhibit different behaviors compared with those shown in Figures 8–11. Low force causes a low ${\lambda }_2$ value regardless of the applied voltage. As the force is increased, the amount of voltage inputted into the system greatly affects the amount of stretch possible in this direction, and not in a linear manner. The highest ${\lambda }_2$ value achieved was with the highest force ($F_n=4$) and highest voltage ($V_n=0.25)$, resulting in a peak ${\lambda }_2$ value of approximately 1.28. This is an increase of 7% over the initial prestretch.

FIGURE 12. Circumferential stretch graphs for various applied dimensionless voltages at various fixed forces. Note that the dimensionless voltage, Vn, is varied from 0 to 0.25

Next, the effect of changing the geometry of the cone was examined. Figures 13,14, and 15 relate to ring-to-disk ratios of 1.5, 4, and 6, respectively. (Here the vertical axes in all subfigures (B) and (C) were kept the same, whereas those for subfigures (A) were varied due to a large range of deformations.) A dimensionless force of 2 was again used for these plots. Considering Figure 13, as the disk radius approaches that of the ring, the edges of the cone become straighter and more in-line with the applied force. This affects the longitudinal stretch, resulting in a ${\lambda }_1$ graph that appears a lot more linear than previous cases. As the tip of the cone becomes narrower (i.e., increase of $b/a$ ratio), the radius of curvature decreases which results in more extreme changes in ${\lambda }_1$ (see Figure 15). As can be seen, Figure 15B looks much more sharply decaying than the nearly linear case seen in Figure 13B even though they all have similar values near $r/a=1$. The ${\lambda }_2$ graphs become less symmetrical with larger spikes both above and below the prestretch level. There is a stretch reduction for lower actuation levels and a stretch peak for higher actuations. Both the maximum and minimum values here occur well toward the rigid disk side. It is interesting to note that the overall deformation in the z-direction greatly increases as the ratio of $b/a$ increases. A ratio of 1.5 has a deformation of $z/a$ less than 1 while a ratio of 6 results in a deformation greater than 4. This is all with the same force. This is due to the fact that the smaller disk is able to exert more pressure than a larger disk with the same applied force. The actual actuation then only causes small changes in the z-direction.

FIGURE 13. Dimensionless deformation and stretch graphs for ring to disk ratio of 1.5. Note that the dimensionless voltage, Vn, is varied from 0 to 0.3

FIGURE 14. Dimensionless deformation and stretch graphs for ring to disk ratio of 4. Note that the dimensionless voltage, Vn, is varied from 0 to 0.3

FIGURE 15. Dimensionless deformation and stretch graphs for ring to disk ratio of 6. Note that the dimensionless voltage, Vn, is varied from 0 to 0.3

Examples

In the sequel, three examples are presented from literature that utilize a variety of elastomeric materials, prestretches, and geometries. These examples show the ability of the presented framework to model the common range of operating parameters found in literature.

Example 1: General material case

The first example is based on data from He et al.³ The three graphs that define the deformation and two principal in-plane stretches can be seen in Figure 16. The following parameters were used to generate the graphs: a rigid disk radius of $a=0.008$ m (this was not given but assumed and is somewhat arbitrary as the nondimensionalized axes cancel out any dependence on a), a rigid ring radius of $b=4a$, and ratios relating the prestretch $A=a/1.2$ and $B=b/1.2$. The nondimensionalized force $F_n$ was set to 2 and the dimensionless voltage values were varied to be $V_n=\{0,0.15,0.2,0.25,0.3\}$. Also, a limiting stretch of 100 was assumed as no material was specified for this theoretical study.

FIGURE 16. Example of possible state based on data from Reference 3, using the following parameters: material not specified, b/a=4, λp=1.2, Fn=2, Jm=100, Vn={0,0.15,0.2,0.25,0.3}

Figure 16A shows the vertical cross section of the conical DEA and how it behaves during actuation. The inner-most profile is the one where no voltage is applied. As voltage is applied the cone elongates and bows outward to become straighter. The elongation that occurs during actuation is due to the increase in ${\lambda }_1$ as can be seen in Figure 16B, and the straightening of the edges is due to the increase of ${\lambda }_2$ at points away from the boundary. Interestingly, the peak ${\lambda }_2$ value shifts from the boundary at zero voltage to about $r/a=1.5$ as the voltage is increased to $V_n=0.25$ despite ${\lambda }_1$ always being strictly decreasing as one moves from the rigid disk to the outer ring ($r/a=1\to r/a=4$). This aptly demonstrates the nonhomogeneity of the stretches within the cone. However, at these forces and voltages, ${\lambda }_1$ clearly dominates as the ${\lambda }_2$ values vary only slightly in comparison. These graphs follow the trend of the results from Reference 3. However, there is some variation, especially in the regions of larger stretches, as the Gent model is being used.

By combining data from Figure 16B,C, Figure 17 can be generated as ${\lambda }_3$ is just $1/{\lambda }_1{\lambda }_2$. As would be expected, the general shape of this graph is a reflection of the ${\lambda }_1$ graph about the horizontal axis as that behavior dominates. The thickness of the membrane is reduced by about half as one moves from the ring to the disk. Experimental work that has been performed in our lab has shown that failure due to dielectric breakdown most often occurs right around the rigid disk ($r/a=1$) where the thickness is at a minimum.

FIGURE 17. Thickness stretch. Note that the dimensionless voltage, Vn, is varied from 0 to 0.3

Example 2: VHB 4910 case

The second example considered data from Guo et al. on antagonistic membranes deforming out of plane.²⁷ Here two membranes are deformed out of plane and attached at the center (where the rigid disk is located). This actuator acts like a pump that mimics the heart. The acrylic based elastomer VHB 4910 was used with an initial thickness $T=1$ mm. The ring-to-disk radius ratio was 2.6667 with $b=20$ mm and $a=7.5$ mm. A prestretch of ${\lambda }_p=2.5$ was used in this case, which is notably greater than in the other graphs so far. A dimensionless force of $F_n=1.5$ was chosen as no force was directly reported. The dimensionless voltage values chosen exceed those used by Guo et al.²⁷ and were chosen to be $V_n=\{0,0.1,0.15,0.2,0.25\}$ The value of the small strain shear modulus for VHB 4910 was reported as 40 kPa with the absolute permittivity being $3.98\times 10^{-11}$ F/m. The resulting graphs can be seen in Figure 18. The lower stiffness of the VHB 4910, as evidenced by the higher value of the limiting stretch, allows for greater actuation response.

FIGURE 18. Example of possible state based on data from Guo et al.27 The following parameters were used: VHB 4910, T=1 mm, b/a=2.667, λp=2.5, Fn=1.5, Vn={0,0.1,0.15,0.2,0.25}, Jm=120, μ=40 kPa, ε=3.98×10−11 F/m

Example 3: Silicone rubber material case

Finally, a much thinner and stiffer membrane material is considered for the last example. Cao et al.²⁸ utilized a silicone membrane material in the creation of a magnetically coupled actuator with two membranes deforming out of plane. The material used was a $100\mathrm{\mu }$m Elastosil with a small strain shear modulus of 431.5 kPa and a limiting stretch of 11.35. Also, the absolute permittivity of the material is $2.479\times 10^{-11}$ F/m. The following geometry was used: $b=15$ mm, $a=7.5$ mm, $b/a=2$. A prestretch of 1.2 was used, and a dimensionless force of $F_n=3$ and dimensionless voltages of $V_n=\{0,0.15,0.2,0.25,0.3\}$ were applied. This exceeds the dimensionless voltage applied by Reference 28 of 0.263. However, even with these large forces and voltages, only small deformation is achieved as can be seen in Figure 19. This is due to a combination of factors. One, the material is extremely stiff with a low limiting stretch. This results in the strain stiffening effects to influence the amount of deformation for even small stretches, showing the importance of using a model like the Gent model for these types of materials. Two, the geometry used in this article is not conducive to large deformations. As was observed in Figure 13A, for small $b/a$ ratios the amount of deformation is less than geometries with larger $b/a$ ratios even with reasonably large forces and voltages.

FIGURE 19. Example of possible state based on data from Cao et al.28 The following parameters were used: Elastosil, T=100μm, b/a=2, λp=1.2, Fn=3, Vn={0,0.15,0.2,0.25,0.3}, Jm=11.35, μ=431.5 kPa, ε=2.479×10−11 F/m

Despite having different values for T, $F_n$, and $V_n$ (as compared against the previous examples), the model is still able to predict accurately the behavior of the actuator. However, the modeling framework only seems to break down at values for which the electromechanical material model itself breaks down, as well. As it turns out, at high voltage values, the ratio of $s_2$ to $s_1$ can diverge for higher $\theta$ values (near the rigid disk) and become a double-valued function along ${\lambda }_1$-direction. Also, negative stress values occur for $s_2$ near the rigid disk at higher voltages. Interestingly, this is the location where electromechanical instabilities such as wrinkling has been observed during experimental work carried out by the authors.

It is worth noting that different values of terms can still result in similar behavior in the graphs. That is due to the fact that the dimensionless terms of force and voltage can mask what is going on behind the scenes. For instance, for a given dimensionless force, the radius of the disk could be doubled and the original thickness of the elastomer cut in half with no noticeable effect on $F_n$. Effort has been made to encompass as large of a working space as possible to make the results useful to the largest possible audience. However, it is possible that extreme cases would not be modeled without modification to the code.

CONCLUSIONS

The various applications of DEAs as promising transducers for robots, energy harvesting, medical devices, and so on, make their study relevant and important. The axisymmetric OPD for a circular DEA membrane has been studied. Based on theoretical work previously derived, a detailed numerical solution framework has been developed, expanded to include a more robust hyperelastic model, and presented in this work. A similar approach can be used to incorporate other material models based on specific applications. Implementation of the proposed framework was shown in the form of practical examples. A range of parameters were used based on reported literature and practical applications. A parametric and comparative analysis was also presented.

The resulting framework has some limitations. The nature of the DAEs rendered the BVP problem solution susceptible to hyperparameters. Therefore, care must be taken when implementing the numerical solution presented. Depending on the initial guessed ${\lambda }_1$, integration steps, and other parameters, the ${\lambda }_1-\theta$ relationship can become double-valued, which can lead to nonconverging solutions. Furthermore, constant permittivity assumptions were used in the modeling approach as this has little effect on the overall behavior. However, it has been shown that permittivity can be stretch dependent.^25,26

Future work (in progress) will expand the proposed framework to include a material model that includes stretch-dependent permittivity as well as viscoelastic effects. The latter is very important for applications that involve out-of-plane motion at various speeds. In addition, future work could include analysis of more degrees of freedom such as cone-tip rotation.

NOMENCLATURE

• a
• rigid disk radius
• A
• position of point on undeformed membrane that will end up at the edge of the rigid disk
• b
• rigid ring inner radius
• B
• position of point on undeformed membrane that will end up at the inner edge of the rigid ring
• D, $\tilde{D}$
• true, nominal electric displacement
• $dl$
• the straight-line distance between R and $R+dR$ in deformed membrane
• E, $\tilde{E}$
• true, nominal electric field
• $\mathbf{f}$
• represents a function
• F
• out-of-plane applied force
• $F_n$
• nondimensionalized force
• $\mathbf{g}$
• represents a function
• $\mathscr{H}$
• Helmholtz free energy of the system
• $J_m$
• limiting stretch in Gent hyperelastic model
• M
• mass matrix
• O
• origin of cone
• Q
• total electric charge on membrane
• r
• radial position
• R
• point of interest on undeformed membrane
• $s_1$
• nominal stress in ${\lambda }_1$-direction
• $s_2$
• nominal stress in ${\lambda }_2$-direction
• $\mathrm{ð’\circledR }$
• plane that coincides with the origin of the cone and the rigid disk
• t
• independent variable in a semiexplicit DAE system
• T
• unstretched, undeformed thickness of elastomer
• u
• overall vertical displacement of cone
• V
• voltage applied across the membrane
• $V_n$
• nondimensionalized voltage
• W
• Helmholtz free energy density
• $W_{DE}$
• dielectric energy density
• $W_M$
• material model contribution to the free energy density
• $\mathbf{x}$
• dependent variables that have a derivative present a semiexplicit DAE system
• $\mathbf{y}$
• dependent variables that do not have a derivative present a semiexplicit DAE system
• z
• vertical position
• $\delta (·)$
• represents a small change in the variable, such as for r, z, $\tilde{D}$, and so on
• $\epsilon$
• absolute permittivity of the elastomeric material
• $\theta$
• angle between the horizontal axis and the tangent line at point R on the deformed membrane
• ${\lambda }_1$
• stretch in the longitudinal direction
• ${\lambda }_2$
• stretch in the latitudinal (circumferential) direction
• ${\lambda }_3$
• stretch in the thickness direction
• ${\lambda }_p$
• prestretch of the membrane
• $\mu$
• small strain shear modulus
• ${\sigma }_1$
• true stress in ${\lambda }_1$-direction
• $\phi$
• azimuthal angle (sweeps around z-axis)
• $\mathit{\chi }$
• matrix containing all of the dependent variables of the system

Conflict of Interest

The authors declare no conflict of interest relevant to this article.

AUTHOR CONTRIBUTIONS

Daniel Korn: conceptualization (supporting); investigation (equal); methodology (supporting); software (lead); writing-original draft (equal); writing-review and editing (supporting). Carson Farmer: investigation (supporting); software (supporting); writing-review and editing (supporting). Hector Medina: conceptualization (lead); investigation (equal); methodology (equal); project administration (lead); supervision (lead); writing-original draft (equal); writing-review and editing (lead).