1 Introduction

Maintained in an aqueous environment, hydrogels are water-swollen networks of hydrophilic polymers, which are capable of imbibing a vast amount of water in response to an alteration of the external stimuli. The applications of various environmental stimuli, e.g., pH^[1-7], temperature^[8-10], light^[11-12], and electric and magnetic fields^[13-14], affect the amount of the absorbed water and the swelling ratio of hydrogel structures. The special properties, e.g., the responsiveness to the stimuli and the reversibility, play a vital role in the applications of smart hydrogels, e.g., responsive sensors and actuators^[15-16], drug delivery^[5], and tissue engineering^[17]. Therefore, hydrogels have drawn great attention in recent years^{[5, 18-20]}.

The mechanism of solvent uptake stems from the stimuli-sensitive properties and the composition of the polymer chains of hydrogels. In most of the constitutive models, the Flory statistical-mechanical model^[21] is used to additively decompose the free energy density into mixing and stretching parts with the consideration of thermodynamic equilibrium. Due to the hydrogel sensitivity, the additional energy density parts are added to the former phenomena, e.g., pH^[1-2], temperature^{[9, 22]}, light^[12], sensitivities, and neutral gels^[23-24]. Among the aforementioned stimulus-responsive hydrogels, there are some acidic groups attached to the polymeric-chains in the pH-sensitive hydrogels ionized gradually through the pH variation of the environmental solvent. The vast applications of the special behavior of pH-responsive hydrogels promote researchers to investigate within the theory framework. Particularly, Marcombe et al.^[1] proposed a theory for the pH-sensitive behavior, where the free energy density was an additional contribution to network stretching, component mixing, and ionizable group dissociating. Some cases have been set up to demonstrate the accuracy of this model^[25-27].

Among the vast applications, pH-responsive sensors and actuators are of significant importance for biological, environmental, and medical applications. The application of bending structures is a regular strategy in the design of sensors^[28-30]. The development of responsive bending structures stems from the applications of inhomogeneous external stimuli or the utilization of layered structures, and the implementation of multilayered construction is a more conventional strategy in the design of sensors. Smart multilayered structures usually bend, since the swelling ratios of different layers are different^[15]. Timoshenko^[31] was a pioneer of performing the research on multilayered structures. He probed the bi-metal sensor. Hu et al.^[30] fabricated another generation of bilayer sensors. They constructed the bilayer through the interpenetration of an N-isopropyl acrylamide (PNIPA) layer with a polyacrylamide (PAAM) layer in some parts. The bilayer deformed in a circular shape in response to the temperature change. Furthermore, similar bending behaviors contained elastomers and hydrogels were used as sensors or actuators, e.g., bilayers and trilayers^{[28, 32-33]}, where the semi-analytical method was used to investigate the reversible finite bending of the responsive multilayer, and the finite element (FE) analysis was performed to validate the semi-analytical results.

Recently, a series of laboratory experiments have been conducted by Guvendiren et al.^[34-36] on the depth-wise structure with cross-link gradients, which was composed of hydrogels, to study the instability patterns of the surface. Polyhydroxyethyl methacrylate (PHEMA) hydrogels were fabricated with an ultraviolet (UV)-curable precursor solution composed of polymerized PHEMA, crosslinker, and a photo-initiator during the diffusion of oxygen. They observed that the critical condition of the surface instabilities depended on the length scale of the surface patterns on the cross-link gradient. Furthermore, they expressed that the cross-link density gradient along the film thickness direction was controllable by considering the factors such as the initiator and concentration of the crosslinker, the intensity of the UV exposure, the time, and the thickness of the strip. Later, using various methods, Wu et al.^[37-39] conducted some stability analyses of materials with cross-link gradients along the thickness direction. With the state space method, they predicted the onset of the surface instability of the structure in the linear elasticity for the small-strain assumption^[39] and hydrogels^[38]. Similar to smart multilayer structures, a functionally graded (FG) hydrogel has the capability of regenerating bending in response to the exterior stimuli. With the FG hydrogel strip, of which the material properties altered along the thickness direction, the interfacial adhesive forces vary continuously and makes the structure suitable enough to be used as sensors and actuators.

In order to figure out the underlying swelling-induced large deformation of hydrogels, recently, some analytical methods have been developed for multilayer responsive smart materials. Regarding the advantages of analytical methods and the applications of smart materials, the behavior prediction of hydrogels has drawn a great deal of interest. Roccabianca et al.^[40] presented the finite bending solution for thick multilayer elastomers under the plane strain condition. The analytical method for incompressible Neo-Hookean elastomer was evaluated through the experimental data, and the existence of neutral axes was studied, especially the geometries. Lucantonio et al.^[41] used analytical and numerical methods to describe the swelling-induced deformation of beams composed of neutral hydrogels. Morimoto and Ashida^[10] analytically studied the bending of a bilayer made of temperature-responsive hydrogels under the plane strain condition. They predicted that the bending in the finite model was larger than that in the linear model. However, the results were not validated. Later, Abdolahi et al.^[32] developed a semi-analytical method for a temperature-responsive hydrogel, and validated the obtained results with the results from the finite element method (FEM). They showed that the assumption of an intermediary virtual state over-constrained the deformations.

In this paper, a semi-analytical approach based on the swelling theory presented by Marcombe et al.^[1] is presented to solve the swelling induced finite bending of FG pH-sensitive hydrogel strips. The hydrogel strip assumed in this article is a single smart layer, whose the cross-link density distribution is altered along the thickness direction. A Cartesian coordinate system is defined at the initial reference. The deformation gradient tensor is introduced without considering any intermediary virtual state. According to the variation of the pH applied to the aqueous environment, the FG pH-responsive layer is deformed into a circular shape. Thus, the final configuration is defined in the polar coordinate system. The obtained results of the semi-analytical approach are in excellent agreement with the FEM outcomes. Some applied case studies are set up to investigate the robustness of the semi-analytical approach and to clarify the various parameter effects on the bending behavior. In contrary to the discontinuity in the deformation and stress fields of the multilayer structures, the FG layer regenerates proper finite bending with continuous trends.

The paper is arranged as follows. In Section 2, the finite bending problem and the geometry are introduced. After dividing this section into two subsections, the deformation kinematics and theory of the swelling induced finite bending of a pH-sensitive FG hydrogel layer are fully discussed. Later, considering the balance equation under the plane strain condition, the stress relations are introduced. Then, in Section 3, the FEM is used to illustrate the accuracy of the presented semi-analytical approach. Finally, in Section 4, we present the conclusions with a brief summary.

2 Finite bending of a pH-sensitive FG hydrogel strip

As illustrated in Fig. 1, a single pH-sensitive layer is assumed, whose cross-link density distribution varies along the thickness direction. Under the free stress state, in the initial configuration, the FG strip is exposed to a constant pH with the rectangular shape. Since the hydrogel strip is immersed in a bath with higher pH, and the smart FG strip begins to bend regarding various cross-link density distributions. The non-uniform cross-link density along the thickness makes the hydrogel layer deform dissimilarly in all directions due to the various swelling ratios in different parts of the strip. As depicted in Fig. 1, the upper part of an FG layer bears a lower cross-link density. The cross-link density can be decreased with various approaches^[37]. Considering the plane strain condition, the out-of-plane dimension is assumed to be infinite.

2.1 Kinematics

Regarding the stress-free condition in the initial state, the Cartesian coordinate system framework is introduced to illustrate the dry hydrogel under the plane strain condition. Referred to Fig. 1, e₁ and e₂ are assumed to be the in-plane axes, and e₃ is defined as the out-plane axis. The position vector in the initial state configuration is given as follows:

(1)

where

and L and H are the length and thickness of the hydrogel layer in the undeformed state, respectively. The total deformation gradient maps the initial state into the current configuration with a semi-annular shape. This process occurs by controlling the solvent pH from the initial one to the desirable one. The variation of pH makes the pH-responsive hydrogel swell. However, the distribution of the cross-link density N=N(X₂) is the main reason for the bending in the semi-annular shape, which is a desirable feature that can be utilized in various applications such as sensors and actuators. The deformed configuration is introduced in the polar coordinate system as follows:

(2)

where

θ is the semi-angle, and r₁ and r₂ are the inner and outer radii of the FG hydrogel layer, respectively. Assume

Then, the total deformation gradient tensor correspondent to the hydrogel layer can be denoted by

(3)

where ⊗ stands for the dyadic product operator, and r(X₂) is the current radius of the FG layer at X₂. Considering the linear form for the tangential component of the polar coordinate position vector with respect to X₁, we can calculate the principal stretches in the FG hydrogel strip as follows:

(4)

2.2 Free energy density function and stress distribution

In this section, we will use the constitutive model developed by Marcombe et al.^[1] to study the swelling behavior of pH-sensitive hydrogels. The Helmholtz free energy density of pH responsive hydrogels is assumed to be a function of the nominal concentrations of the ions and the deformation gradient. The free energy density are decomposed additively into four parts as follows:

(5)

where C denotes the nominal concentration of the particles, and the subscripts H⁺, +, and - stand for the hydrogen ion, the counter ion, and the co-ion, respectively. In addition, W_net, W_sol, W_ion, and W_dis are parts of the free energy density, which are the results of the network stretching, the mixing of the solvent with the network, the mixing of mobile ions with the solvent, and the acidic dissociating groups, respectively. The free energy of the network stretching and the mixing of the solvent with the network are defined as follows:

(6)

where G is the elastic modulus of the hydrogel strip defined by

in which N, K, and T are the density of the polymer chains, the Boltzmann constant, and the absolute temperature, respectively. I and J are the first invariant of the right Cauchy-Green tensor and the determinant of the total deformation gradient tensor, respectively. The parameter v represents the solvent molecule volume, and χ introduces the dimensionless mixing parameter. The contribution due to the mixing of the mobile ions and solvent is defined as follows:

(7)

where the subscript m denotes the mobile ions, i.e., the hydrogen ions, the counter ions, and the co-ions. c_m^ref is the reference value for the concentration of mobile ions. The contribution of the acidic dissociating groups is

(8)

where the subscripts AH and A are for the associated acidic groups and fixed charges, respectively. The parameter γ is the value of the enthalpy energy increase during the dissociation of an acidic group. The following relations among the particles are confirmed due to the prevailing electro-neutrality in the gel and exterior solution and the consistency of the numbers of the acidic groups:

(9)

where f is the number of the existing acidic groups attached to the polymer networks divided by the monomer number. To find out the principal stresses for the FG hydrogel layer, Eqs. (6), (7), and (8) are substituted into Eq. (5) to recast the total free energy density. The Cauchy stress and total free energy density of a hydrogel strip have the relation as follows:

Utilizing the total free energy density, the principal stresses for the pH-sensitive hydrogel strip are recast as follows:

(10)

where c_m denotes the nominal concentration of species in the current state, and it is defined by c_m=C_m/J.

In order to simplify the total free energy density and principal stresses, the above relations are required to be rewritten as a function of just C_H⁺ and C₊. Thus, the ionic equilibrium condition, known as the Donnan equations, is required to be considered as follows:

(11)

and the equilibrium of the chemical dissociation should be

(12)

where N_A=6.023×10²³ is the Avogadro number, and K_a represents the constant of the acidic dissociation whose dimension is similar to the particle concentration. Recasting Eqs. (11) and (12), we receive at a cubic equation for C_H⁺ as follows:

(13)

where

(14)

We determine the principal stresses as a function of the radius, its derivative with respect to X₂, and the semi-angle, and substitute the tangential and radial stretch equation (4) into the radial and tangential stress equation (10). Then, we have

(15)

(16)

where

(17)

and c_H⁺ is calculated by solving the cubic equation (13) as a function of r(X₂), r'(X₂), and θ. For the sake of briefness, the derivatives of variables with respect to X₂ are indicated with a superscript '. In order to determine the variables, in light of symmetry, the equilibrium equation is defined as follows:

(18)

Substituting the radial and tangential stresses (see Eqs. (15) and (16)) in the equilibrium equation (18), we have

(19)

where

(20)

Summarizing the above equations, c_H⁺ is a function of r(X₂) and r'(X₂), and θ is not substituted into Eqs. (19) and (20) but is considered in this approach. Considering the FG hydrogel, the derivative of the cross-link density is observed with respect to the X₂ in Eq. (20).

To solve the nonlinear second-order differential equation (18), the Donnan equations and the equilibrium of the chemical dissociation must be applied to the stress relations to recast the equations as a function of the hydrogen concentration. It is noted that, the hydrogen concentration depends on J. Next, the cubic equation is required to be solved to find the hydrogen concentration as a function of the radial and tangential stretches. The cubic solution is substituted into the principal stress equations (15) and (16). The stress equations are provided to be substituted into the equilibrium equation. Employing the three stage Lobatto IIIA formula^[42] results in the nonlinear second-order differential equation (18), which needs boundary conditions to be solved. The upper and bottom surfaces of an FG hydrogel layer are assumed to be traction free, i.e.,

(21)

This approach is given to be applied through an iterative technique. However, in the first iteration, the semi-angle is assumed by an initial guess. After solving the first iteration, the net tangential force and the net tangential moment components are considered to be zero for any strip cross section, i.e.,

(22)

(23)

Considering the above equations, the initial guess for the semi-angle in any iteration can be calculated with a optimization technique^[43] through minimizing the residual tolerance of the equation system. The net tangential forces and moment are calculated by utilizing the trapezoidal numerical integration. The iterative method proceeds until the residual value becomes desirable. Then, the final value of the semi-angle is used to calculate r(X₂) and r'(X₂).

3 Results and discussion

In this section, to illustrate the robustness and accuracy of the proposed solution method, some problems are analyzed analytically and numerically for the finite bending of an FG pH-sensitive hydrogel. A rectangular FG pH-sensitive hydrogel layer is assumed with the thickness of 0.01 m and the length of 0.02 m under the plane strain condition. The material parameters are set as follows^[1]:

The cross-link density of the FG strip varies along the thickness direction with an exponential form similar to the approach of Almasi et al.^[44], i.e.,

(24)

where the subscripts B and T indicate the properties of the bottom and top faces of the hydrogel layer, respectively. The linear cross-density distribution is applied through setting k=0, and the cross-link densities of the top and bottom faces are 0.001 and 0.01^[37], respectively. The pH value of the reference state is 2.

In order to analyze the same problem in the FEM, a rectangular hydrogel layer is divided into several strips with constant properties, i.e., 20 strips, 40 strips, and 80 strips. Due to the symmetry assumption, one half of the FG hydrogel layer is modeled with proper boundary conditions. The 2D plane strain deformable part is considered with a 4-node bilinear plane strain quadratic element type mesh, which is named as CPE4. Furthermore, the swelling behavior of the pH-sensitive hydrogel developed by Marcombe et al.^[1] is implemented as a user-defined material subroutine UHYPER of the commercial FE software ABAQUS. The UHYPER subroutine is used to study the large deformation of the hyperplastic materials, and the free energy density and its derivatives with respect to the first and third invariants of the deformation gradient are required to define the subroutine. Neglecting the end effect, an extra length is considered in the FE solution, which is ignored in the reporting results. The analytical and numerical results are compared to clarify the capability of the presented semi-analytical scheme.

As illustrated in Fig. 2(a), the particular values of the exponential coefficient k are assumed to define the cross-link density distribution along the thickness direction using Eq. (24), and k=0, ±5, and ±10 indicate the linear and exponential distributions similar to Ref. [38]. The free swelling ratio of the FG pH-sensitive hydrogel layer is shown in Fig. 2(b).

In order to evaluate the presented semi-analytical approach, the results of the swelling induced finite bending of an FG pH-sensitive hydrogel layer obtained with the semi-analytical method and FEM are compared. Regarding the distribution of the cross-link density, the distributions of the current radius and the radial and tangential stresses are investigated. In order to demonstrate the convergence of the results, the cross-link density of the following analysis is assumed in the linear pattern (k=0), and the FE analysis is conducted with a rectangular hydrogel layer for 20 strips, 40 strips, and 80 strips. The bending behavior of the hydrogel layer is studied after applying pH variations to the external solvent. The results of both the methods are presented at pH=5, 6, and 8.

As illustrated in Fig. 3, the distribution of the current radius along the thickness direction is demonstrated with 20 strips, 40 strips, and 80 strips in the FE analysis at pH = 8, which agrees well with those of the proposed semi-analytical approach.

The radial and tangential stresses are calculated under the same condition. From Fig. 4, it can be seen that, when the number of hydrogel divisions in the FE analysis increases, the outcomes of the radial and tangential stress fields will converge to the semi-analytical results. The FE results are accurately the same as the semi-analytical approach when the number of strips is 80. Thus, the following FE results are reported for 80 strips in a hydrogel layer. As depicted in Fig. 4, the radial stresses at the inner and outer radii vanish, which fulfills the boundary condition in both the semi-analytical method and the FEM. In addition, the FE tangential stress is chattering through the thickness direction, which stems from the divisions of the pH sensitive hydrogel layer with different material properties.

To complete the verification, as shown in Fig. 5 the above case study is carried out at various pH including 3.5, 4.5, 5.5, and 8, and the stress field results of the FEM and the semi-analytical method are compared. As mentioned in Refs. [1] and [28], most of the swelling of finite bending occurs.

As shown in Fig. 6, the normalized radial and tangential stress contours in both the semi-analytical method and the FEM are presented in detail at pH = 8 on the deformed configuration.

The distributions of the stress fields are perfectly identical in both the semi-analytical method and the FEM. It is observed that the stress values do not vary along the tangential direction. Therefore, the aspect ratio has no influence on the stresses as mentioned in Refs. [29] and [32]. The stress fields of the FG hydrogel layer in the semi-analytical method agree with those of the FEM, which shows the accuracy and capability of the proposed method.

As illustrated in Fig. 6, the bending direction of the FG pH-sensitive hydrogel layer represents that the swelling ratio of the top strip is higher than that of the bottom strip. In addition, the continuous stress fields are observed in the FG hydrogel layer. However, the multilayer hydrogels undergo discontinuous ones^{[29, 32]}.

After validating the results of the semi-analytical method, the current radius distributions of the FG pH-sensitive hydrogel with different k are depicted in Fig. 7. It is shown that the patterns of the current radius variations are similar for different k.

The effects of the cross-link density on the radial and tangential stress fields are investigated at pH = 8, as shown in Fig. 8, where the radial stress is ignored in the small-strain Timoshenko bilayer beam theory^[31].Having inhomogeneous stress fields along the thickness direction, the swelling ratio is inhomogeneous due to the direct effects of the mechanical stresses on the swelling ratio. From Fig. 8, we can see that the tangential stress is about two orders of magnitude higher than the radial stress, which causes the tangential stress to be dominant. Therefore, the free stress location is a decisive factor for the design of smart stretchable electronics with such materials. As shown in Fig. 8(b), there are two points around 0.2 and 0.8, of which the tangential stress vanishes through finite bending.

When the external solvent pH increases from 2 to 8, the FG pH-responsive hydrogel layer will swell in a circular shape due to the cross-link density distribution. The deformed configuration of the FG hydrogel layer is introduced by the semi-angle θ and the bending curvature 1/r of the FG layer, which are considered as the vital design parameters of hydrogel sensors and actuators^{[29, 32]}. The bending curvature defined for the inner radius normalized by the FG hydrogel initial length and the semi-angle is compared with respect to the cross-link distribution, as shown in Fig. 9. It is shown that, when the the hydrogel stiffness decreases, the semi-angle increases while the curvature bends. As demonstrated in Fig. 9, the hydrogel layer with the exponential coefficient k=-5 reaches the maximum semi-angle and bending curvature. As reported in Refs. [1] and [28], the major part of the swelling of a pH-sensitive hydrogel takes place between pH = 3.5 and pH = 5.5. In addition, before pH = 3 and after pH = 6, the swelling ratio does not vary significantly. It is noteworthy that the acidic groups will be forced to be associated when pK_a is greater than pH while dissociated when pH ≫ pK_a.

Due to the importance of the aspect ratio in the design procedure of sensors and actuators, a parametric study is conducted. As discussed earlier, the radial and tangential stresses do not change along the tangential direction (see Fig. 6), which causes the aspect ratio to have no effect on the stress fields. In addition, when the aspect ratio changes, the bending curvature is constant while the semi-angle and the bending of the hydrogel layer both change. According to Fig. 10, the semi-angle has a linear relation with the aspect ratio, while the effect of the aspect ratio on the bending of the hydrogel layer is obvious, which shows the importance of the aspect ratio in the design of proper sensors or actuators. Furthermore, the negative coefficient of the cross-density distribution causes the material to behave much more similarly to the hydrogels with Nv=0.001, which is much softer than the hydrogels with Nv=0.01. The cross-link density is the stiffness representation of hydrogels, and the hydrogel with a high cross-link density has a small swelling ratio. Therefore, the negative coefficient of the cross-density distribution makes the hydrogel swell more than the positive one (see Fig. 10).

4 Summary and conclusions

In this paper, a semi-analytical solution for finite bending is developed for FG pH-sensitive hydrogels, which can be used in various applications such as microfluidic and biomedical actuators and sensors. A pH-sensitive hydrogel layer with linear or exponential cross-link density distributions is investigated in response to the external stimuli. Considering the free energy density, which is additively decomposed into several contributions including network stretching, the solvent mixing with networks and ions, and the dissociation of acidic groups, the behaviors of pH-responsive hydrogels are predicted. Assuming the total deformation gradient tensor for mapping the initial state into the deformed one, we develop a second-order nonlinear differential equation through satisfying the equilibrium equation. Regarding the presented boundary condition, we solve the system of equations with an iterative procedure. The optimization approach is used in each iteration to calculate the semi-angle through vanishing the net tangential force and moment component. In order to investigate the accuracy and performance of the proposed method, the stress and deformation fields are compared with those of the FEM in some cases. The effects of the cross-link density variation are studied on some design factors, e.g., the deformation and stress fields, the semi-angle, the bending curvature, the aspect ratio, and the neutral axes. The results show that the stress, deformation, and bending curvature increase when the pH value of the exterior solvent increases. Continuous deformation and stress fields are observed due to the continuous distribution of the cross-link density. However, the multi-layer hydrogels are discontinuous. In addition, the aspect ratio of an FG layer affects the semi-angle dramatically, but almost has no influence on the bending curvature, stress, and deformation distribution. The proposed semi-analytical method can be used in the applications such as sensors and actuators, and the semi-analytical solution can be used in the design and optimization process of smart structures.