A Computing Method to Determine the Performance of an Ionic Liquid Gel Soft Actuator

A new type of soft actuator material—an ionic liquid gel (ILG) that consists of BMIMBF4, HEMA, DEAP, and ZrO2—is polymerized into a gel state under ultraviolet (UV) light irradiation. In this paper, we first propose that the ILG conforms to the assumptions of hyperelastic theory and that the Mooney-Rivlin model can be used to study the properties of the ILG. Under the five-parameter and nine-parameter Mooney-Rivlin models, the formulas for the calculation of the uniaxial tensile stress, plane uniform tensile stress, and 3D directional stress are deduced. The five-parameter and nine-parameter Mooney-Rivlin models of the ILG with a ZrO2 content of 3 wt% were obtained by uniaxial tensile testing, and the parameters are denoted as c10, c01, c20, c11, and c02 and c10, c01, c20, c11, c02, c30, c21, c12, and c03, respectively. Through the analysis and comparison of the uniaxial tensile stress between the calculated and experimental data, the error between the stress data calculated from the five-parameter Mooney-Rivlin model and the experimental data is less than 0.51%, and the error between the stress data calculated from the nine-parameter Mooney-Rivlin model and the experimental data is no more than 8.87%. Hence, our work presents a feasible and credible formula for the calculation of the stress of the ILG. This work opens a new path to assess the performance of a soft actuator composed of an ILG and will contribute to the optimized design of soft robots.


Introduction
At room temperature, the ionogel is a polymerized gelatinous mixture from the ionic liquid and polymer matrix under UV irradiation [1,2]. The high conductivity and stability of ionic liquids and the good mechanical properties of polymers make ionic liquid gels an ideal replacement for the traditional electroactive polymers [3]. Due to their high environmental adaptability and low-pressure impedance characteristics, soft robots have shown broad application potential in the fields of biology, medicine, agriculture, and so on. The adoption of electroactive polymer (EAP) materials for soft robots has become a hot research topic in recent years. Progress in electrochemical actuators has been made over the past few decades due to their desirable mechanical properties for intelligent robots, which are an alternative for air-and fluid-derived equipment [4][5][6][7]. Ionic liquid gels (ILG) are suitable building blocks for advanced actuators due to their tunable ionic conductivity, chemical stability, thermal stability, and simple ion transport [8][9].
Electrochemical actuators have been further developed over the past few decades for their desirable mechanical properties in intelligent robots, which are alternatives to air-and fluid-derived devices. The flexible ionic conductivity of ILG is better suited to the evolution of building blocks due to the simpler ion transport actuators [10].
Noncovalent interactions provide the gels with the very high mechanical strength and excellent self-healing ability of supramolecular materials [11,12]. Based on these studies, we used ZrO 2 to fabricate supramolecular nanocomposites with the electrochemical behavior of an ionic liquid and mechanical strength of an ionogel polymer.
Among numerical approaches, finite element analysis is one of the most effective and most common information extraction methods to evaluate and optimize robot designs. By using this method, analytical models can be greatly simplified, which greatly increases the computational efficiency. The largest drawback is that ignorance of the nonlinear and constitutive model and the simplification of the computational model often lead to coarse solutions [13,14].
Numerical simulations offer sufficient insight for each case during general soft robot design.
Lee et al. used the finite element method (FEM) to successfully predict the mechanical behavior of an ionic polymer-metal composite (IPMC) actuator [15]. Wang et al. established a model based on the FEM to determine the electromechanical bending behavior of photocurable ionogel actuators (PIA) [16]. He et al. used the common large-scale finite element analysis software ANSYS to simulate an ILG, which is based on the SOLID186 element and the nonlinear hyperelastic Mooney-Rivlin model [17].
Because the use of an ILG actuator requires deep understanding of the mechanical properties of the soft robot, it is necessary to establish a theoretical model to obtain its performance index. To reduce the experimental cost and time, we systematically analyzed the ILG via numerical simulations and verified the accuracy of the calculation, which contributed to the development of the ILG in the soft robot. The numerical simulation results matched the corresponding experiments, proving the validity of the model [18,19].

Fabrication of the Ionic Liquid Gel
In the experiments, the ILG was composed of 1-butyl-3methylimidazolium tetrafluoroborate (BMIMBF 4 ), hydroxyethyl methacrylate (HEMA), 2-diethoxyacetophenone (DEAP), and ZrO 2 , with masses of 900 mg, 68.6 mg, 1.4 mg, and 30 mg, respectively. The mixed solution was then placed into a magnetic stirrer to form a suspension. Following that, the sample was placed on an ML-3500C Maxima-type cold light source for ultra-high-intensity UV radiation curing to induce polymerization. The UV intensity is 90000uW / cm 2 (15 ″ / 380 mm distance) [20]. Figure 1(a) shows a comparison of the morphology of the ionic liquid before and after gel formation: the left figure shows the liquid state, while the right shows the solid state formed after polymerization. The schematic diagram of the principle behind gel formation is illustrated in Figures 1(b) and 1(c). Under UV light irradiation and the action of the DEAP catalyst, the polymer matrix was cross-linked into a porous network structure.

Material Nonlinearity and Parameters
The morphological analysis of the freeze-dried sample by scanning electron microscopy (SEM) showed that a porous microstructure was present throughout the ionogel. Distilled water was used to replace the internal ionic liquid in the ILG after freeze-drying treatment, and then, an S4800 Hitachi high-resolution field-emission scanning electron microscope was used to scan the sample. Figure 2 shows the spatial structure of the ionic liquid carrier HEMA for the typical 3D porous structure with a 5000x magnification, in which its matrix is cross-linked to generate a 3D support skeleton, offering good mechanical strength and self-repairing performance. Due to the porous 3D network in the ionogel, gelated BMIMBF 4 retains a relatively high ionic conductivity.

Hyperelastic Hypothesis.
Assuming that the ILG is an isotropic, incompressible, hyperelastic body, the following  assumptions can be made based on the theory of continuum mechanics to study its mechanical properties.
(1) The strain energy function W of a unit mass of material is an analytic function of the strain tensor of the natural state, termed the hyperelastic hypothesis. If the rate of change of W is equal to the power of the stress, then the material is a hyperelastic material.
The mechanical properties of a hyperelastic material are described by the strain energy density function W, which has many functions.
(3) The volume of the material before and after deformation can be assumed to be the same. λ 1 , λ 2 , and λ 3 are set as the x, y, and z directions of the main (extension) deformation rate, respectively, given by 1 where x, y, and d are the length, width, and thickness, respectively, and x 0 , y 0 , and d 0 are the corresponding initial values before deformation. Because the material is incompressible, its volume is the same before and after deformation, giving

Hyperelastic Stress.
Considering the mechanical performance requirements for potential applications, three kinds of deformation states should be examined: Figure 3(a) shows the uniaxial tensile stress, Figure 3(b) shows uniform prestretching in the X and Y directions, and Figure 3(c) shows the Maxwell stress increased along the thickness.
The physical properties are mainly expressed by the strain energy function, and each model is a special form of this function [21][22][23]. Once the form of the strain energy function W is determined, the Cauchy stress tensor P can be given by where I is the unit tensor, which is the left Gauss deformation tensor, and p is the hydrostatic pressure resulting from the assumption of incompressibility.
where B is the component of the Green strain tensor. The relationship between the invariants and the principal elongation is of a function of B.
Hence, the isotropic and incompressible deformation process of an ILG is given as According to (4) and (6), we can obtain where I 1 , I 2 , and I 3 are the relative changes in the length, surface area, and volume of the elastomer, respectively.

Mooney-Rivlin Model.
The Mooney-Rivlin model was chosen after comparing various hyperelastic constitutive models. The mechanical properties of ionic gel materials can be studied by using the Mooney-Rivlin formula, which is considered to be a nonlinear finite element of ionic gels in this study [24][25][26]. The strain energy function in the Mooney-Rivlin model equation is as follows: where c ij is a constant.
The Mooney-Rivlin model is the most widely used strain energy function in the finite element method. It assumes that the strain energy density is a first-order function of the principal strain constant. Under large deformations, the 50.0 휇m S4800 3.0kV 8.2mm × 700 SE(M) Figure 2: SEM image of a freeze-dried BMIMBF 4 -based gel after replacing the ionic liquid with water [17]. mechanical properties of the ILG, as an incompressible hyperelastic material, are described.
3.4. Selection of the Constitutive Model for the ILG. Yeoh noted that the practical value of the higher-order strain energy function is small because the reproducibility of the ILG material is not sufficient and does not allow accurate estimation of a large number of parameters [27].
Due to its simplicity and practicality, the Mooney-Rivlin model is widely used in finite element analysis. The firstorder Mooney-Rivlin model describes the mechanical properties of incompressible hyperelastic materials under large deformation. The higher-order Mooney-Rivlin model can obtain a good approximation for the solution of large strain.
The mechanical response of the hyperelastic material model is determined by the strain energy density function. The Mooney-Rivlin constant of the material must be accurately evaluated to obtain a reliable result from the hyperelastic analysis. In finite element analysis, the hyperelastic material is generally assumed to be a homogeneous isotropic material whose elastic modulus (Young's modulus) E, initial shear modulus G 0 , and Poisson's ratio v satisfy the following relationship [28,29].
From (8), we get The stresses in all directions are where ∂W/∂I 1 and ∂W/∂I 2 are the partial differentials in the strain energy function W for I 1 and I 2 , respectively, and σ 1 , σ 2 , and σ 3 are the stresses in the x, y, and z directions, respectively. When uniformly stretched in the X direction, λ 2 = λ 3 , and from (3), we get From (6) and (16), we obtain I 1 = λ 2 1 + 2 λ 1 , Because only axial tensile deformation is considered, the stress in the other two directions is zero.
From (14) or (15) and (16), we obtain Substituting (19) into (13), For (11), the partial differentials of the strain energy function W for I 1 and I 2 are given by ∂W ∂I 1 = c 10 + 2c 02 I 1 − 3 + c 11 I 2 − 3 , Therefore, For (23), the partial differentials of the strain energy function W for I 1 and I 2 are given by ∂W ∂I 1 = c 10 + 2c 20 Therefore, From (8), we get The stresses in all directions are When uniformly stretched in the X and Y directions, λ 1 = λ 2 , and from (3), we obtain From (6) and (31), we get Because tensile deformations are only considered in the X and Y directions, the stress in the thickness direction is zero.

∂W ∂I 2 34
Substituting (34) into (28) or (29), For (26), the partial differentials of the strain energy function W for I 1 and I 2 are given by where ς is the insulation constant, ς 0 is the vacuum dielectric constant (8.85 × 10 −12 F/m), V is the electric field intensity, and u is the voltage. (9), the strain energy equation in the five-parameter Mooney-Rivlin model is

Production of Tensile Test
Sample. The prepared solution was poured into a mold and polymerized into a gel under the irradiation of a UV lamp. The gel is white in color and is filled in the transparent glass mold as shown in Figure 4. Then, the sample was cut to the size shown in Figure 5, with a thickness of 3.4 mm. As shown in Figure 6, the experimental instrument was a UTM2502 electronic universal testing machine. The mechanical sensor on the testing machine can achieve a precision of 0.1 mN, the displacement sensor has a precision of 0.001 mm, and the stretching rate is 500 mm/min.

Experimental
Results. As seen from Figure 6, the ILG becomes longer and thinner as the load increases, which is consistent with the assumption that the material is incompressible. The tensile stress-strain curve of the ILG is shown in Figure 7, and the average tensile strength (Young's modulus) of the material obtained from the tensile stress-strain curve is 7.6 kPa.
In the ionogel, BMIMBF 4 exhibits a high level of hyperelastic toughness when the tensile deformation reaches 360%. The tensile test showed that the tensile properties of the gel increased with an increase in ZrO 2 content. The increase in ZrO 2 content should generate more cross-linking sites and higher conversion rates (shown in Figure 2), which contributes to the overall mechanical properties. As the content of  ZrO 2 increases, the tensile strength of the ILG increases, while the elongation rate decreases. Considering the above data, the optimum amount of ZrO 2 to provide a large tensile strain and tensile strength is 3 wt%.

Stress Calculation.
The relationship between the real stress σ i and engineering stress σ E is where σ i is the true stress and σ E is the engineering stress. The calculated stress given in Table 1 is the engineering stress. We next tested whether the above-derived engineering stress expression for the ILG is accurate and whether it can be used for the design of an ILG actuator or sensor. To verify the accuracy and practicability of the induced stress expression of the ILG, the deformation rate at each point in the uniaxial tensile test was substituted into the engineering stress expression to calculate the engineering stress. A comparison of the calculated uniaxial tensile stress with the experimental data is given in Table 1.
For the five-parameter Mooney-Rivlin model, the parameters c 10 , c 01 , c 20 , c 11 , and c 02 are as follows. As seen from Table 1 in the comparisons of the calculated tensile stress with the experimental data, the relative error of the five-parameter Mooney-Rivlin model is less than 0.51% and the relative error of the nine-parameter Mooney-Rivlin model is no more than 8.87%.   Figure 8 shows that the stress-strain curve calculated by the five-parameter model almost coincides with the experimental curve. The first half of the stress-strain curve calculated by the nine-parameter model almost coincides with the experimental curve, while in the second half, the error between the stress-strain curve calculated by the nineparameter model and the experimental curve becomes increasingly larger, but the error remains small.
The above analyses indicate that the simulated values are consistent with the experimental values, making our derivation a feasible and credible stress formula for the calculation of the ILG properties.
Due to our current laboratory conditions, only the uniaxial tensile test was performed. The next step is to improve the laboratory conditions in order to carry out the plane uniform tensile test and the Maxwell stress experiments. Further studies of the five-parameter and nine-parameter Mooney-Rivlin models will be carried out.

Conclusions
In this paper, an ILG is modeled by the hyperelastic nonlinear finite element model. The simulation results show that the Mooney-Rivlin model can well adapt to the constitutive relation of the material [30,31]. The main advantage of the ILG is that the stress-strain curve can be obtained by the performance parameters of the material in a relatively short time, which provides a theoretical basis for the optimal design of a soft robot.
The simulation results show that the average error between the calculated data and the experimental data is small, and the model has a good correlation with the experimental data. The model requires the input of the ILG material parameters. A standard uniaxial stretching method is used to obtain the desired ILG material parameters.
In the future, we will seek a generalized algorithm for identifying the ILG mechanical properties. Notably, all the results of this study show that there is a good correlation between the 3D theoretical assumptions and the experimental conditions, which proves that our method can be used to optimize the design of a soft robot. This work opens a new path to study the performance of ILG soft actuators, which will be the direction of future work.