1 Introduction

Compared with traditional rigid robotics historically built for manufacturing product lines, soft robots, an emerging class of machines made from soft materials, offer unprecedented benefits of safety, adaptability, and complex motion^[1-14]. Recent advances in materials and manufacturing technologies have boosted the fast development of soft robots. In particular, smart soft materials that change their shape and size upon being exposed to external stimuli offer a promising avenue to soft robots. Examples include dielectric elastomers (DEs)^[3], electroactive polymers (EAPs)^[4-6], shape memory alloys (SMAs)^[7], shape memory polymers (SMPs)^{[6, 8]}, hydrogels^[9-10], and liquid crystal elastomers (LCEs)^[11-13]. For instance, Wang et al.^[11] developed an adaptive soft robot based on a thermal-responsive carbon black doped LCEs (LCE-CB) bimorph to mimic the locomotion behaviors of an inchworm. In addition, such an artificial soft robotic inchworm is able to sense the environmental change and to react to stimulation correspondingly. The thermal-responsive LCE-CB deforms with the change of its temperature due to phase transition^[14-16]. As temperature is heated higher than the critical temperature of phase transition, the thermal-responsive LCEs shrink. The main components of the adaptive fully soft bodied robot^[11] presented in Fig. 1(a) include actuators of ultra-thin stretchable heaters with the open-mesh shape^[17-21] and thermally actuated artificial muscles based on a bimorph structure made up of LCE-CB and polyimide (PI, Kapton tape). The open-mesh shaped heater is electrically connected through an anisotropic conductive film (ACF) cable to supply power. Joule heat generates with the current injection into the ultra-thin deformable heater. The phase transition of the LCE-CB is then induced as the temperature increases, which causes shrinkage of the LCE-CB layer and thermal expansion of the PI layer in the bimorph structure, respectively. Since the thermal expansion of PI is pretty small (≈ 0.01%) while the shrinking deformation of LCE-CB is relatively large (≈ 24%) when the temperature reaches 100 ℃ after phase transition^{[11, 22]}, a bending deformation of the bimorph structure shown in Fig. 1(b) will then be induced by the interface strain mismatch, which recovers its original shape after cooling.

The thermo-mechanical behavior of soft robots based on the LCE bimorph structure is important to their operation. Experimental characterization of the thermo-mechanical response time to reach the steady deformation is about 30 s^[11]. To better understand the thermo-mechanical responses and to achieve improved designs for soft robots, an analytical model for describing the transient thermo-mechanical behavior is critically necessary to be established. To this end, this paper focuses on developing a transient thermo-mechanical analytical model to precisely capture both thermal and mechanical responses of the bimorph structure, which is further validated with the finite element analysis (FEA). The predictions of the temperature distribution and the curvature radius response induced by the heat generation from ultra-thin stretchable heaters are obtained from the model.

2 Mathematical modeling

Due to the two reasons as follows: (ⅰ) the thickness of the LCE-CB/PI bimorph structure (500 μm) larger than that of the stretchable heater composed of PI and gold (2.8 μm) by 2 orders of magnitude, and (ⅱ) the significantly reduced bending stiffness from the open-mesh structure of the heater^[11], it is rational to model the heater as a planar heat source which has few impacts on the bending of the bimorph structure. Considering that the stretchable heater's in-plane dimension (16.6 mm×7.9 mm) is almost the same as that of the bimorph structure (17.0 mm×9.4 mm)^[11], the heat transfer of the structure can be approximately treated as a one-dimensional problem. As shown in Fig. 1(c), the soft robot^[11] is then simplified as a bimorph structure with air convection boundary conditions on both sides (T_air and h denote the ambient temperature and the natural convection coefficient, respectively). The origin of the axis z along the thickness direction is located at the bottom surface of PI. The thickness of the PI layer is z₁ while the thickness of the LCE-CB layer is (z₂-z₁). Then, the Fourier heat conduction equations for this one-dimensional problem can be given as

(1)

where ΔT=T-T_air defines the temperature increase, which is the difference between temperature and T_air, and α and k are the thermal diffusivity and thermal conductivity of PI and LCE-CB with the subscripts 1 and 2, respectively. Q(z, t)=Q₀Q₁(z)Q₂(t) is the time-dependent planar heat source^[23], where Q₀ is the amplitude of the input power, Q₂(t) represents the time-variation of the unit power, while the δ function is used to describe Q₁(z) as

where .

The natural convection boundary conditions at both top (z=z₂) and bottom surfaces (z=0) of the bimorph structure and the temperature/heat flux continuous conditions at the LCE-CB/PI interface (z=z₁) can be then given as follows:

(2)

The initial conditions are

(3)

As a general and effective method, Green's function method^[24] is used to solve Eq. (1). In order to obtain the appropriate Green's functions for the nonhomogeneous partial differential equations in Eq. (1), the homogeneous version of Eq. (1) is considered for the same region, which gives

(4)

and the corresponding boundary/initial conditions can be written as

(5)

(6)

Then, the method of separation of variables θ_i(z, t)=Z_i(z)Γ(t), i=1, 2 is used, and Eq. (4) turns to be

(7)

with the following solutions:

(8)

where

(9)

and the eigenvalues λ_1n and λ_2n satisfy

(10)

with

(11a)

(11b)

Combining the initial conditions given in Eq. (6), the complete solutions to Eq. (4) are obtained by the superposition principle as follows (z₀=0):

(12)

with and , which can be written in terms of Green's function as

(13)

and the forms of Green's functions of this problem are obtained by comparing the solutions in Eqs. (12) and (13) as

(14)

The desired Green's function can be derived by replacing t in Eq. (14) into (t-τ) as

(15)

Then, the solutions to Eq. (1) can be obtained according to Green's function method as

(16)

where z₀=0, and g_j(z', τ) denotes the heat source term in the jth layer. Combining initial conditions and the heat source term together with Eq. (16), we can obtain

(17)

Then, the analytical temperature predictions are presented by substituting Eq. (15) into Eq. (17) as follows:

(18)

where

(19)

As the temperature increases with time, the phase transition of LCEs will lead to the gradual shrinkage of the LCE-CB layer. At the same time, the thermal expansion of the PI layer will occur. Then, the bimorph structure will start to bend, i.e., induce the locomotion of the soft robot. Such kind of bending can be described by a bilayer beam model with initial strains in both layers^[25], which are provided by thermal expansion of the PI layer and shrinking of the LCE-CB layer, respectively. Considering that the shrinking deformation of the LCE-CB layer (≈ 24%) is much larger than the thermal expansion of the PI layer (≈ 0.01%), the bending from mechanical mismatch is mainly attributed to the shrinkage of the LCE-CB layer^[11]. Thus, only the initial thermal strain in the LCE-CB layer is taken into consideration in the following analysis. The bending normal stress σ_i of both layers can be obtained from Hooke's law,

(20)

where c denotes the uniform strain, z_b represents the neutral axis location, R is the bimorph structure' curvature radius, and E_i is Young's modulus of the ith layer. The shrinking strain ε₂⁰(T₂) in the LCE-CB layer can be obtained from the fitting functions^[26] of the experimental data of shrinking strain's dependence on the temperature^[11].

Since Young's modulus of the PI (≈ 2 GPa) is 4 orders of magnitude larger than that of LCE-CB (≈ 0.5 MPa), c is small enough to be negligible. An approximation based on the equilibrium of the bimorph structure can be given as follows:

(21)

3 Results and discussion

The validation of the analytical results is carried out with the FEA based on the continuum element C3D8RT of the commercial software ABAQUS without any assumption. The stretchable heater is a PI/Au/PI sandwich structure with the thicknesses of 1.3 μm/200 nm/1.3 μm^[11], respectively. The thicknesses of the PI and LCE-CB are 0.05 mm and 0.5 mm, respectively^[11]. The thermal conductivity, thermal diffusivity, Young's modulus, and Poisson's ratio are, respectively, 317 W·m^-1·K^-1, 1.54 ×10^-6 m²·s^-1, 79 GPa, and 0.44 for gold^[27]; 0.4 W·m^-1·K^-1, 1.27 ×10^-4 m²·s^-1, 0.5 MPa, and 0.5 for LCE-CB^{[11, 28]}; 0.12 W·m^-1·K^-1, 7.75×10^-8 m²·s^-1, 2.5 GPa, and 0.34 for PI^[29]. The shrinking deformation dependence on the temperature of the LCE-CB from experiments is fitted with a piecewise linear function^[26]. The thermal expansion coefficients for gold and PI are 15.4 ×10^-6 K^-1 and 17.0 ×10^-6 K^-1, respectively^[26]. The natural convection coefficient is set as h=12 W· m^-2·K^-1[26] while T_air = 25℃, and the input power from heat source is 237 mW^[11].

Figure 2(a) illustrates the condition of continuous constant heating input with Q₂(t) = 1. The maximum temperature in the bimorph structure occurs at the interface of the LCE-CB and PI (z=z₁), and the time-variation of the maximum temperature derived from Eq. (18) is shown in Fig. 2(b), which agrees well with the FEA results. The maximum temperature increases from the initial ambient temperature of 25℃ to the steady-state temperature of 86.8℃. It is shown that the temperature takes almost 30 s to reach the steady state, which indicates that the continuous constant heating source can affect the bending behavior for a relatively long time. The analytical predictions on the curvature radius from Eq. (21) are presented in Fig. 2(c). The results from the FEA agree well with the analytical predictions. The asymptote with the value of 3.921 mm represents the curvature radius under the steady-state temperature, which agrees well with the reported experimental result^{[11, 26]}. The curvature radius drops from 102.8 mm at 4 s to 4.7 mm after 20 s of heating. With the increase in the heating time, the changing trend of the curvature radius tends to be gentle, and the curvature radius gets closer and closer to the asymptote. Considering the situation that the motion of the soft robot needs to be controlled, a pulsed operation of the heater can be used to solve this problem. As shown in Fig. 2(d), the pulsed injected power can be modeled as

(22)

Then, the time-variation of the maximum temperature with the first three periods is presented in Fig. 2(e) with t₀=25 s and t_p=25 s. The temperature is almost stable after 3 periods as well as indicated by the curvature radius shown in Fig. 2(f). Through such a type of pulsed heating operation, the bending behavior of the soft robot can be controlled in a certain range (4.2 mm-27.4 mm) as shown in Fig. 2(f), instead of constant heating and complete cooling, which is an applicable method to control the soft robot.

Figure 3(a) presents the LCE-CB layer thickness influence on the bending curvature radius of the bimorph structure. As the LCE-CB thickness increases from 0.5 mm to 1.5 mm, the bimorph structure's bending curvature radius increases sharply at a certain time. Figure 3(b) shows the changing trend of the curvature radius versus the heating time under different modulus ratios E₁/E₂ of the bimorph structure. The results show that a larger modulus ratio will lead to a larger curvature radius at a certain time. These results indicate that, a thinner LCE-CB layer, the smaller modulus ratio E₁/E₂, and longer heating time before the steady state are helpful to generate a smaller curvature radius, i.e., significant movement of the soft robot.

Another important parameter of concern is the response time, which is defined as the heating time when the curvature of the structure reaches 80% of the curvature at the steady state in this paper. The response time can be used to predict and design the motion period of the soft robot. The effects of the thicknesses of the LCE-CB and PI on the response time are shown in Figs. 4(a) and 4(b), respectively. The results suggest that a thicker structure will increase the response time under the same input power of 237 mW, which is acceptable. Since the bending of a thicker structure needs more energy, the response time will increase due to the constant input power. Figure 4(c) illustrates the influence of the input power on the response time. The response time tends to first increase and then decrease with the increase in the input power, which is mainly indicated by the strain-temperature curve of LCE-CB^[11].

4 Conclusions

In summary, we have established a transient thermo-mechanical analytical model to predict the temperature field and locomotion of the LCE-CB bimorph based soft robot. This model precisely captures the thermal-induced bending process of the soft robot due to the Joule heating from the ultra-thin stretchable heater. Both analytical temperature and bending curvature predictions are in good agreement with those from the FEA. The key design parameters of the soft robot including the LCE-CB layer's thickness and the modulus ratio which affect the thermal response of the bimorph structure are systematically investigated. It is concluded that a thinner LCE-CB layer and a lower modulus ratio could yield more bending deformations. The response time and pulsed joule heat loading operation are also studied to shed light on the control and design of the soft robot. Overall, these results, from an analytical perspective, provide important insights and guidance for the design, fabrication, and control of the thermal responsive LCE based bimorph structured soft robots.