Thermally controlled large deformation in temperature-sensitive hydrogels bilayers

ABSTRACT The present work investigates the thermally controlled deformation characteristics in temperature-sensitive hydrogels bilayers. The free energy density for temperature-sensitive hydrogels is modified, upon which the finite element model is developed and implemented through user-defined material subroutine UHYPER in the commercial software ABAQUS. The modified UHYPER implementation allows for more vividly depicting the continuous deformation in phase temperature region for temperature-sensitive hydrogels. Several thermally controlled cases of temperature-sensitive hydrogel including grippers, self-folding boxes, thermally driven origami are presented to illustrate a wide array of complex interesting applications or phenomena. Furthermore, we develop a simple model to theoretically calculate the bending angle of the temperature-sensitive hydrogel bilayers, which has been validated by the finite element simulation results. Our study can provide more insights for optimal design in thermally controlled hydrogels structures. Graphical abstract


Introduction
Gels are created when flexible long-chain cross-linked elastic polymer networks absorb solvents and swell. When water is absorbed as a solvent, the resulting gel is called a hydrogel. Hydrogels have good reversible deformability and the degree of swelling is influenced by various factors such as mechanical forces [1], pH [2], temperature [3], light [4], electric field [5,6], etc. The expansion for hydrogel is usually large and reversible [7,8], thus developing it in promising applications for actuators [9], sensors [10], drug delivery [11], and bioengineering [12,13].
Since Flory and Rehner [14,15] presented a statistical mechanical model of gels in 1943, there has been a great deal of research work so far on the deformation theory of hydrogels [16][17][18][19][20][21]. Huang et al. [22] summarize constitutive models of main types of hydrogels for various stimuli. Lei et al. [23] mainly introduce the theoretical and numerical researches on the hydrogel network models, aiming to bridge the gap between the micro-structure and the macroscopic mechanical responses for hydrogels. Among them, the free energy function for gels is described as the sum of the free energy due to the stretching network and the free energy due to the mixing of the fluid with the polymer, together with the free energy contributions due to factors for certain types of gels [24]. Based on the Flory and Rehner model, Cai and Suo [25] adopt a theoretical approach to describe the deformation of Poly (N-isopropyl acrylamide) (NIPA) temperature-sensitive hydrogels by defining the interaction parameters as a function of temperature and polymer volume fraction. Chester and Anand [26] develop a continuum-level theory to describe the coupled mechanical deformation, fluid permeation, and heat transfer of temperature-sensitive hydrogel. Ding et al. [27] formulate the state of mechanical and chemical equilibrium via a variation approach, upon which the finite element model is developed and implemented. In addition to theoretically and numerically exploring the deformation behavior of temperature sensitive hydrogel, lots efforts have been paid on the experiments, which can clearly visualize the swelling or shrinking deformation of temperature-sensitive gels during the phase transition period. Liu Group quantitatively investigate the peculiar behavior on the NIPA hydrogel by using traditional mechanic testing and dynamic mechanical analysis [28,29].
Meanwhile, when exposed to the stimuli, a hydrogel is commonly attached or bonded to another material which does not swell in the solvent. The mismatch in deformation of two materials can be utilized in application of soft sensing and actuators. Hu et al. [30] propose a double-layered strip of temperature-sensitive hydrogels to model a thermally controlled gel hand. Yoon et al. [31] use a thin and gradient cross-linked hinge structure made from two layers of temperature-sensitive hydrogels with different swelling ratios to make a self-folding microscale hydrogel device by connecting thicker and fully crosslinked panels. Fan et al. [32] design dual-gradient structures of hydrogel sheets that can exhibit ultrafast inverse snapping deformation. Lucantonio et al. [33] describe the swelling-driven curving of bilayer gel beams. Morimoto and Ashida [34] investigate the finite bending response to temperature of a bilayer gel under plane strain conditions. Recently, Shojaeifard et al. [35] develop a new semi-analytical method to predict the bending of hydrogels and investigated design factors such as semi-angle, bending curvature, and aspect ratio. Tang et al. [36] realize the shape morphing in a multi-material system for programming complex soft structures. Dortdivanlioglu et al. [37] tune swelling-induced interface crease instabilities at hydrogel bilayers by controlling the material properties. With a sloped thickness design, Huang et al. [38] theoretically and experimentally investigate the reconfigurable bi-material beam structures, which could be transformed into some other 2D self-scroll and 3D self-helical configurations.
Here, we conduct a theoretical and numerical study on thermally controlled large deformation in temperature-sensitive hydrogels bilayers, thus can illustrate various complex interesting applications or phenomena.
The work in this paper is developed based on the single-phase theory of temperaturesensitive hydrogels. We used the commercial finite element software ABAQUS and a user subroutine to simulate several cases of non-uniform large deformations of temperaturesensitive hydrogels, including: thermally controlled grippers, self-folding box made of temperature-sensitive hydrogel, and some folding phenomena caused by temperature variations. Afterward, we develop a simple model to theoretically calculate the bending angle of the temperature-sensitive gel structures, which has been validated by the finite element simulation results.

Equilibrium conditions and governing equations
We specify the reference state as the network condition of the hydrogel in the dry state and denote the coordinates of the network cells in the reference state by X K . Under different external loads, such as contact with the solvent or subject to geometrical constraints or mechanical forces, the network of the hydrogel deformed; and we define the deformation occurring in this reference state as a function of the reference coordinates, noted as x i (X). Hence, the state of deformation is characterized by the gradient of deformation, which is defined as: When the hydrogel is in the deformed state, the state of the hydrogel is then described by the distribution of solvent molecules C(X) and x i (X). We assume that the liquid solvent molecules are incompressible, denoting the volume of each solvent molecule by v and the chemical potential of the solvent molecules by μ. The small deformation of hydrogel is denoted as δx i (X), defining the body force as B i (X) and surface traction as T i (X). According to thermodynamic regulations, the change in the free energy of the gel should be equal to the sum of the work done by the external mechanical forces and the external solvent, i.e.
This equation holds for any arbitrary variables δx and δC in equilibrium state. We only consider states of equilibrium for temperature sensitive hydrogels, in which the chemical potential of solvent μ is uniform throughout the external solution and the gel. Using Legendre transform, the new free energy function Ŵ is derived by: Once the function W (F, C) is determined and differentiated, Eq. (2) is an algebraic equation. C can be expressed as a function of F and µ. Thus, we can define Ŵ as a function of the network deformation gradient and the chemical potential of the solvent molecule, Ŵ (F,μ). The primary objective of this transformation, as detailed elaborated in previous work on hydrogel [17,19,27], is to convert the equilibrium Eq. (3) into a hyperelastic solid equilibrium equation, We define the nominal stress jointly by the chemical potential and the deformation gradient, i.e.

Flory-rehner free energy function
As the previous theoretical work on hydrogel [19,25], we also adopt the Flory-Rehner free energy function.
In the above expression, N represents the number of polymeric chains per unit reference volume, k B T is the temperature in energy form, v represents the volume of each solvent molecule, F iK is the deformation gradient with respect to the reference state, and C is the nominal concentration of water molecules. For temperature-sensitive hydrogels, the Flory interaction parameter χ is a function of temperature and is shown below [39]: Here, The symbol ϕ here indicates the volume fraction of the polymer in the temperature-sensitive hydrogel, A 0 , A 1 , B 0 , and B 1 are material parameters obtained from experiments. In the subsequent finite element simulations, we will use the material parameters of the previously mentioned poly(Nisopropyl acrylamide) (PNIPAM) gel. According to data provided by Afroze et al. [40], the values of its relevant parameters are A 0 = −12.947, A 1 = 17.92, B 0 = 0.04496, and B 1 = −0.0569/K.

Constitutive equations
We assume that all molecules in the hydrogel are incompressible, i.e. that the volume of the hydrogel is equal to the sum of the volume of the network in the dry case and the volume of the pure liquid solvent [17].
Eq. (8) represents constraint equation of incompressibility. Applying Eqs. (3) and (6), we can obtain an alternative form of the free energy function for temperature-sensitive hydrogels, namely: Substituting the above Eq. (9) into the previously stress expression Eq. (5) gives: Here, I = F iK F iK and J = detF are the invariants in the deformation gradient. The volume per molecule is represented by v = 10 −28 m 3 , k B T = 4 x 10 −21 J at room temperature and k B T/ v = 4 × 10 7 Pa. Nv is a dimensionless measure of the crosslink density of the polymer in the drying network and its value measures the long chain density of the polymer. H is the transpose matrix of the inverse matrix of the deformation gradient matrix F.

Implementation and improvement of the finite element method
When using the previous UHYPER subroutines for temperature-sensitive hydrogels [19,27], there will be the snap-through instability and the multiple-solution problems in the calculation of the free energy. Generally, the above problems appear in the region of lower critical solubility temperature (LCST). If the numerical calculation we model involves such phase transformation temperature state, the errors occur in the simulation. The simulations have the tendency to terminate abruptly when the phase transition temperature reaches. The phase transition temperature becomes the dividing line. In order to avoid the above problems in the calculation, the logarithmic term of the free energy expression can be replaced by a polynomial expression [41], and the first three terms of the expansion by Taylor's series were chosen as follows: Actually, the term log 1 À 1 J À � originally satisfies the well-known Boltzmann relationship in statistical methods [42,43], which is adopted to represent the entropy change in mixing long polymer with the solvent. In order to avoid the instability of the calculation generated by more terms and ensure the accuracy of calculation, the first three terms were chosen to eliminate the multiple solutions. After introducing the modifications, the previously mentioned free energy function Eq. (9) can be written as: To have a better idea of the modified free energy, we have plotted the free energy of Eqs. (9) and (12) at phase transition temperature region as a function of volume ratio at three critical temperatures, respectively. LCST for PNIPAM gel is near T = 305 K, so we choose T = 304 K, T = 305 K, and T = 306 K. As shown in Figure 1, at T = 305 K, the original free energy function derived from Eq. (9) has two local minimum locations, corresponding to two stable states. However, the free energy function derived from Eq. (12) has only one local minimum at a specific volume ratio, corresponding to a stable swollen state, thus snap-through instability and the multiple-solution problems are eliminated. We rewrite the UHYPER subroutine based on the modified free energy function (list in the Appendix).
A free swelling cube model is proposed to verify the modified subroutine, which is shown in Figure 2. We constrain the swelling of three of the faces along their normal and allow the other three faces to swell freely. The parameters we chosen are Nv = 0.01, μ/ k B T = 0 and a starting temperature of 290 K.
In order to verify the modified FEM model, the numerical solutions calculated by the improved subroutine and the previous subroutine, analytical calculations and the experimental data are all shown in Figure 3.
As it can be observed from Figure 3, simulation calculation by the former subroutine is discontinuous and the discontinuity occurs when the temperature reaches around 305 K. We need to split the simulation into two stages, i.e. first from 290 K to a singularity temperature of about 305 K, and then from 305 to 315 K. With our subroutine, the phenomenon of multiple solutions in the simulation is avoided and we can directly achieve the continuous calculation through phase transition temperature region. It is effective to depict the volume-temperature behavior for temperature-sensitive hydrogels, which change significantly during the phase transition temperature region.

Numerical modeling
In this section, we present several finite element models incorporating temperaturesensitive hydrogel applications and simulate them using the commercial finite element software ABAQUS and its subroutines UHYPER. We use thermally controlled bending of bilayers consists of temperature-sensitive hydrogel and polydimethylsiloxane (PDMS) Figure 1. The free energy as function of volume ratio at various temperatures corresponding the analytical calculations using Eqs. (9) and (12). materials, to realize the grippers, self-folding boxes, origami driven by temperaturesensitive hydrogel.  The blue diamond squares and the yellow dots are corresponding the numerical prediction using UHYPER by Ding et al. [1] and the our modified UHYPER, respectively. The red dots represent experimental data from Oh et al. [2], and the light green curves represent the analytical calculations.

Hydrogel grippers
Soft grippers, for their working flexibility and endurance, are widely popular in smart manufacturing, medical, and bioengineering. By utilizing the reversible deformation of temperature-sensitive hydrogels, we have simulated a thermally controlled gel gripper device formed by a temperature-sensitive hydrogel material and two kinds of PDMS materials with different Young's moduli, which are shown in Figure 4.
As shown in the model in Figure 4, the green part is the temperature-sensitive hydrogel material, the red part and the white part are two PDMS materials with different Young's moduli, among which the white part is the harder. The Young's moduli for the two PDMS are settled as E h = 100 Nk B T 0 and E s = 0.2 Nk B T 0 , respectively, and Poisson's ratio ν is 0.33. We fixed the bottom surface at the center of the grippers and specified the parameters of the temperature-sensitive hydrogel as Nv = 0.1, μ/ k B T = 0, λ 0 = 2.87372. The initial temperature was set in 290 K and allowed to increase gradually.
In dividing the elements, we assigned the element type as the 8-node linear brick element (C3D8). As the temperature increases, the temperature-sensitive hydrogel blocks which are fixed to the model shrink, thus driving the deformation of the attached elastic PDMS material. The deformation process is shown in Figure 5.

Self-folding box with temperature-sensitive hydrogel
When temperature-sensitive hydrogels are applied to the joints of components, the shape of the device can be changed by controlling the temperature. Various structures made by hydrogel that can open and close automatically in response to the change of stimuli have been researched. Through our subroutine, we also display a tetrahedral self-folding box model driven by the change of temperature. As is shown in Figure 6, the finite element model is simulated using a combination of temperature-sensitive hydrogel material and PDMS materials In the self-folding box model, we use the same material parameters and definitions as in previous models, and their corresponding colors in the model are the same as those we defined before, respectively. The connecting domain between the faces consists of a temperature-sensitive hydrogel and the elastic material PDMS. The bottom surface is constrained in three translational degrees of freedom and the initial temperature applied at the temperature-sensitive hydrogel is 290 K.
The element type is defined as C3D8 (8-node linear brick element). As the temperature rises, the temperature-sensitive hydrogel at the joints shrinks gradually, driving the movement of the attached regions made of elastic PMDS, and the box folds up automatically.
The deformation patterns at different temperatures for previous subroutine and our subroutine are provided in Figure 7. It is obvious to see that the self-folding box model, computed using our subroutine eventually folds together.

Origami driven by temperature-sensitive hydrogel
We arrange the distribution of the temperature-sensitive hydrogel to achieve the process of origami, forming a variety of interesting shapes. As the temperature change, some interesting phenomena happen that flat structures fold up into a paper aeroplane model and a windmill model.

Paper aeroplane model
The finite element model is shown in Figure 8. We utilize the temperature-sensitive gel material in green, the PDMS material with a smaller Young's modulus in red and the PDMS material with a larger Young's modulus in white. The Young's moduli are consistent with the previous model, E h = 100 and E s = 0.2 N k B T 0 , respectively, and Poisson's ratio ν = 0.33. The relevant parameters for the temperature-sensitive hydrogel were also chosen in accordance with the previous models. We designed this model to be axisymmetric. We have imposed constraints of six degrees of freedom on two points of the model, both of  which marked in red circles. In order to obtain the desired folding structures, we have placed temperature-sensitive hydrogel material in different areas on both surfaces of the model as required. The initial temperature of the temperature-sensitive hydrogel is set to 290 K for all locations. As the bending angle of temperature-sensitive hydrogel depends on the variation of the temperature, the final temperature is settled as required.
The 8-node linear brick element (C3D8) is chosen as the element type for this model, which was divided into hexahedral elements. With the temperature rising, the hydrogel shrinks, driving the entire model to fold up, as shown in Figure 9.

Windmill model
In the windmill model, the gradual curling of the windmill blades is achieved by the shrinkage of the temperature-sensitive hydrogel during temperature changes. The Windmill model is shown in Figure 10.
Considering that the four blades of the windmill are identical, we have selected one piece of the blades for modeling calculations, i.e. the dashed part in Figure 10 (a). The detailed model information is shown in Figure 10 (b). The material in the green parts of the long strip is the temperature-sensitive hydrogel material, the red parts are the PDMS material with a smaller Young's modulus, and the remaining white parts are the PDMS material with a larger Young's modulus, and their values and the relevant parameters of the temperature-sensitive hydrogel are the same as those we have chosen earlier. In the model, we fix one of the edges of the triangular model, limiting the rigid body displacement of this edge in three degrees of freedom. Starting from a starting temperature of 290 K, the temperature changed differently for the temperature-sensitive hydrogels according to the extent of bending. In addition, the type of element is the 8-node linear brick (C3D8) hexahedral element.
We finally mirror the blade calculation to generate the entire windmill model, so that the deformation pattern for the entire windmill model visually. The deformation patterns at various stages are shown in Figure 11.

Theoretical prediction for thermally controlled bending
In the above numerical cases, the bending angle of structures depends deeply on the temperature. Motivated by the thermal bending of bimetallic strips proposed by Timoshenko [44], the hydrogel gripper is analogized to a bi-layered cantilever beam model, which is shown in Figure 12. As the temperature of the temperature-sensitive hydrogel increases and the volume shrinks, the deformation of the grippers and the internal loads applied are shown in Figure 13.
Note that in following descriptions, the letters labeled '1' represent the relevant variables for the elastic temperature-sensitive hydrogel material, while the letters labeled with the subscript '2' represent the relevant variables for the soft PDMS. P 1 and M 1 represent the axial tension and bending moment applied to the temperature-sensitive hydrogel cross-section, respectively, and P 2 and M 2 represent the axial tension and bending moment applied to the soft PDMS material cross-section, respectively. Since there are no external forces applied to the bilayer beam model, all forces and moments at any cross-section must be in equilibrium, so that we can obtain Furthermore, according to Judy et al [45], we use an equivalent flexural stiffness EI ð Þ eq to equivalently represent the moment-bending relationship between the two materials. The equivalent flexural stiffness EI ð Þ eq is written as In Eq.ð 15), the expressions for K, m, and n are Then we derive the relationship between the internal moment M and the radius of curvature ρ, i.e.
In the beam model, because the longitudinal unit elongation of the temperature-sensitive hydrogel and the PDMS material must be equal at the surfaces where the two materials are attached to each other, we are able to obtain the following relationship.
where ε T denotes the strain of the temperature-sensitive hydrogel with the variation of temperature. Through geometric compatibility, and λ(T) is a function that describe the relationship between swelling ratio of temperature-sensitive hydrogel and temperature T. Combining Eqs. (13), (14) and (20), it is obtained that According to Eq. (21), we can express the curvature κ of the temperature-sensitive hydrogel bilayers as Figure 13. Temperature-sensitive hydrogel grippers with deformation and internal loading.
where E denotes Young's modulus, h denotes thickness, and w denotes width. The bending angle θ, defined as the angle between the tangents of the two ends of the gel model. The angle of bending θ can be expressed as In our calculation, the Young's modulus of the elastic PDMS material is settled as E 2 = 0.2N k B T 0 . Furthermore, based on the illustration of Liu et al [27], the modulus of hydrogels varies during the deformation of hydrogels. Thus, the reduced incremental modulus of a temperature-sensitive hydrogel is derived. E ¼ @s We define the extension of the temperature-sensitive hydrogel along each direction as λ 0 and λ 1 , respectively, and Eq. (10) yields Nv We assume the initial state of temperature-sensitive hydrogel is in free-swelling, i.e. λ 1 = λ 0 , so E f = 2. Meanwhile, we set h 1 = h 2 = 0.5, w 1 = w 2 = 1, the parameters of the temperaturesensitive hydrogel to Nv = 0.1, μ/k B T = 0 and ramped up from the initial temperature of 290 to 340 K. In what follows, we would focus on the variation of the model turning angle θ with length L. Plots of the temperature-bending angle relationship are drawn for L = 7, L = 9, and L = 11. The temperature-bending angle relationship for temperature-sensitive hydrogel bilayers is plotted and compared with the finite element model results in Figure  14. The deformation pattern for various temperatures is shown in Figure 15.
It can be seen from Figure 14 that the theoretical solution and the finite element results are in good agreement. As for the phenomenon that the value of the bending radius of the temperature-sensitive grippers is greater than the value of π, the reason is that the temperature-sensitive hydrogel is deformed significantly as the temperature increases and the bilayers can roll up several laps during the process of crossing the phase change point [46]. The above theory aims at developing a reliable but simple procedure to specific bending angles for temperature-sensitive hydrogel bilayers.

Conclusion
The prevalent theory of hydrogels based on Flory-Rehner's free energy can well describe the large deformation of temperature-sensitive hydrogels. However, this approach comes with snap-through instability and the multiple-solution problems for finite element implementation. To overcome these problems, the logarithmic part in the energy equation is replaced by its multinomial expansion. Meanwhile, numerical accuracy by this approach is verified with experimental results, as well as analytical solutions.
Using the above developments, we have provided several numerical examples on thermally controlled bending in temperature-sensitive hydrogels bilayers such as grippers, self-folding box devices, and flat folding phenomena.
In addition, based on Timoshenko's theory, we derive a bending theory for temperature-sensitive hydrogel bilayer. By compared the theoretical solution with the finite element results for a temperature-sensitive hydrogel gripper, we have validated .We  hope that the provided subroutines and the proposed method will provide a useful reference for future numerical simulations and experimental work on temperature-sensitive hydrogels.

Disclosure statement
No potential conflict of interest was reported by the author(s).