Dynamic modeling of dielectric elastomer actuator with conical shape

With desirable physical performances of impressive actuation strain, high energy density, high degree of electromechanical coupling and high mechanical compliance, dielectric elastomer actuators (DEAs) are widely employed to actuate the soft robots. However, there are many challenges to establish the dynamic models for DEAs, such as their inherent nonlinearity, complex electromechanical coupling, and time-dependent viscoelastic behavior. Moreover, most previous studies concentrated on the planar DEAs, but the studies on DEAs with some other functional shapes are insufficient. In this paper, by investigating a conical DEA with the material of polydimethylsiloxane and considering the influence of inertia, we propose a dynamic model based on the principles of nonequilibrium thermodynamics. This dynamic model can describe the complex motion characteristics of the conical DEA. Based on the experimental data, the differential evolution algorithm is employed to identify the undetermined parameters of the developed dynamic model. The result of the model validation demonstrates the effectiveness of the model.


Introduction
Soft robots, a kind of flexible machinery, aim at operating in natural environments and realizing complex functions [1]. Although conventional rigid robots have made great progress in the field of automation manufacturing, soft robots are more flexible and provide great potential applications [2]. In addition, soft robots are mostly made of soft materials, and they are capable of deforming greatly and rather adopt to the complex external environments [3,4].
Traditional robots usually take electric motors, hydraulic motors and cylinders as their actuators. However, soft robots mostly employ the flexible actuators fabricated by soft materials [5]. The pneumatic actuator is a typical flexible actuator. In [6], a soft gripper is fabricated by full multi material 3D printing technology, which can freely deform and grip various objects. Moreover, a climbing robot designed in [7] is capable of performing 3D climbing locomotion using two suction cups. Different from the pneumatic actuator, the soft actuator based on smart materials is another typically flexible actuator. In [8], a jellyfish robot powered by the voltage, whose amplitudes and frequencies can be set as different values within one period, is applied to the electrodes. Based on the experimental data, the undetermined parameters in the dynamic model are identified by the differential evolution algorithm. The model validation indicates that the dynamic model contributes to describing the viscoelastic behavior and electromechanical response of the conical DEA.

DEA modeling
In this section, a dynamic model of a DEA with conical shape is developed. For ease of presentation, three different states of the DEA are declared in advance. The first state is called undeformed state, the second state is called pre-stretched state, and the third state is called electro-deformed state, whose diagrams are shown in Fig 1(a), 1(b) and 1(c), respectively.

(A) Un-deformed state
A DE membrane with thickness d 0 is clamped by a frame with inner circle radius R. A loadbearing plate with radius R 0 is placed on the center of the DE membrane. Two sides of the DE membrane, which are two annular regions, are coated with the compliant electrodes. Thus, the radial length of the DEA is L 0 = R − R 0 .

(B) Pre-stretched state
A weight with the mass m is placed on the center of the load-bearing plate. Subjected to the gravity P, the weight will move down a distance z 1 to reach the equilibrium position. As a result, the DE membrane is pre-stretched as a conical shape. As shown in Fig 1(b), L 1 , d 1 and h 1 are the dimensions of the DEA corresponding to the pre-stretched state, where L 1 is the generatrix length, d 1 is the thickness, and h 1 is the height difference between the upper surface and the lower surface.

(C) Electro-deformed state
When a driving voltage F is applied to the electrodes, the DE membrane reduces in thickness and expands in area. Thus, the weight will move down a displacement z 2 . As shown in Fig 1(c), L 2 , d 2 and h 2 are the dimensions corresponding to the electro-deformed state.

PLOS ONE
Dynamic modeling of dielectric elastomer actuator The volumes of the DEA for the un-deformed state, the pre-stretched state and the electrodeformed state are: Strictly speaking, the deformation of the DEA with conical shape is inhomogeneous [26,27]. However, to simplify the dynamic modeling, the inhomogeneity of the deformation is ignored in the following development [24,28]. Since the DEA is incompressible [29], the volume of the DEA is constant. Thus, V 0 = V 1 = V 2 . From (1), we can get According to (2), the relationships among z 1 , z 2 , d 1 and d 2 are ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi The DEA studied in this paper is the conical shape. For ease of description, the generatrix, thickness and circumferential stretches are employed to describe the states of the DEA. In the pre-stretched state, the pre-stretches of the DEA are λ pre,L , λ pre,d and λ pre,C , respectively. In the electro-deformed state, the stretches of the DEA are λ 1 , λ 2 , and λ 3 , respectively. According to Fig 1, the following equations hold: According to (2)-(5), the following equation is established: The relationship between the charge Q and the voltage F is where ε and C are the permittivity and the capacitance of the DE material, respectively. According to (3)-(6), the relationship between δλ 1 and δz 2 is dz 2 dl 1 ¼ L 2 L 0 ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi From (6) and (7), the charge on the electrode varies by To calculate the work of the inertial forces during the electromechanical deformation, we consider the cylindrical coordinates shown in Fig 2, where O, r, φ and z represent the coordinate origin, the radial distance, the azimuth angle and the height of the cylindrical coordinates, respectively.
As shown in Fig 2(b), an infinitesimal element with inner radius r 1 and outside radius r 1 + dr 1 is investigated. In the electro-deformed state, the displacement of the element along the r-direction, φ-direction, and z-direction are 0, 0 and z r1 , respectively. So, the relationship between z r1 and z 2 is The inertial forces in each material element along the r-direction, φ-direction, and z-direction are 0, 0 and dF r1 , respectively. According to D'Alembert's principle, we can get where ρ is the density of the DE material. Thus, the changes of works done by the inertial forces are 0, 0 and δH I,z , respectively. According to (10) and (11), the work done by the inertial force dF r1 is The change of the free energy of the DEA is equal to the sum of the works done by the driving voltage, the gravity and the inertial forces. That is, where W is the free energy density of the DEA, and δW represents the change of W. By submitting (9) and (12) into (13), the free energy density W varies by Submitting (8) into (14), we can get where In order to describe the viscoelasticity of the DE material, the rheological model with two parallel units (as shown in Fig 3) is adopted [30]. The part A only consists of a spring α 0 , while the part B consists of four parallel formations and each formation consists of a spring α i (i = 1, 2, 3, 4, . . ., n) with a series-wound dashpot. In this paper, we suppose each dashpot to be a Newtonian fluid with viscosity η i . Let ξ ij (j = 1, 2) be the stretches due to the dashpots, the stretches of the spring α i are determined by multiplication rules l e i1 ¼ l 1 =x i1 , and l e i2 ¼ l 2 =x i2 ¼ l À 1 1 x À 1 i2 . The free energy density W of the DEA [29] can be described as where W s is the Helmholtz free energy associated with the stretching of the elastomer, and D is the electric displacement. The electric displacement D is equal to In this paper, we choose Gent model [17], [30] to describe the elastic energy density of the DEA. Therefore, the elastic energy density of the DEA is where W a i ela are elastic energy densities of the spring α i ; μ i are shear modulus of the spring α i , respectively; J i are deformation limits of the spring α i , respectively.
According to (5)- (7) and (17)- (19), the free energy density of the DEA is According to Newton's third law of motion, the stresses of the spring α i (i = 1, 2, 3, 4, . . ., n) are equal to the corresponding stresses of the dashpot. So, From (19) and (21), the strain rates of the dashpots can be expressed as The viscoelastic relaxation time T i (i = 1, 2, . . ., n) of the DEA is defined as the ratio of η i to μ i . So, Submitting (20) into (15), and combining the result with (22), the dynamic model of the conical DEA can be described as ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi So far, we have developed a dynamic model to describe the inherent nonlinearity, complex electromechanical coupling and time-dependent viscoelastic behavior of the conical DEA. In the following works, we conduct experiment to acquire experimental data of the conical DEA, and then employ differential evolution algorithm to identify the undetermined parameters in (24) based on these data.

System description
In this section, we first introduce the fabrication of the conical DEA briefly. Then, the experimental platform is described. It is worth pointing out that the DE film has been wrinkled before applying voltage, and such wrinkles are difficult to be eliminated completely. However, we repeatedly adjusted the pose of the load-bearing plate to eliminate the wrinkles of the DE film as much as possible. Meanwhile, the DEA was left standing for a long time to make the wrinkles eliminate gradually before each experiment. Through the above measures, the DE film has fewer wrinkles in the experiment, the effect of wrinkling on the electrical deformation is minimized, and the availabilities of experimental results are ensured.

Experimental platform
The experimental platform (see Fig 5)  The function of the I/O module is to output an original voltage signal for the high voltage amplifier, and capture the real-time displacement data from the laser sensor. The high voltage amplifier is used to amplify the original voltage signal by 1000 times and apply it to the electrodes of the DEA.

Model identification
In this section, we first introduce the driving voltage applied in the experiment. Then, the undetermined parameters are identified based on the differential evolution algorithm. Considering the precision and the hardware capabilities, in the dynamic model (24), four springdashpot units are employed to describe the viscoelasticity of the DEA.

Driving voltage
To facilitate the acquisitions of the experimental data, the following driving voltage is applied.
where a i is amplitude; f i is the frequency; t is the time; rem(α/β) is the remainder of α divided by β. By letting t m ¼ remðt; 1=f i Þ, the periodic driving voltage in t 2 [0, + 1) is generated.
By setting different values of a i and f i , the driving voltages with different amplitudes and different frequencies are generated within one period.

Parameters identification
In the pre-stretched state, the vertial displacement of the weight is measured to be z 1 = 1.26 (cm). The sampling period of the experiment is set as T = 0.01 (s). When a i = 5.5+ 0.5i (kV) (i = 1, 2, . . ., 5) and f i = 0.2i (Hz), the diagram of the driving voltage is shown in Fig 6. To avoid the negative displacement of the weight, the maximum frequency of the driving voltage is limited to 1.0 (Hz) in all experiments [31]. The differential evolution algorithm for the parameters identification is briefly listed in the Fig 7. Considering that we do not have any prior knowledge about the values of J i , μ i and T i , we set the large enough search ranges to ensure that the differential evolution algorithm could find out the optimal solution. That is, the search range of J i is (0, 9 × 10 8 ], the search range of μ i is (0, 8 × 10 6 ] and the search range of T i is (0, 3 × 10 6 ]. For conveniently describing the performance of the model prediction, the root-meansquare error e rms and the maximum tracking error e m are introduced. e rms ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi 1 n where z ei and z mi represent the experimental data and the model predicted value of the displacement in the vertical direction; n is the sampling quantity within the sampling time.  Table 1 lists the identified parameters of the dynamic model (24). The root-mean-square error e rms is 0.69% and the maximum tracking error e m is 1.60%.

Model validation
The input of the dynamic model (24) is the voltage shown in (25). By setting different values of a i and f i , in this section, the generalization ability of the proposed dynamic model of the conical DEA is validated.

Model validation with different driving voltage amplitudes
The amplitudes of the driving voltage are set to be a i = 5.5+ 0.5i (kV) (i = 1, 2, . . ., 5). Moreover, the frequencies are set to be f i = 0.2, 0.4, 0.6, 0.8, 1.0 (Hz), respectively. So, the driving voltage has various amplitudes but single frequency in each test experiment.
Applied by the driving voltage with single frequency and multi amplitudes, the comparisons of the model prediction and the experimental result in each test experiment are shown in Fig  10. The modeling error for all test experiments are shown in Table 2. According to the above results, the root-mean-square error of the modeling for any test experiment is less than 3%, and the maximum modeling error for any test experiment is less than 6%. Therefore, the generalization ability of the proposed dynamic model of the DEA is fairly good. https://doi.org/10.1371/journal.pone.0235229.g008    Table 3.
According to the above results, the root-mean-square error of the modeling for any test experiment is less than 2%, and the maximum modeling error for any test experiment is less than 3%. Therefore, the developed dynamic model has excellent performance in the generalization ability.

Model validation corresponding to force analysis
To further verify the validation of the proposed model, the force versus displacement and force versus voltage tests for the dynamic response are performed. The amplitudes and the frequencies of the driving voltage are chosen to be a i = 5.5+ 0.5i(i = 1, 2, . . ., 5) (kV) and f i = 0.2i (Hz), respectively. Based on the real-time displacement data measured by the laser sensor, the accelerated velocity of the weight is calculated by adopting the third-order differentiator. Thus, the output force of the DEA for the experiment can be obtained according to Newton's second law. Moreover, the output force with respect to the model prediction can be calculated according to the proposed model (24).
In this way, the comparisons of the model prediction and the experimental result corresponding to the force versus time, force versus displacement and force versus voltage are given in Fig 12. The root-mean-square error and the maximum tracking error for all tests are 0.0028% and 6.0011%, respectively. Therefore, the validation of the proposed dynamic model is further verified.
In the above works, we verify the validity of the model driving by the voltage with different amplitudes and frequencies, respectively. Meanwhile, the force versus displacement and force versus voltage analyses are conducted. According to the comparison results, the proposed dynamic model is valid.
Next, to further reflect the value of the model, the amplitude-frequency response analysis is developed. The sinusoidal voltages with frequencies 0.01 Hz to 10 Hz (spacing 0.01 Hz) are employed in the theoretical calculations. The amplitude-frequency response curve is shown in Fig 13. With the increase of the frequency of the driving voltage, the amplitude of the conical DEA reduces continuously. This may originate from the viscoelasticity of the DE material [23].

Conclusion
In this paper, the dynamic model of the conical DEA is proposed based on the theory of nonequilibrium thermodynamics. First, three different states of the DEA are declared and its deformation mechanism is analyzed. Then, the infinitesimal element with conical shape in cylindrical coordinates is used to calculate the work done by the inertial force. To describe the elastic energy and the viscoelasticity of the DEA, Gent model and the rheological model are employed respectively. Next, the undetermined parameters in the dynamic model of the DEA are identified by using the differential evolution algorithm. Finally, the comparisons of the experimental result and the model prediction output demonstrate that the proposed dynamic model can describe the inherent nonlinearity, complex electromechanical coupling and timedependent viscoelastic behavior of the conical DEA. In addition, we find that the DEA shows the obvious hysteresis behavior, creep behavior, and even rate-dependence hysteresis behavior during the experiments. The proposed model can still handle the above behaviors accurately. So, the dynamic model contributes to understanding the complex motion characteristics of the conical DEA.