Dynamics of a Magnetically Coupled Dielectric Elastomer Actuator.

A magnetically coupled dielectric elastomer actuator (MCDEA) is an emerging two-degree-of-freedom DEA system that demonstrates rich dynamic behaviors and promising potential applications in robotics, energy harvesting, and smart structures. However, the three-dimensional geometry, nonlinearity, and viscoelasticity in the system lead to complex nonlinear dynamic behavior that is challenging to predict. In this work, we develop a numerical model that can accurately characterize its dynamic responses. This model, along with experimental results, demonstrates the complex dynamic phenomena of this MCDEA system, including superharmonic, primary-harmonic, and subharmonic resonances, bifurcations, and multiple modes of oscillation. Dynamic control strategies including phase tuning and frequency control are proposed in this work, allowing control of the amplitude and the appearance of a speciﬁc resonance and the occurrence of emerging dynamic behaviors, such as beat phenomena. These ﬁndings have numerous additional applications in areas including active vibrational control, energy harvesting, and programmable soft motors for on-demand locomotion.


I. INTRODUCTION
Along with the rapid development of soft robotics, the soft actuation technology of dielectric elastomer actuators (DEAs) has gained increasing attention in the past decade because of their large actuation strain, inherent compliance, and low cost [1]. The basic unit of a DEA consists of a dielectric elastomer membrane sandwiched between compliant electrodes; when the membrane is subjected to an electric field, the generated Maxwell pressure results in an in-plane expansion and transverse compression of the membrane. On the basis of this transduction mechanism, numerous DEA-based actuators have been proposed (e.g., minimum-energy DEAs [2][3][4][5], stacked DEAs [6,7], roll DEAs [8,9], balloon DEAs [10,11], and cone DEAs [12][13][14][15]). Among these, the simple structure and high force or stroke actuation of the cone DEA are particularly attractive [12,13]. A cone DEA is named after its conical geometry, where an elastomer membrane with a rigid circular frame is deformed out of plane by a protrusion force, which is generated by a central biasing element. The output performance of a cone DEA is highly dependent on its biasing mechanisms, such as a bistable mechanism [14][15][16][17], a deadweight [15,[18][19][20][21], a linear compression spring [13], and an antagonistic mechanism [12,[22][23][24][25][26][27][28]. When a DEA membrane pair is coupled in an antagonistic manner, the biasing force shapes a double conical configuration where the two DEA membranes can be actuated independently to achieve antagonistic actuation (known as a "double-cone DEA"). However, the antagonistic mechanism is conventionally achieved by a rigid attachment between the two DEA membranes such as a centrally bonded region [14] or a rigid rod [12], and that restricts the actuation of the two membranes to be in phase.
A magnetically coupled dielectric elastomer actuator (MCDEA) [27,29] [a schematic diagram of which is shown in Fig. 1(a)] differs from a conventional DEA as the compliant magnetic coupling introduces an additional degree of freedom (DOF) to the system [27]. A pneumatic pump was developed on the basis of this MCDEA configuration and silicone dielectric elastomer and exhibited low viscous loss [29]. The pumping performance of this design was maximized by taking advantage of the compliant magnetic coupling so that one membrane could oscillate freely at resonance to greatly improve the stroke of the pump diaphragm, despite the presence of damping from the compressed air. The resonant actuation of MCDEAs can amplify the stroke or power output and increase the energy efficiency, which opens the potential for many applications that exploit this behavior, such as dynamic locomotion [28,30], active vibrational control [31], energy harvesting [32], and programmable soft motors [33,34]. These attractive features of MCDEAs can be realized solely by the embedded compliant coupling and simple inputvoltage programming. Despite these promising potential applications, the active dynamic behavior of a MCDEA is poorly understood and difficult to characterize due to the complex interaction between the nonlinear DEA forces, electromechanical coupling, and magnetic repulsion. Thus, a comprehensive study focusing on the active dynamics of an MCDEA that develops an electromechanically coupled dynamic model to characterize the complex actuation behaviors is essential. Such a modeling framework could also benefit the study of other multi-DOF electromechanically coupled systems.
Herein, by adopting the novel MCDEA configuration, we focus on the nonlinear dynamics of this system under parametric excitations to understand how the active dynamics of this nonlinear system can be controlled and exploited. One clear advantage of MCDEAs over other DEA designs is that the two actuation signals, one for each membrane, can be tuned freely to have different voltage amplitudes, phase difference, and frequency difference, which, has not been demonstrated or analyzed systemically in other DEA systems. On the basis of the feature of free input tuning, we conduct a comprehensive study on the effects of the control signals on the two-DOF system. This actuation tuning, electromechanical coupling, together with the compliant interaction between the two membranes is believed to produce new physical insights into a two-DOF DEA system. The rest of this paper is structured as follows. The dynamic model of the MCDEA is developed in Sec. 2 and verified against experimental findings in Sec. 3. Using this model, in Sec. 4, we conduct an in-depth analysis of the nonlinear dynamics of the MCDEA subjected to different actuation signals, including various voltage amplitudes, phase difference, and frequency difference. Finally, in Sec. 5, we discuss and summarize the key findings of this work.

II. DYNAMIC MODELING OF A MCDEA SYSTEM
The thermodynamic equilibrium framework developed by Suo and coworkers [35][36][37][38] is one of the most-widelyused modeling frameworks for analyzing DEA dynamics. Thermodynamic equilibrium requires that, under an isothermal process, the increase in Helmholtz free energy of a DEA should equal the sum of the work done by the external forces and the actuation voltages. The thermodynamic equilibrium approach was widely adopted in studies of DEA configurations such as idealized rectangle DEAs [39][40][41] and pure-shear DEAs [42][43][44][45]. Dissipation in the dielectric elastomer, such as viscoelasticity, can also be considered by use of a nonequilibrium-thermodynamics framework [46][47][48][49]. When the geometry becomes complex and three-dimensional, the strain-stress distribution can be inhomogeneous, and such a framework usually requires a very high computational workload. For dynamic analyses, a simplified kinetic model is typically used to derive the equations of motion of the system. For example, by use of such a kinetic model, the dynamics of a conical DEA [13,18,50,51] were characterized by both modeling and experiments. However, despite the advancements in the dynamic modeling of DEAs, the existing models consider only single-DOF systems. In contrast, the MCDEA, which has two degrees of freedom with a strongly nonlinear coupling mechanism, has not been investigated previously. In this section, we develop a dynamic model that characterizes the complex dynamic response of this MCDEA system.

A. Model overview
As illustrated in Fig. 1(a), the MCDEA has two circular DEA membranes with radius b coupled by the repulsion force from a pair of permanent magnets with radius a. In its passive equilibrium, the magnetic repulsion, F mag , forces the membranes to deform out of plane until it is balanced by the tension of the membranes, F DEA (ignoring gravitational force). The two identical membranes are deformed by d 0 in opposite directions and the distance between the two magnets is s 0 in this equilibrium. When a voltage, , is applied across one membrane, an electrically induced Maxwell pressure causes a force imbalance and the two membranes move toward one side until a new equilibrium state is obtained between the forces F DEA and F mag . Figure 1(b) illustrates the simplified schematic structure of this MCDEA system, where m is the mass of the magnets, d is the displacement, the spring with an arrow represents a nonlinear elastic element, and the dashpot with an arrow represents a nonlinear viscous element. Because of the compliant-coupling mechanism, the deformations of the two membranes d I and d II during this process might not be identical [ Fig. 1(b)] and then the distance between the two magnets s varies. The total out-ofplane deformation of the two membranes is d I = d 0 + d I and d II = −d 0 + d II , respectively (note that an upward displacement is defined as positive).
The three-dimensional geometry, nonlinearity, viscoelasticity, and electromechanical coupling in the system lead to a complex nonlinear dynamic modeling problem. To reduce the complexity of the model, the following assumptions are proposed (after Refs. [13,18,50,51]): (i) the magnets translate only along the vertical axis (i.e., this is a two-DOF system); (ii) the membranes are truncated cone shapes and the circumferential strain is constant; (iii) the radial strain is homogeneous along the radial axis; (iv) the mass of the membrane and the electrode is negligible compared with the mass of the coupling magnet; (v) the electrical response of this system is not considered in this model since the estimated RC constant is less than 10 −4 s [50].

B. Kinetic model
In its reference state, the dielectric elastomer, with initial thickness H 0 , is given a biaxial prestretch of λ p . During an out-of-plane deformation, d n , the radial stretch of the membrane, λ 1_n becomes where n = I or II for the top or bottom membrane respectively.
The angle between the membrane and the horizontal plane during out-of-plane deformation is denoted by α n , where sin α n = d n d 2 n + (b − a) 2 by our assuming the elastomer is incompressible (i.e., fixed volume). The membrane thickness at deformation d n is where λ 2 n is the circumferential stretch and is equal to λ p on the basis of assumption (ii). The equations of motion of the two magnets in the vertical axis yield where m is the mass of the magnet andd n (n = I, II) is the vertical acceleration of the magnets. The vertical force of each DEA membrane, F DEA n , can be presented as where σ 1_n is the true radial stress of membrane n.

C. Material model
The Gent model [52] is used in this work to characterize the strain-stress relationship, and the Kelvin-Voigt viscoelastic model (shown in Fig. 2) is used to describe the viscoelasticity of the elastomer. This viscoelastic model consists of two branches: the first is the nonlinear spring described by the Gent model and the second is the dashpot element, which characterizes the strain-rate-dependent hysteresis in the strain-stress function. The use of the Kelvin-Voigt model has been proven to be valid in characterizing the viscoelasticity of silicone elastomers (see, e.g., Refs. [5,53]) due to the significantly reduced stress relaxation and creep in comparison with the commonly used VHB 4905/4910 tape from 3M [53], where the Maxwell viscoelastic model [37,38] and a fusion of the two (i.e., Kelvin-Voigt-Maxwell model) [51,54] are more popular. The small hysteresis loop in Fig. 3 also indicates low viscoelastic behavior of silicone elastomer. To verify the accuracy of the Kelvin-Voigt model in describing the viscoelasticity of silicone elastomers, a free oscillation experiment is conducted; details of this are described in the next section.
The radial stresses caused by the nonlinear springs, σ 1s_n , are expressed as where μ is the shear modulus of the elastomer, J is the constant of the limiting stretch, E n = n /H n is the electric field, where n is the voltage across the membranes, and 0 and r are the absolute permittivity of a vacuum and the relative permittivity of the dielectric elastomer, respectively. The dashpot in the Kelvin-Voigt model is modeled as a Newtonian fluid (following Refs. [35][36][37][38]). The dashpot in the rheological model (as illustrated in Fig. 2) shares the same strain as the Gent spring, and the radial stresses caused by the deformation of the dashpots are proportional to the rate of deformation and can be described as where η is the viscosity of the dashpot and is greater than zero.
The total radial stresses can be given as

D. Magnetic repulsion model
The magnetic repulsion force, F mag , involves threedimensional interactions of the magnetic fields of the two magnets, which can be extremely complex to model. For example, a cuboidal magnet model describing the magnetic force in the z direction (the same poles facing each other) includes 256 terms [55]. A comprehensive study on the magnetic repulsion of two disk magnets is beyond the scope of this work and, as a result, a simplified magnetic repulsion model [56] is used here to ensure that when the two magnets are infinitely close to each other, the force becomes infinite, and when the distance between the two magnets tends to infinity, the force tends to zero. This model can be written as where k and z are to be determined from experimental testing, and s is the distance between the center of the two magnets.

E. Model summary
With the excitation voltages, I and II , defined, the dynamic response of the MCDEA can be characterized on the basis of the dynamic model, summarized as follows.
The equation of motions of this system are given as F mag = k s z , where n = I or II for the top or bottom membrane, respectively.

III. DYNAMIC MODEL VALIDATION
In this section, the process for fabrication the MCDEA is briefly illustrated and the experimental method for model validation is described. The modeled results are compared with the experimental results to assess the accuracy of the dynamic model.

A. Fabrication and experimental characterization
The MCDEA prototype is fabricated by the following steps. First, silicone elastomer (ELASTOSIL, thickness 100 µm, Wacker Chemie AG) is prestretched biaxially by a factor of 1.2 × 1.2. It is then bonded, with use of silicone transfer tape (ARclear 93495, Adhesives Research), to an acrylic frame with an inner radius of 15 mm. Two 7.5-mm-radius, 0.5-mm-thick disk magnets (0.28 × 2-kg pull force, weight 1.35 g, First4Magnets) are attached to the center of the membrane by the same method. Custom carbon grease [20-wt % carbon black powder (1333-86-4, Cabot Corporation, USA) with 80-wt % vegetable oil] is used as the compliant electrode; in comparison with other formulations, this is found to eliminate the swelling effect, as was found in a previous study [50]. Two DEA frames are connected by bolts and fasteners with 2-mm gaps in between. To compliantly couple the two DEA membranes, the magnets are positioned with the same poles facing each other such that the membranes are deformed out of plane by repulsion.

Quasistatic force-displacement tests
Quasistatic force-displacement tests are performed to identify the quasistatic model parameters of the MCDEA and the magnetic repulsion with the experimental setup shown in Fig. 4(a). The detailed experimental setup is described as follows. The DEA frame is fixed to the testing rig and a linear actuator (X-LSQ150B-E01, ZABER) deforms the center of the DEA membrane out of plane at a low velocity of 0.05 mm/s to ensure negligible viscoelasticity. A constant voltage is generated by a highvoltage amplifier (5HV23-BP1, Ultravolt) and is applied to the DEA during deformation to analyze the effect of the electric field on the force-displacement relationship. The voltage amplitude is determined to be 3.47 kV (equivalent to an electric field of 50 V/µm). A load cell (NO.1004, TEDEA) measures the reaction force of the membrane. The magnetic repulsion is measured in the same way, where one set of magnets is attached to the linear actuator while the other set is fixed to one end of the load cell. With the same poles facing each other, the linear actuator moves one set of magnets relative to the other and the magnetic repulsion is measured by the load cell. All signals are collected by a data-acquisition device (BNC-2111, National Instruments) at a sampling frequency of 5000 Hz and controlled by MATLAB (The MathWorks).

Frequency-sweep tests
Frequency-sweep tests are conducted to investigate the active dynamic performance of the MCDEA, with the experimental setup shown in Fig. 4(b). A sinusoidal voltage signal consisting of an alternating-current (ac) voltage with a biasing direct-current (dc) voltage, I = dc + ac cos I t, is applied to the first DEA membrane via a high-voltage amplifier (5HC23-BP1, Ultravolt), while the second membrane is left to oscillate passively. Two laser displacement sensors measure the displacements of the two membranes at a sampling frequency of 40 000 Hz. The excitation frequency is swept forward from 0 to 120 Hz (120 to 0 Hz for a backward sweep) at a rate of 1 Hz/s, generated by MATLAB using the "chirp" function. ac = dc = 1.74 kV is used in the frequency sweeps. As can be seen, despite the simplifications made in this model, good accuracy is still achieved. Note that four peaks are observed in both the forward sweeps and the backward sweeps, which, have not been observed in previous rigidly coupled DEA systems. The distinguishing four resonant peaks are believed to be the result of the complex interactions between the two DEA membranes and the magnets and are investigated in depth in the next section.

B. Model validation
The value determined for the Kelvin-Voigt model is then verified against the value obtained from a passive freeoscillation test on a different sample, and good agreement is found between the modeling and experimental results, as illustrated in Fig. 6. A passive free-oscillation test is performed to verify the Kelvin-Voigt viscoelastic model by eliminating the effects of magnetic interaction and parametric forcing by the voltage inputs. A single-cone DEA is fixed to the testing rig, and the central disk is deformed out of plane by 4 mm by use of a string and then fixed. Differing from the quasistatic and frequency-sweep tests, the DEA dimensions here are an frame inner radius of 0 mm, a central disk radius of 4 mm and mass of 0.13 g, and a membrane thickness of 50 µm. The same carbon electrode is applied to the membrane. This passive oscillation test aims at validating the accuracy of the Kelvin-Voigt viscoelastic model with silicone materials by use of the model parameters determined from frequency-sweep tests (described above) to predict the free-oscillation response of a different sample. After the out-of-plane deformation of the DEA is fixed for a period of 120 s (to eliminate the effects of stress relaxation), the string is cut and the central mass, together with the silicone membrane, is released and allowed to oscillate freely. The same laser displacement sensor is used to measure the passive step response at a sampling frequency of 20 000 Hz. An illustration of the setup of the free oscillation test is shown in

IV. DYNAMIC STUDY OF A MCDEA SYSTEM UNDER PARAMETRIC EXCITATION
In the previous section, the proposed numerical model was validated against experimental results with excellent accuracy. The dynamic response of the MCDEA under an ac actuation voltage was obtained by our computing the time-series response of the system using a direct numerical simulation, sometimes referred to as a "brute-force numerical approach." Despite its ease of use, such a method can be computationally costly, especially when short time steps are required, and several parameters are varied. Additionally, this approach is unable to measure unstable responses 044033-7 or detect bifurcations. Hence, in this section, before conducting the comprehensive dynamic analysis of MCDEAs, we begin by introducing an advanced numerical simulation framework using the MATLAB-based continuation package COMPUTATIONAL CONTINUATION CORE (COCO) (a more detailed introduction to this software package can be found in Ref. [57]). With the same numerical model and parameters as developed in Sec. 3, COCO can predict a range of dynamic phenomena, including bifurcations and unstable solutions.
One of the most significant advantages of the MCDEAs is that the emerging complex active dynamic phenomena can be controlled by simple voltage programming. To understand this principle in depth, in this section, by using the electromechanically coupled dynamic model developed, we analyze the active dynamic responses of the MCDEA system with the following four input-voltage cases: (a) I = dc + ac , II = 0. This is the fundamental case, where one membrane is actuated by an alternatingcurrent voltage with a biasing direct-current voltage, while the other membrane remains passive.
(b) I = dc + ac with various ac / dc ratios, II = 0. This case investigates the effects of the amplitudes of the voltage components (i.e., ac and dc biasing-voltage components) on the dynamic behaviors of the MCDEA. Cases (a) and (b) with one membrane remaining passive represent an "isolated" actuation strategy (i.e., the passive membrane may act as a galvanically isolated end effector since energy is transferred from the active membrane via contactless magnetic repulsion force). Here the stroke is controllable via the excitation frequency and ac / dc ratio of the active membrane.
(c) I = dc + ac , II = dc + ac with equal excitation frequencies (i.e., I = II = ) but with a different phase difference between the ac signals. As highlighted in Sec. 1, a distinguishing feature of the compliantly coupled MCDEA is that the two input signals can be tuned freely, which can lead to interesting dynamic phenomena. Hence, in this case, we investigate the effects of two input signals with the same frequency but different phases, between 0 and 2π . The phase tuning can be a potential alternative dynamic control strategy in addition to conventional voltage-amplitude control [27,29].
(d) I = dc + ac , II = dc + ac with different excitation frequencies (i.e., I = II ), which causes the emergent beat phenomena of the MCDEA. Such phenomena can be useful in acoustics and enable the potential for DEA-driven soft loudspeakers and active vibrational control.

A. Case (a): I = dc + ac , II = 0
In this first case, only one DEA membrane is actuated, while the other membrane can oscillate passively, as in the dynamic study in Sec. 3. A sinusoidal actuation voltage with the same dc and ac amplitudes as used in Sec. 3 ( ac = dc = 1.74 kV) is used here.
The frequency-domain response of the MCDEA simulated with COCO is shown in Fig. 5, where |X I | is the oscillation amplitude of d I . It is noteworthy that, along with the four resonant peaks observed in experiments and brute-force numerical simulations, an additional peak near 124 Hz is predicted in the COCO simulation. The blue dots at the root of the high-amplitude branch of peak 5 and the region highlighted by the blue box in Fig. 7 represent the period-doubling bifurcations where the system switches to an emerging behavior with the oscillation period twice that of the original response. These period-doubled responses are not observed during the relatively fast frequency sweep [1 Hz/s in Figs. 4(c)-4(f), where a steady-state response is never reached] but are able to be captured by COCO, showing the clear advantage of using this simulation tool for DEA dynamic studies.
For peak 4, the resonant peak is heavily distorted to the right, leading to two fold bifurcations, denoted by two red dots. For excitation frequencies between points A and B, three possible periodic solutions exist, in which the two solid curves represent stable solutions and the dashed curve indicates an unstable solution. As the excitation frequency increases from below the resonance (e.g., from below 80 Hz), the response of the system follows the upper stable-solution path with increasing amplitude. When bifurcation point A is reached, the response jumps 044033-8 down to the other stable solution with lower amplitude. In contrast, when the excitation frequency is decreased from above the resonance (e.g., from more than 100 Hz), the response follows the lower stable branch until reaching bifurcation point B, where the response jumps up to the upper stable branch, as illustrated by the enlarged section highlighted by the red box in Fig. 7. The region between the two points is referred to as a "hysteresis region." The detailed steady-state time series, phase portraits and Poincaré maps, and discrete Fourier transforms (DFTs) of the response for the five resonant peaks marked in Fig. 7 are shown in Fig. 8. First, from the time-series results [Figs. 8(a)-8(e)], it can be noted that the two outputs are in phase for peaks 1, 3, and 5 and are antiphase for peaks 2 and 4. These correspond to the responses of the two underlying linear modes, as is commonly found in two-DOF oscillatory systems [58]. When a component of the excitation frequency is close to the natural frequency of one of the linear modes, the response of the system will be dominated by this specific mode.
Second, in can be seen from both the time-series results and the DFTs that for peaks 1 and 2 the two outputs have a fundamental response frequency (which is defined as the frequency of the component with the largest amplitude in a DFT), ω n , twice that of the driving frequency, I , demonstrating a superharmonic response. For peaks 3 and 4, where the oscillation amplitudes are higher, the fundamental response frequency equals the driving frequency, showing a harmonic resonance. For peak 5, the fundamental response frequency is half the driving frequency, showing a clear subharmonic resonance.
The phase portraits and Poincaré maps in 50 cycles overlap for all five peaks [Figs. 8(f)-8(j)], demonstrating steady periodic oscillations of the MCDEA. Two Poincaré points emerge for peak 5 [ Fig. 8(j)], which is due to the period doubling of the system.
To understand the occurrence of multiple resonant peaks in this system, recall [from Eq. (5)] that the excitation force is given by p = ε 0 ε r ( 2 /H 2 ), where p is the Maxwell pressure (i.e., the force experienced by the mass is a function of the square of the driving voltage signal, , which also adds nonlinearity to this MCDEA system). Therefore, recalling that the voltage signal for the first mass is given by I = dc + ac cos I t, we find the force experienced by the first mass is proportional to where This demonstrates that the forces experienced by the mass contain two time-dependent components: one at frequency I and the other at frequency 2 I . As a result, when the excitation frequency, I , is close to either r /2 or r (where r is one of the resonant frequencies of this system), resonance could occur. The superharmonic response (corresponding to the case where I is close to r /2) of this MCDEA system can be potentially useful for engineering systems where high-frequency oscillations can be generated by lower-frequency signals, which could increase the energy efficiency of the systems.
It also can be noted from the Maxwell pressure equation that the excitation of the DEA is a function not only of the actuation voltage but also of the displacement of the DEAs [recall from Eqs. (1) and (2), H = H 0 /λ 1 λ 2 , and λ 1 = f (d)]; hence, this system is excited parametrically. Subharmonic resonances (such as peak 5) are very common in parametrically excited systems, as demonstrated in Refs. [59][60][61][62].

B. Case (b): I = dc + ac with different ac / dc ratio, II = 0
In this study, we investigate the effects of different ac / dc ratio on the dynamic responses of the MCDEA by fixing the biasing excitation term, E a , from Eqs. (9) and (10) (i.e.,E a = 4.52 kV 2 ) [note the value is obtained from case (a), where ac = dc = 1.74 kV]. The role of 0 ≤ ac / dc ≤ 1 is followed. One DEA membrane is excited by the ac voltage with various ac / dc ratios, while the other membrane remains unactuated but allows passive oscillation. The voltage-tuning strategy used in this study differs from that used in previous studies; that is, instead of fixing the biasing dc voltage, dc , and varying the ac component, ac , we fix the biasing excitation term here. As the biasing excitation term is a product of both dc and ac , fixing dc and varying ac can cause a change in the biasing excitation term, which leads to a shift in the resonant behavior. Tuning the ac / dc ratio while keeping the biasing excitation term constant can ensure a more-precise control of the amplitude with a fixed resonant behavior. Figure 9 shows the frequency response of the MCDEA with ac / dc ratio ranging from 0 to 1. It can be clearly noticed that as the ac / dc ratio increases, the amplitudes of the resonant peaks increase. When the ac / dc ratio is greater than 0.8, a new subharmonic peak (peak 5) emerges, which is due to the parametric excitation of this system (also shown in Fig. 10). As can be seen from Fig. 10, at low ac / dc ratios, peak 4 does not exhibit a bifurcation (i.e., the response follows a single stable path); however, as the ac / dc ratio increases to approximately 0.3, bifurcations occur, leading to the hysteresis found in case (a). This is due to the lower amplitudes that are achieved when the ac component of the excitation is lower.

044033-9
(a) (f) (k)  Figure 11 shows the amplitudes of resonant peaks 1-4 against the ac / dc ratio. The amplitudes of the I excitation component, E b , and 2 I excitation component, E c , are also shown in Figs. 11(a) and 11(b), respectively. The increase of oscillation amplitudes of peaks 1 and 2 follows the same trend as their excitation amplitude, E c . The same is true for peaks 3 and 4 and E b . The greater amplitude of E b results in larger oscillation amplitudes of peaks 3 and 4. In Fig. 11(b), the curve for peak 4 is branched as the ac / dc ratio increases, which represents the two possible stable solutions.  C. Case (c): I = dc + ac , II = dc + ac , different relative phases In the last two case studies, only one DEA membrane is actuated, while the other membrane oscillates passively. In this case, both DEA membranes are actuated by an ac voltage signal with the same excitation frequency (i.e., I = II = ) but different relative phases. Conventional rigidly coupled, antagonistic DEAs allow only antiphase voltage signals to achieve bidirectional actuation. However, because of their compliant coupling, MCDEAs allow voltage signals with a relative phase ranging from 0 to 2π . Different relative phase values lead to emerging dynamic responses that were not observed in previous DEA studies [27].
The square of the two actuation voltage signals in this case study is written as where θ is the phase difference between I and II and is within the range from 0 to 2π and E a = 4.54 kV 2 , E b = 3.02 kV 2 , and E c = 1.51 kV 2 , following case (a). The frequency responses of the MCDEA driven by two actuation signals, for different relative phases, θ, are shown in Fig. 12. It can be seen that the change in phase not only affects the amplitudes of the resonant peaks but can also determine the existence of all the superharmonic, primaryharmonic, and subharmonic resonances. For example, at θ = 0, peaks 1 and 3 (i.e., the superharmonic and primaryharmonic resonances of the first mode, where d I and d II oscillate in phase) completely vanish, while peaks 2 and 4 (i.e., the superharmonic and primary-harmonic resonances of the second mode, where d I and d II oscillate in antiphase) reach their maximum amplitudes. This can be explained by the fact that the and 2 components of the forcing ( 2 I and 2 II ) are both multiplied by coefficients that are opposite in sign. As a result, when θ = 0, the excitation forces 044033-11 are in antiphase and hence only excite the antiphase resonances. When θ = π /2, the phase of the 2 component of the excitation is 2θ = π [Eqs. (13) and (14)], which leads to an in-phase excitation, and hence the first superharmonic resonance, peak 1 (exhibiting an in-phase response), reaches its maximum amplitude. Similarly, the case θ = π strongly excites the in-phase primary response, peak 3, and the antiphase superharmonic resonance, peak 2 (as the phase of the 2 component is 2θ = 2π , which generates an antiphase excitation force). A new subharmonic resonance, peak 6, emerges near 170 Hz and reaches its highest amplitude when θ = 0 (where peak 5 also reaches its maximum). The two subharmonic peaks 5 and 6 are not triggered at θ = π .
Detailed oscillation amplitudes of all six resonant peaks as a function of the relative phase difference are shown in Fig. 13. For peaks 1 and 2, which are the superharmonic resonances, the amplitude-phase response repeats itself when θ varies from 0 to 2π , since they are excited by the 2 excitation component. For primary harmonic peaks 3 and 4, which are driven by the excitation component, the response repeats only once between 0 and 2π . For subharmonic peaks 5 and 6, high amplitudes are achieved only when θ is close to 0 and reduce rapidly as θ increases.
The results are validated against experimental results by forward sweeps from 0 to 120 Hz at 1 Hz/s, as shown in Fig. 14. Brute-force numerical simulation is used as a comparison. The experimental results agree well with the model simulation and the occurrences of the first four resonant peaks are as demonstrated in the numerical analysis above. A resonant peak (peak 4) with a lower amplitude is observed in the experiment at θ = 0, which could be due to the reduced stability caused by the slight asymmetry in the two masses and the two DEA membranes (i.e., prestretch, RC constant). The two subharmonic peaks 5 and 6 are experimentally observed by a slow-rate frequency sweep from 110 to 130 Hz and from 150 to 170 Hz, respectively, at 0.02 Hz/s with θ = 0, where the two peaks reach their maximum [Figs. 14(g) and 14(h)]. The measured period-doubling bifurcation points (120.6 and 164.8 Hz) are very close to the model prediction (120.5 and 165 Hz). However, the experimental results do not follow the highamplitude branch until their maximum, which could also be due to the reduced stability as described above. It is noteworthy that the period-doubling bifurcation occurred in this case differs from that in case (a) in the "delayed" transition. Figure 15(a) shows the simulated transient dynamic response near the first period-doubling bifurcation for peak 5 with θ = 0 and /2π = 124 Hz. Even though the two DEA membranes are actuated with two signals at /2π = 124 Hz when t = 0, a perioddoubling bifurcation does not occur until t ∼ 8 s. Between t = 0 s and t = 8 s, the fundamental response frequency equals the driving frequency. Only from t = 8 and 9 s, does the fundamental response frequency gradually decay to ω I = 1/2 , as can been seen from the enlarged subplots in Fig. 15(a). The oscillation amplitude increases dramatically, overshoots, and reaches a steady state with an oscillation frequency half the driving frequency at t = 11 s. This transition process is also demonstrated in the phase portraits and Poincaré maps in Fig. 15(b), where the Poincaré points are first concentrated in the center (before bifurcation) and then spread out toward the two sides due to period doubling and finally overlap, demonstrating a highamplitude steady-state response after the period-doubling bifurcation.
The phase tuning demonstrated in this case represents an approach beyond voltage amplitude tuning to control the dynamic behavior of a two-DOF DEA. Adjusting the phase difference between the two actuation signals controls not only the amplitudes but also the occurrence of a specific resonant peak. This phase-tuning strategy is 044033-13 proposed to be advantageous for future vibration-control and energy-harvesting applications, where a specific mode of oscillation and resonance is required.

D. Case (d): I = dc + ac , II = dc + ac , different frequencies
In the final case, we examine the dynamics of the MCDEA driven by two ac voltage signals with the same voltage amplitudes [ ac = dc = 1.74 kV, following case (a)] but different excitation frequencies (i.e., I = II ). COCO is ill-suited to such simulations due to its requirement for responses to be of a finite period (which is not met when the ratio between the excitation frequencies is irrational). Because of this, brute-force numerical simulation is used instead.
In the first study, we simulate the steady-state response of the MCDEA with II /2π fixed at 10 Hz while I /2π varies from 9.5 to 10.5 Hz in increments of 0.2 Hz. The steady-state results obtained with six sets of frequency combinations are shown in Fig. 16. Twenty seconds is allowed for the system to reach a steady state in the simulation. The system demonstrates a clear beating phenomenon under two different actuation frequencies. The envelope of the displacement d I shows a periodic behavior with frequency equal to the absolute value of the frequency difference | I -II |. The response frequencies of d I and d II contain a mixture of I and II . For example, Fig. 17 shows the fast-Fourier-transform result for d I and d II with excitation frequencies I /2π = 9.5 Hz and II /2π = 10 Hz. It can be seen that d I has a strong frequency component of its excitation frequency of 9.5 Hz, but also a weaker frequency component from the other excitation frequency of 10 Hz due to the interaction of the two masses. The same is true for d II .
In the study described above, the two actuation frequencies are chosen to be away from any resonant frequencies.
A beating phenomenon can also be observed when the actuation frequencies are close to the resonance. Figure 18  The frequency difference for all three cases is set at 0.5 Hz and, as can be seen in Fig. 18, the displacements in all three cases demonstrate a beat frequency of 0.5 Hz. It is noteworthy that as the actuation frequencies are close to the resonances, the beat amplitudes are significantly increased, as shown in Figs. 18(b) and 18(c) and the envelopes are less harmonic than in Fig. 18 Figure 19 shows the measured beat phenomena with frequency differences of 1, 0.2, and 0.1 Hz. In all these cases, the excitation frequencies are close to the primary resonance of the first mode, and hence the two displacements are in phase.
In this study we demonstrate that by introducing a frequency difference in the actuation voltage signals a beating phenomenon can be achieved where the beat frequency is simply the difference between the two actuation frequencies. In practical applications, the high programmability of the beat phenomenon of the MCDEA can be desirable, where the beat frequency, the beat amplitude, and the mode can be controlled separately. This is because the beat frequency is determined solely by the voltage-frequency difference, while the beat amplitude can be tuned by controlling the actuation frequencies to be close to, or away from, resonances, and the beat mode (in phase or antiphase) can be determined by having the voltage frequencies matching the resonance with the desired mode.

V. DISCUSSION AND CONCLUSION
The magnetically coupled double-cone DEA shown in this work represents an emerging two-DOF nonlinear system that demonstrates rich nonlinear dynamic behavior and promising potential applications in robotics, energy harvesting, and smart structures. However, the threedimensional geometry, nonlinearity in the coupling mechanism, electromechanical coupling, and the elastomers in the system lead to a complex nonlinear dynamic modeling problem. A numerical model is developed in this work to characterize its active dynamic response and it is verified against experimental results with excellent accuracy. One of the most significant advantages of MCDEAs is the emerging complex active dynamic phenomena can be controlled by simple voltage programming, which features the advantage of an advanced soft active material. With the model developed, the dynamics of this system are characterized for four different actuation cases: (i) only the first membrane is excited; that is, I = dc + ac , II = 0; (ii) I = dc + ac with different ac / dc ratio, II = 0; (iii) both membranes are excited; that is, I = dc + ac , II = dc + ac of the same frequency as I but different relative phases; and (iv) both membranes are excited; that is, I = dc + ac , II = dc + ac with different frequencies. The key findings from these four case studies can be summarized as follows: (a) With only one membrane actuated, this MCDEA system exhibits multiple complex nonlinear dynamic phenomena, including superharmonic, primary-harmonic, and subharmonic resonances, bifurcations and multiple oscillation modes (in phase and antiphase).
(b) By analyzing the dynamics of the MCDEA with different ac / dc component ratio of the actuation voltage, while keeping the biasing excitation amplitude constant, we find the oscillation amplitudes of the resonance peaks are strongly correlated to the corresponding excitation amplitudes.
(c) Adjusting the phase difference between the two actuation signals not only controls the amplitudes of the superharmonic, primary-harmonic, and subharmonic resonances but can directly determine the existence of a specific resonance. This finding suggests that phase tuning can offer a simple yet powerful control of the dynamic response of a two-DOF DEA system, along with the conventional voltage-amplitude adjustment.
(d) Apart from voltage regulation and phase tuning, the dynamics of the proposed MCDEA system can be controlled by varying the frequency difference between the two voltages. Beat phenomena, with high programmability, are demonstrated in this study. The beat frequency can be determined solely by the voltage-frequency difference, the beat amplitude can be adjusted by tuning the actuation frequencies to be close to, or away from, resonances, and the beat mode can be determined by letting the voltage frequencies match the resonance with the desired mode.
The dynamic model developed in this work offers a simple yet powerful tool for the analysis of nonlinear DEA systems and can help in the design and optimization of DEAs in dynamic applications. The voltage-control strategies demonstrated in this work, including phase tuning and frequency control, generate multiple controllable dynamic behaviors, including controlled resonance amplitude, mode, and frequency response (i.e., control of the relationship between the driving frequency and the response frequency) and beat phenomena. These emerging dynamic phenomena generated by these control strategies are believed to be useful in robotics, such as for active vibrational control, and energy harvesting. This design can potentially be used as a highly programmable soft motor in robotic locomotion to enable multiple locomotion modes by a simple voltage control of the actuator rather than reprogramming the whole robot structure.