Coupling behavior of the pH/temperature sensitive hydrogels for the inhomogeneous and homogeneous swelling

Smart hydrogels consist of covalently cross-linked networks which are capable of responding to the environmental conditions through the swelling phenomenon. The major reaction of the network against the environmental conditions is change in the swelling amount. Some of the effective parameters in the amount of swelling are temperature [1–3], pH [4, 5], mechanical load [6, 7], ionic and salt concentration [8], light [9, 10] and electric field [11, 12]. Due to this interesting behavior, these materials have found novel applications in tissue engineering and drug delivery systems [13], sensors and actuators [14–18], microfluidics and micro-valves [19, 20]. In most of applications, a constitutive model is required to describe the hydrogel response under complicated thermo-chemo-mechanical conditions.

The coupled nonlinear field theory of hydrogels has attracted a great deal of interest, especially in recent years [4, 6, 8, 21–26]. A basic work in this field is that of Hong et al [23] in which a constitutive model was presented for neutral hydrogels. This work provided a framework for other researcher study in this field. On the other hand, in line with their previous work, Hong et al [27] developed a FEM tools for modeling the swelling behavior of neutral hydrogels, successfully. Following these pioneering works, the behavior of the neutral hydrogels was studied by other researcher. The coupled thermo-mechanical behavior of the temperature-sensitive hydrogels has been investigated in recent years. Ji et al [28] developed a constitutive model in the scope of finite-strain kinematics for thermally induced volume transitions in stimulus–responsive hydrogels. Employing a biphasic approach, Birgersson et al [29] presented a model for transient deformation of a neutral hydrogel due to the temperature gradient. Chester and Anand [30] introduced a large deformation constitutive model for the coupled thermo-mechanical behavior of thermally responsive elastomeric gels. Validating with the experimental data, Cai and Suo [31] established a constitutive model for the temperature sensitive poly-(N-isopropylacrylamide) hydrogel (PNIPAM hydrogel). However, in the vicinity of the phase-transition temperature, multiple-solutions and snap-through instability was observed. Following Cai and Suo [31], Mazaheri et al [32] presented a modified model for the PNIPAM hydrogel which eliminates the instability and successfully implemented the modified model in inhomogeneous swelling problems [20].

Furthermore, complicated behavior of the pH-sensitive hydrogels has been studied in recent years [4, 5, 33, 34]. Marcombe et al [4] developed a model for the pH-sensitive hydrogels at the equilibrium state. They successfully employed electro-neutrality assumption for the hydrogel network and the external solution. Comparing to the experimental results, they validated their model. Implementing the model in the finite element framework, they discussed on some examples. Also, Kurnia et al [35] studied the behavior of pH-sensitive hydrogels and presented the solution for the swelling behavior of a cylinder made of the pH-sensitive hydrogel. They employed conservation laws for mass, momentum and ion, and validated the results with the available experimental data. Li et al [36] investigated the constitutive modeling of polyelectrolyte gels and verified the model using a polyacrylamide-co-acrylic acid hydrogel. They utilized the Gent model for capturing the non-Gaussian chain effect and validated the model through comparing with the experiments. Following Marcombe et al [4], Yan et al [34] developed a field theory for the behavior of a network of polyampholytic pH-sensitive hydrogel. They studied the influence of pH, salt concentration, geometric constraint and the effects of the electrical double layer on the swelling properties.

In this work following Guo and Gao [37] and Mazaheri et al [32], we develop a field theory for the swelling behavior of the pH/temperature-responsive carboxymethyl chitosan/poly(N-isopropylacrylamide)semi-IPN hydrogels. We consider this material as a PNIPAM hydrogel with some acidic groups attached to it. The interaction parameter of the network is assumed to be a function of both the temperature and polymer volume fraction as reported for the PNIPAM hydrogels by Afroze et al [38]. Also, we employ the electro-neutrality assumption for the electric charge distribution. The presented model is validated with the experimental results available in the literature. It is noteworthy to mention that the multiple-solutions and instability problems in the vicinity of PTT are also avoided.

The paper is organized as follows. Firstly, the model is explained in detail and governing formulations are presented in section 2. Thereafter, homogeneous swelling of the pH/temperature sensitive hydrogels is investigated in section 3.1. In section 3.1.1, a comprehensive study on the free swelling problem is presented and the model is validated comparing to the experimental data reported by Guo and Gao [37]. As another example of the homogeneous deformation, the biaxial constrained swelling problem of the hydrogel is studied in section 3.1.2. Then, the inhomogeneous swelling of the spherical hydrogel shells on a fixed rigid core is analytically studied for pH/temperature sensitive hydrogel and the results are presented in section 3.2.1. Also, numerical implementation of the presented model in ABAQUS is performed. With the aid of this numerical tool, considering the contact effects, the inhomogeneous swelling of the spherical hydrogel shells on a fixed rigid core is investigated in section 3.2.2. Finally, we present a summary and draw conclusions in section 4.

In this paper, the pH/temperature sensitive hydrogel network is composed of polymeric chains bearing acidic groups $AH.$ Immersing in an aqueous solution, the network imbibes the solution and swells. Due to existence of ions in the solutions and acidic groups on the network, the swelling of the network is affected by the external solutions conditions, specially its pH and ionic strength. Also, due to existence of semi-NIPAM chain in the network, the network is capable of responding to the temperature changes. A schematics draw of this material structure beside the external solution is shown in figure 1.

Zoom In Zoom Out Reset image size

Diffusion of the external solvent molecules and ions through the hydrogel, results in deformation in the network. Considering an element of the network in the reference (dry) and current (deformed) states, the coordinates are denoted by vectors ${\bf{X}}$ and ${\bf{x}}({\bf{X}}),$ respectively. Employing the Lagrangian approach, the current coordinate, ${\bf{x}}({\bf{X}}),$ is a function of its reference one which describes the network deformation. In terms of ${\bf{x}}({\bf{X}}),$ the deformation gradient tensor, ${\bf{F}},$ and consequently, the right Cauchy-Green deformation tensor, ${\bf{C}},$ can be defined as:

$$\begin{eqnarray}&&{\bf{F}}=\displaystyle \frac{\partial {\bf{x}}}{\partial {\bf{X}}},\;\;\;\;{\bf{C}}={{\bf{F}}}^{T}{\bf{F}},\end{eqnarray} \tag{ 1 }$$

where the superscript 'T' denotes the transpose of the tensor. The region occupied by the network in the references state and its boundary are denoted by $V$ and ${\rm{\Gamma }},$ respectively. As shown in figure 1, in the equilibrium state, the external solution is composed of four species, i.e. solvent molecules, counter-ion, co-ion and hydrogen ion which their number in the external solution are ${\bar{n}}_{s},$ ${\bar{n}}_{+},$ ${\bar{n}}_{-},$ and ${\bar{n}}_{{H}^{+}},$ respectively. Assuming small variations in the species numbers, the change in the free energy of the external solution at a specified temperature, $T,$ is:

$$\begin{eqnarray}&&{\mu }_{s}\delta {\bar{n}}_{s}+{\mu }_{{H}^{+}}\delta {\bar{n}}_{{H}^{+}}+{\mu }_{+}\delta {\bar{n}}_{+}+{\mu }_{-}\delta {\bar{n}}_{-},\end{eqnarray} \tag{ 2 }$$

where ${\mu }_{s},$ ${\mu }_{+},$ ${\mu }_{-}$ and ${\mu }_{{H}^{+}}$ are electro-chemical potential of the related species, respectively and $\delta$ denotes the variational operator. Also, neglecting body forces, the amount of potential energy change in the network due to the mechanical force is equal to $-{\int }_{\Gamma }{t}_{i}\delta {x}_{i}{\rm{d}}{\rm{\Gamma }}$ where ${\bf{t}}={\bf{P}}\hat{{\bf{n}}}$ is the traction force on the boundary of the network. ${\bf{P}}$ and $\hat{{\bf{n}}}$ are the nominal stress tensor and the boundary unit normal vector, respectively. The change in the free energy of the network is equal to ${\int }_{V}\delta W{\rm{d}}V$ in which $W$ is the free energy density of the network expressed in the reference state. Considering both the hydrogel and the external solution as a system, the condition of the equilibrium is [4]:

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta W{\rm{d}}V+{\mu }_{s}\delta {\bar{n}}_{s}+{\mu }_{{H}^{+}}\delta {\bar{n}}_{{H}^{+}}+{\mu }_{+}\delta {\bar{n}}_{+}+{\mu }_{-}\delta {\bar{n}}_{+}\\ &&\quad -\displaystyle {\int }_{\Gamma }{t}_{i}\delta {x}_{i}{\rm{d}}{\rm{\Gamma }}=0.\end{eqnarray} \tag{ 3 }$$

On account of conservation of any species in the system, we may write:

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta {C}_{s}({\bf{X}}){\rm{d}}V+\delta {\bar{n}}_{s}=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta {C}_{+}({\bf{X}}){\rm{d}}V+\delta {\bar{n}}_{+}=0,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta {C}_{-}({\bf{X}}){\rm{d}}V+\delta {\bar{n}}_{-}=0,\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta {C}_{{H}^{+}}({\bf{X}}){\rm{d}}V+\delta {\bar{n}}_{{H}^{+}}=\displaystyle {\int }_{V}\delta {C}_{{A}^{-}}({\bf{X}}){\rm{d}}V,\end{eqnarray} \tag{ 7 }$$

where ${C}_{i}$ is nominal concentration of species $i$ in the network and is related to the true concentration, ${c}_{i},$ via ${C}_{i}={c}_{i}\mathrm{det}({\bf{F}}).$ Also, employing the electro-neutrality both in the network and the external solution we arrive at:

$$\begin{eqnarray}&&{C}_{{H}^{+}}({\bf{X}})+{C}_{+}({\bf{X}})={C}_{{A}^{-}}({\bf{X}})+{C}_{-}({\bf{X}}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\bar{n}}_{+}+{\bar{n}}_{{H}^{+}}={\bar{n}}_{-}.\end{eqnarray} \tag{ 9 }$$

Substituting (4)–(9) in (3), we attain a variational statement for the free energy change of the system as:

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\delta W{\rm{d}}V-\displaystyle {\int }_{\Gamma }{t}_{i}\delta {x}_{i}{\rm{d}}A-\displaystyle {\int }_{V}\delta {C}_{s}({\bf{X}}){\mu }_{s}{\rm{d}}V\\ &&\quad -\displaystyle {\int }_{V}\delta {C}_{-}({\bf{X}})\{{\mu }_{-}+{\mu }_{{H}^{+}}\}{\rm{d}}V\\ &&\quad -\displaystyle {\int }_{V}\delta {C}_{+}({\bf{X}})\{{\mu }_{+}-{\mu }_{{H}^{+}}\}{\rm{d}}V=0.\end{eqnarray} \tag{ 10 }$$

During the acidic dissociation, the number of the fixed acidic groups on the network remains constant and only the hydrogen ions can freely move between the external solution and the network. Thus, we have [4]:

$$\begin{eqnarray}&&{C}_{{A}^{-}}({\bf{X}})+{C}_{AH}({\bf{X}})=f/\nu ,\end{eqnarray} \tag{ 11 }$$

where $f$ and $\nu$ are number of the acidic groups attached to the network per monomer and volume of a monomer, respectively. It is concluded from (8) and (11) that both the ${C}_{{A}^{-}}({\bf{X}})$ and ${C}_{AH}({\bf{X}})$ are not independent variables and can be stated in terms of ${C}_{{H}^{+}}({\bf{X}}),$ ${C}_{-}({\bf{X}})$ and ${C}_{+}({\bf{X}}).$ Thus, considering the free energy function as $W=W({\bf{F}},{C}_{s},{C}_{{H}^{+}},{C}_{-},{C}_{+}),$ its variation statement is calculated as $\delta W=\tfrac{\partial W}{\partial {\bf{F}}}{\boldsymbol{:}}\delta {\bf{F}}+\tfrac{\partial W}{\partial {C}_{s}}\delta {C}_{s}+\tfrac{\partial W}{\partial {C}_{{H}^{+}}}\delta {C}_{{H}^{+}}\;+\tfrac{\partial W}{\partial {C}_{-}}\delta {C}_{-}+\tfrac{\partial W}{\partial {C}_{+}}\delta {C}_{+}.$ Now, considering the free energy variation and employing divergence theorem for the second term of (10), we arrive at:

$$\begin{eqnarray}&&\displaystyle {\int }_{V}\left\{\left(\displaystyle \frac{\partial W}{\partial {\bf{F}}}-{\bf{P}}\right){\boldsymbol{:}}\delta {\bf{F}}+\left(\displaystyle \frac{\partial W}{\partial {C}_{s}}-{\mu }_{s}\right)\delta {C}_{s}\right.\\ &&\quad +\displaystyle \frac{\partial W}{\partial {C}_{{H}^{+}}}\delta {C}_{{H}^{+}}+\left(\displaystyle \frac{\partial W}{\partial {C}_{-}}-({\mu }_{-}+{\mu }_{{H}^{+}})\right)\delta {C}_{-}\\ &&\left.\quad +\left(\displaystyle \frac{\partial W}{\partial {C}_{+}}-({\mu }_{+}-{\mu }_{{H}^{+}})\right)\delta {C}_{+}\right\}{\rm{d}}V=0.\end{eqnarray} \tag{ 12 }$$

Therefore, following the standard arguments, the constitutive laws for the network are expressed as:

$$\begin{eqnarray}&&{\bf{P}}=\displaystyle \frac{\partial W}{\partial {\bf{F}}},\;\;\;{\mu }_{s}=\displaystyle \frac{\partial W}{\partial {C}_{s}},\;\;\;{\mu }_{-}+{\mu }_{{H}^{+}}=\displaystyle \frac{\partial W}{\partial {C}_{-}},\\ &&{\mu }_{+}-{\mu }_{{H}^{+}}=\displaystyle \frac{\partial W}{\partial {C}_{+}},\;\;\;\displaystyle \frac{\partial W}{\partial {C}_{{H}^{+}}}=0.\end{eqnarray} \tag{ 13 }$$

Neglecting the volume of the other species in the network due to their low concentration, the change in the volume of the network originates from the fluid molecules diffusion through the network. Thus, we have [23]:

$$\begin{eqnarray}&&1+{\nu }_{s}{C}_{s}=\mathrm{det}({\bf{F}})=J,\end{eqnarray} \tag{ 14 }$$

where ${\nu }_{s}$ and $J$ are the volume of a molecule of the fluid and deformation gradient determinant, respectively. To implement this constraint on the system, employing the Lagrange method, we should modify the free energy as $\tilde{W}=W+{\rm{\Pi }}(1+{\nu }_{s}{C}_{s}-J).$ In terms of the modified free energy density, the constitutive laws for the nominal stress and chemical potential of solvent molecules are recalculated as:

$$\begin{eqnarray}&&{\bf{P}}=\displaystyle \frac{\partial \tilde{W}}{\partial {\bf{F}}}=\displaystyle \frac{\partial W}{\partial {\bf{F}}}-{\rm{\Pi }}\displaystyle \frac{\partial J}{\partial {\bf{F}}},\;\;\;{\mu }_{s}=\displaystyle \frac{\partial W}{\partial {C}_{s}}+{\rm{\Pi }}{\nu }_{s},\end{eqnarray} \tag{ 15 }$$

where ${\rm{\Pi }}$ is the Lagrange multiplier and its elimination from (15) results in:

$$\begin{eqnarray}&&{\bf{P}}=\displaystyle \frac{\partial \tilde{W}}{\partial {\bf{F}}}=\displaystyle \frac{\partial W}{\partial {\bf{F}}}-\displaystyle \frac{1}{{\nu }_{s}}\left({\mu }_{s}-\displaystyle \frac{\partial W}{\partial {C}_{s}}\right)\displaystyle \frac{\partial J}{\partial {\bf{F}}}.\end{eqnarray} \tag{ 16 }$$

Considering a small density of crosslinking for the hydrogel, we assume an additive decomposition of the free energy density in the form of (see e.g., [4, 23, 27, 30, 31, 39] among others):

$$\begin{eqnarray}&&W={W}_{{\rm{elastic}}}({\bf{F}},T)+{W}_{{\rm{mix}}}({C}_{s},T)+{W}_{{\rm{ion}}}({\bf{F}},{C}_{+},{C}_{-},{C}_{{H}^{+}},T)\\ &&\quad \quad +\;{W}_{{\rm{diss}}}({C}_{AH},{C}_{{A}^{-}},T),\end{eqnarray} \tag{ 17 }$$

where ${W}_{{\rm{elastic}}},$ ${W}_{{\rm{mix}}},$ ${W}_{{\rm{ion}}}$ and ${W}_{{\rm{diss}}}$ are the free energy density due to the elastic deformation, mixing between the solvent and the network chain, mixing of ions and dissociation of the acidic groups, respectively. We use a neo-Hookean model to describe the elastic deformation of the network where it is expressed in terms of the deformation as:

$$\begin{eqnarray}&&{W}_{{\rm{elastic}}}=\displaystyle \frac{1}{2}NKT({I}_{1}-3\;-\;2\;\mathrm{log}(J)),\end{eqnarray} \tag{ 18 }$$

in which $N$ and $K$ are density of the polymer chains in the reference state and Boltzmann constant, respectively. In addition, ${I}_{1}$ is the first invariant of ${\bf{C}}.$ Also, based on the Florry–Huggins theory, the mixing part of the free energy density is [40, 41]:

$$\begin{eqnarray}&&{W}_{{\rm{mix}}}=\displaystyle \frac{KT}{{\nu }_{s}}\left({\nu }_{s}{C}_{s}\;\mathrm{ln}\left(\displaystyle \frac{{\nu }_{s}{C}_{s}}{1+{\nu }_{s}{C}_{s}}\right)+\displaystyle \frac{{\nu }_{s}{C}_{s}\chi }{1+{\nu }_{s}{C}_{s}}\right),\end{eqnarray} \tag{ 19 }$$

where $\chi ={\chi }_{0}+\phi {\chi }_{1}={A}_{0}+{B}_{0}T+\phi ({A}_{1}+{B}_{1}T)$ is the interaction parameter of the network assumed to be a function of the temperature and the polymer volume fraction, $\phi =1/J,$ for the temperature sensitive PNIPAM hydrogel. In addition, ${A}_{0},$ ${B}_{0},$ ${A}_{1}$ and ${B}_{1}$ are material parameters adopted for the experiments reported by Afroze et al [38]. In light of the swelling constraint, (14), we rewrite (19) as:

$$\begin{eqnarray}&&{W}_{{\rm{mix}}}=\displaystyle \frac{KT}{{\nu }_{s}}(J-1)\left(\mathrm{ln}\left(\displaystyle \frac{J-1}{J}\right)+\displaystyle \frac{\chi }{J}\right).\end{eqnarray} \tag{ 20 }$$

However, as discussed by Mazaheri et al [32] this part of the free energy density encounters a numerical instability in the vicinity of phase transition temperature due to multiple solution problems. It is worth to mention that some techniques such as Riks methods can be used to overcome this instability however, this is computationally expensive in a numerical scheme [27, 42]. To avoid such numerical difficulties, we should modify the mixing energy statement so that the instability is removed while the model accuracy is preserved. To this end the method presented by Mazaheri et al [32] is utilized. This instability is eliminated via polynomial expanding of the ill-posed logarithmic term as [32]:

$$\begin{eqnarray}&&{W}_{{\rm{mix}}}=\displaystyle \frac{KT}{{\nu }_{s}}(J-1)\left(\mathrm{ln}\left(\displaystyle \frac{J-1}{J}\right)+\displaystyle \frac{\chi }{J}\right)\simeq \displaystyle \frac{KT}{{\nu }_{s}}(J-1)\\ &&\quad \quad \quad \times \left(-\displaystyle \frac{1}{J}-\displaystyle \frac{1}{2{J}^{2}}-\displaystyle \frac{1}{3{J}^{3}}+\displaystyle \frac{\chi }{J}\right).\end{eqnarray} \tag{ 21 }$$

Employing this free energy, a stable model is developed which is capable of predicting the transition behavior of the hydrogel for temperature changes. As a results, we can use this model in the inhomogeneous problems without encountering any inherent instability of the model with a lower computational cost. Considering a low concentration of the ions, the change in the free energy density of the system due to the ions mixing is entropic and expressed as [4]:

$$\begin{eqnarray}{W}_{{\rm{ion}}} & = & KT\left({C}_{{H}^{+}}\left(\mathrm{ln}\left(\displaystyle \frac{{C}_{{H}^{+}}}{{c}_{{H}^{+}}^{{\rm{ref}}}\mathrm{det}({\bf{F}})}\right)-1\right)\right.\\ & & +{C}_{-}\left(\mathrm{ln}\left(\displaystyle \frac{{C}_{-}}{{c}_{-}^{{\rm{ref}}}\mathrm{det}({\bf{F}})}\right)-1\right)\\ & & \left.+{C}_{+}\left(\mathrm{ln}\left(\displaystyle \frac{{C}_{+}}{{c}_{+}^{{\rm{ref}}}\mathrm{det}({\bf{F}})}\right)-1\right)\right),\end{eqnarray} \tag{ 22 }$$

where ${c}_{d}^{{\rm{ref}}}$ is a reference concentration for the species $d.$ Also, the free energy change of the acidic dissociation which originates from the entropy change due to mixing and enthalpy change for the acidic dissociation, is stated as [4]:

$$\begin{eqnarray}&&{W}_{{\rm{diss}}}=KT\left({C}_{{A}^{-}}\;\mathrm{ln}\left(\displaystyle \frac{{C}_{{A}^{-}}}{{C}_{{A}^{-}}+{C}_{AH}}\right)+{C}_{AH}\;\mathrm{ln}\right.\\ &&\quad \quad \quad \;\times \left.\left(\displaystyle \frac{{C}_{AH}}{{C}_{{A}^{-}}+{C}_{AH}}\right)\right)+\gamma {C}_{{A}^{-}}\\ &&\quad \quad =KT\left({C}_{{A}^{-}}\;\mathrm{ln}\left(\displaystyle \frac{{C}_{{A}^{-}}}{f/\nu }\right)+{C}_{AH}\;\mathrm{ln}\right.\\ &&\quad \quad \quad \;\times \left.\left(\displaystyle \frac{{C}_{AH}}{f/\nu }\right)\right)+\gamma {C}_{{A}^{-}},\end{eqnarray} \tag{ 23 }$$

where $\gamma$ is the enthalpy change due to the dissociation of one acidic group. Employing the defined free energy density and constitutive laws, we arrive at:

$$\begin{eqnarray}&&\displaystyle \frac{\partial W}{\partial {C}_{{H}^{+}}}=0=\displaystyle \frac{\partial {W}_{{\rm{ion}}}}{\partial {C}_{{H}^{+}}}+\displaystyle \frac{\partial {W}_{{\rm{diss}}}}{\partial {C}_{{H}^{+}}}.\end{eqnarray} \tag{ 24 }$$

Keeping in mind that both ${C}_{{A}^{-}}({\bf{X}})$ and ${C}_{AH}({\bf{X}})$ are not independent, through applying the chain rule in (24), we have:

$$\begin{eqnarray}&&{N}_{A}{K}_{a}={c}_{{H}^{+}}^{{\rm{ref}}}\mathrm{exp}\left(-\displaystyle \frac{\gamma }{KT}\right).\end{eqnarray} \tag{ 25 }$$

Now, based on (13), the chemical potential of the co-ion is:

$$\begin{eqnarray}&&{\mu }_{-}+{\mu }_{{H}^{+}}=\displaystyle \frac{\partial W}{\partial {C}_{-}}=\displaystyle \frac{\partial {W}_{{\rm{ion}}}}{\partial {C}_{-}}+\displaystyle \frac{\partial {W}_{{\rm{diss}}}}{\partial {C}_{-}}.\end{eqnarray} \tag{ 26 }$$

Substituting (25) in (26) and after some mathematical manipulation we obtain:

$$\begin{eqnarray}&&{\mu }_{-}+{\mu }_{{H}^{+}}=KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{{H}^{+}}}{{c}_{{H}^{+}}^{{\rm{ref}}}}\right)+KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{-}}{{c}_{-}^{{\rm{ref}}}}\right)\\ &&\quad \quad \quad \quad \quad =\;KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{{H}^{+}}{c}_{-}}{{c}_{{H}^{+}}^{{\rm{ref}}}{c}_{-}^{{\rm{ref}}}}\right).\end{eqnarray} \tag{ 27 }$$

Similarly, for the chemical potential of counter-ion we deduce:

$$\begin{eqnarray}&&{\mu }_{+}-{\mu }_{{H}^{+}}=\displaystyle \frac{\partial W}{\partial {C}_{-}}=KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{+}}{{c}_{+}^{{\rm{ref}}}}\right)-KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{{H}^{+}}}{{c}_{{H}^{+}}^{{\rm{ref}}}}\right)\\ &&\quad \quad \quad \quad \quad =\;KT\;\mathrm{ln}\left(\displaystyle \frac{{c}_{+}{c}_{{H}^{+}}^{{\rm{ref}}}}{{c}_{+}^{{\rm{ref}}}{c}_{{H}^{+}}}\right).\end{eqnarray} \tag{ 28 }$$

Assuming the dilute solution, the chemical potential of species $d$ in the solution is [8]:

$$\begin{eqnarray}&&{\mu }_{d}=KT\;\mathrm{ln}\left(\displaystyle \frac{{\bar{c}}_{d}}{{c}_{d}^{{\rm{ref}}}}\right),\end{eqnarray} \tag{ 29 }$$

in which ${\bar{c}}_{d}$ is the true concentration for the species $d.$ Substituting (29) in (27) and (28) yields:

$$\begin{eqnarray}&&\displaystyle \frac{{c}_{-}}{{\bar{c}}_{-}}=\displaystyle \frac{{\bar{c}}_{{H}^{+}}}{{c}_{{H}^{+}}},\;\;\;\displaystyle \frac{{c}_{+}}{{\bar{c}}_{+}}=\displaystyle \frac{{c}_{{H}^{+}}}{{\bar{c}}_{{H}^{+}}},\end{eqnarray} \tag{ 30 }$$

which are well-known equations named the Donnan equations. Based on (30), it is concluded that the co-ion and counter-ion true concentrations in the network are dependent on the true concentration of the hydrogen ion in the network and external solution conditions. Thus, once the external solution conditions and hydrogen ion concentration in network are known, dependent values of the co-ion and counter-ion true concentrations in the network are determined.

The acidic dissociation obeys a relation as below:

$$\begin{eqnarray}&&\displaystyle \frac{[{H}^{+}][{A}^{-}]}{[AH]}={K}_{a},\end{eqnarray} \tag{ 31 }$$

where $[i]$ represents the molarity of species $i$ in the network which is equal to ${c}_{i}=[i]{N}_{A}$ while ${N}_{A}=6.023\times {10}^{23}$ is the Avogadro number. Now, substituting (8) and (11) in (31) and expressing the concentration in the current state we arrive at:

$$\begin{eqnarray}&&\displaystyle \frac{{c}_{{A}^{-}}{c}_{{H}^{+}}}{{c}_{AH}}={N}_{A}{K}_{a}=\displaystyle \frac{{c}_{{H}^{+}}({c}_{{H}^{+}}+{c}_{+}-{c}_{-})}{\displaystyle \frac{f}{vJ}-({c}_{{H}^{+}}+{c}_{+}-{c}_{-})}.\end{eqnarray} \tag{ 32 }$$

Also, replacing (30) in (32) results in a cubic equation for ${c}_{{H}^{+}}$ as below:

$$\begin{eqnarray}&&{A}_{1}{(\nu {c}_{{H}^{+}})}^{3}+{N}_{A}{K}_{a}\nu {A}_{1}{(\nu {c}_{{H}^{+}})}^{2}\\ &&\quad +\left(-{A}_{2}-{N}_{A}{K}_{a}\displaystyle \frac{{A}_{3}}{J}\right)(\nu {c}_{{H}^{+}})\\ &&\quad -{N}_{A}{K}_{a}\nu {A}_{2}=0,\end{eqnarray} \tag{ 33 }$$

where the coefficients ${A}_{i}s$ are:

$$\begin{eqnarray}&&{A}_{1}=\left(1+\displaystyle \frac{{\bar{c}}_{+}}{{\bar{c}}_{{H}^{+}}}\right),\;\;{A}_{2}={\bar{c}}_{{H}^{+}}{\bar{c}}_{-}{\nu }^{2},\;\;{A}_{3}=f\nu .\end{eqnarray} \tag{ 34 }$$

As concluded from Donnan equations, the co-ion and counter-ion concentrations in the network are dependent on the ${c}_{{H}^{+}}$ which could be determined from (33) and are functions of the deformation due to presence of $J=\mathrm{det}({\bf{F}})$ in the coefficient of this cubic equation. Thus, all variables can be expressed as functions of the deformation field and external solution concentrations. It is assumed that the external solution is large enough so that the salt concentration and pH are assumed to be constant during the deformation of the hydrogel network. The hydrogen ion concentration in the external solution is calculated as: ${\bar{c}}_{{H}^{+}}={N}_{A}[{H}^{+}]={N}_{A}{10}^{-\overline{pH}}.$ Also, having the counter-ion concentration in the external solution, the related co-ion concentration can be computed from the electro-neutrality assumption as: ${\bar{c}}_{-}={\bar{c}}_{+}+{\bar{c}}_{{H}^{+}}.$ Now, having the external concentration, the hydrogen ion, the co-ion and the counter-ion concentrations in the hydrogel network can be determined from the Donnan equation (30) and the cubic equation of (33). Thus, the complicated multi-field swelling problem is reduced to a deformation field problem. In order to investigate the deformation field, a statement for the stress is necessary. In this paper, we use the nominal stress determined from the constitutive laws. After some simplifications we have:

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{ij}\nu }{KT}=N\nu ({F}_{ij}-{H}_{ij})+\left(\displaystyle \frac{-\mathrm{1/2}+({\chi }_{0}-{\chi }_{1})}{{J}^{2}}\right.\\ &&\quad \quad \quad \quad \;+\;\left.\displaystyle \frac{-\mathrm{1/3+2}{\chi }_{1}}{{J}^{3}}-\displaystyle \frac{1}{{J}^{4}}\right)J{H}_{ij}\\ &&\quad \quad \quad \quad \;-\;({c}_{{H}^{+}}+{c}_{+}+{c}_{-}-{\bar{c}}_{{H}^{+}}-{\bar{c}}_{+}-{\bar{c}}_{-})J{H}_{ij}.\end{eqnarray} \tag{ 35 }$$

where ${\bf{H}}=\tfrac{1}{J}\tfrac{\partial J}{\partial {\bf{F}}}={{\bf{F}}}^{-T}.$ It is worth to mention that in this work we present a constitutive model for pH and temperature sensitive hydrogels. To account for the pH sensitivity, we employ the well-known model of Marcombe et al [4]. But to capture the temperature sensitive behavior of the model, we implement the mixing parameter function in the free energy to model the temperature sensitivity. Also, to achieve a stable behavior of the model especially in the vicinity of the phase transition temperature, the mixing function is modified. Thus, the model is capable of simulating the behavior of the pH and temperature sensitive hydrogels.

In this section, the proposed model is employed to investigate some bench-mark problems involving homogeneous and/or inhomogeneous swelling of the pH/temperature sensitive hydrogels. The material parameters of the network are adopted for experiments available in the literature and shown in table 1. Based on the other works [4, 23, 37], we assume that $\nu ={\nu }_{s}={10}^{-28},$ $KT=4\times {10}^{-21}\;J$ and ${K}_{a}={10}^{-4.7}.$ Firstly in section 3.1, some homogeneous examples are solved and the results are compared to the experimental data available in the literature to show the validity of the model. Then in section 3.2, inhomogeneous swelling of the hydrogel is studied for a spherical shell of the pH/temperature sensitive hydrogel on a hard core for variation of the both pH and temperature. This problem is solved both analytically and numerically. To solve boundary value problems the proposed model is implemented in the finite element framework through producing a user defined subroutine in ABAQUS software.

Table 1.  Material parameters for the network [38].

${A}_{0}$	${B}_{0}$	${A}_{1}$	${B}_{1}$
−12.947	0.044 96 K−1	17.92	−0.0569 K−1

3.1. Homogeneous deformations

3.1.1. Free swelling

In this section, we investigate the free swelling of pH/temperature sensitive hydrogels to validate the model with the experimental results available in the literature. Due to the homogenous deformation in the free swelling problem, we have:

$$\begin{eqnarray}&&{\lambda }_{1}={\lambda }_{2}={\lambda }_{3}=\lambda ={\left(\displaystyle \frac{V}{{V}_{0}}\right)}^{\mathrm{1/3}},\end{eqnarray} \tag{ 36 }$$

where $V/{V}_{0}$ is the swelling ratio defined as the ratio of the current volume $V,$ to its value at the dry state, ${V}_{0}.$ In the free swelling problem, the network experiences a free-stress state. Thus, to determine the equilibrium state for the free swelling state, equating the stress equation, namely (35), to zero, the value of swelling ratio at the equilibrium is identified and the results are shown in figure 2 for both the temperature and pH variations. For comparison we present the results of the free swelling of the neutral PNIPAM hydrogel where $Nv=0.04$ and $f=0.00$ in figure 2(a). As depicted in this figure, the neutral PNIPAM has no sensitivity to the pH variation but in contrary is sensitive to the temperature variation. In figures 2(b) and (c), the results are shown for $f=0.05$ and $f=0.10,$ respectively, while $Nv=0.04$ and the salt molarity is $0.001\;M$ for both cases. The results demonstrate that increasing the acidic group density on the network, makes it more sensitive to the pH variations. To investigate the effect of salt concentration, the free swelling results are presented for the external solutions with salt molarity of $0.001\;M$ and $0.050\;M$ in figures 2(b) and (d), respectively, in which $Nv=0.04$ and $f=0.10.$ As observed from figures 2(b) and (d), the swelling ratio decreases by increase in the molarity of the external solution due to the diffusion of more ions in the network space and consequently, repulsion of some solvent molecules from the network.

Zoom In Zoom Out Reset image size

To study the effect of the crosslinking density, namely $Nv,$ the relevant results are shown for $Nv=0.03$ and $Nv=0.02$ in figures 2(e) and (f), respectively, where the salt molarity and acidic group density are constant and equal to $0.001\;M$ and $f=0.10$ for both cases. As illustrated in figures 2(b), (e) and (f), the network is more sensitive to pH/temperature for the smaller values of $Nv$ and larger volume change is observed due to the less resistance of the network against deformation. As depicted in all parts of figure 2, the swelling ratio decreases as the temperature increases due to PNIPAM content of the network. But, increasing the pH value results in a decrease in the hydrogen ion and, as a result, more acidic groups are going to be dissociated and consequently more hydrophilic regions are developed and the chains are stretched due to electrical repulsion of the negative dissociated acidic groups. Thereupon, more fluid molecules are absorbed in the network and larger amounts of the volume change are observed. The attained results are in agreement with the experimental data available in the literature [37].

To validate the model with the experimental data, we consider the experiments reported by Guo and Gao [37]. But the required parameters have not been reported in the mentioned reference. Thus, we should fit the values of $N\nu ,$ $f$ and the external solution molarity for the experiments. First, we determine the value of $N\nu =0.04$ for the swelling behavior of neutral PNIPAM when $f=0.00.$ Then, $f$ and the external solution molarity should be found through fitting for pH values of 2.1 and 7.4 as reported by Guo and Gao [37]. The fitted value of the external solution molarity that is kept constant in all experiments is equal to $0.001\;M.$ The value of $f$ which should be fitted for every experimental data set, is computed as $f=0.06,\;\;0.12$ and $0.18$ for the best-fit.

As shown in figure 3, the presented model predicts qualitatively, the swelling behavior of the pH/temperature sensitive hydrogel. Also, below the phase transition temperature, a good agreement is observed between the present work results with those of the experiments. As one may see from this figure, at the small values of pH, the network is not that sensitive to the acidic group density due to the large amount of hydrogen ion in the external solution which prevents from the acidic groups dissociation. Thus, the acidic group density has almost no effect on the swelling behavior of the network as predicted by the model and confirmed by the experiments. But, for the larger values of pH, scarcity of the hydrogen ion results in the acidic group dissociation and the larger values of swelling ratio. This phenomenon is more prominent for larger amounts of the acidic group density as shown in figure 3(a). In this case, a network with more amounts of the acidic groups produces more hydrogen ions which results in more fluid uptake in the network and consequently more swelling. As depicted in figure 3, the presented model can predict the pH/temperature behavior qualitatively both in the small and large values of pH. Also, the model is capable of predicting the hydrogel behavior for the temperature changes.

Zoom In Zoom Out Reset image size

3.1.2. Biaxial constrained swelling

As investigated by many researchers, the swelling of a hydrogel layer on a hard sub-layer is studied in this section. The hydrogel layer is assumed to be in the biaxial constrained conditions due to its small thickness. The pre-stretches of the hydrogel layer in two constrained directions are fixed and equal to ${\lambda }_{1}={\lambda }_{2}={\lambda }_{0}$ where ${\lambda }_{0}$ is the equilibrium stretch at a reference point in which the hydrogel has a stress-free state. Also, the stretch in the direction normal to the constrained plane is ${\lambda }_{3}=\lambda .$ Substituting the pre-stretch in the stress expression and equating it to zero, the value of swelling ratio at the equilibrium is calculated. The equilibrium point of the hydrogel is found in terms of the pH and temperature changes and shown in figures 4(a) and (b), respectively, in which the free swelling results are available for comparison. The material parameters are selected as $N\nu =0.04$ and $f=0.10$ where the external solution molarity is $0.001\;M.$ We assume the reference point with pH = 2 and T = 290 K and as a result, the corresponding value of ${\lambda }_{0}$ is obtained as 2.190. As depicted in figure 4(a), increasing the pH value leads to increase in the swelling ratio for both the free swelling and biaxial constrained cases. But, for the biaxial constrained the swelling ratio is smaller than the free swelling case due to presence of a constraint that prevents the hydrogel from swelling. The effect of increasing the temperature is shown in figure 4(b). As illustrated in this figure, increasing the temperature results in reduction in volume of the hydrogel both for the free and biaxial constrained swelling. However, the amount of volume reduction is larger for the free swelling state due to absence of any constraint against shrinkage of the network. Thus, despite of figure 4(a) results, for increasing the temperature from the reference point, the swelling ratio for the biaxial constrained swelling is larger than the free swelling case due to preventing the hydrogel from shrinking.

Zoom In Zoom Out Reset image size

3.2. Inhomogeneous deformations

In this paper, a modified version of mixing part of the free energy is used and as a result, a continuous model is developed which is numerically stable and can be implemented for solving inhomogeneous swelling problems [32]. To show the capability of the model in the inhomogeneous swelling, we present the solution of the inhomogeneous swelling of a pH/temperature sensitive hydrogel spherical shell on a rigid core. The hard core constrains the hydrogel shell swelling in its inner radius which results in the inhomogeneous deformation of the shell. Analytical and numerical solutions are derived and compared. Thereafter, the numerical framework which is validated with the analytical results is employed to simulate the complicated contact problem of the mentioned hydrogel shell with a flat channel.

3.2.1. Analytical solution

Motivated from Zhao et al [43], a bench-mark problem in the hydrogel inhomogeneous swelling is a spherical shell on a hard core whose volume changes due to the environmental conditions, e.g. pH and temperature changes. The schematic draw of this problem is shown in figure 5 in which the shell is shown at a stress-free state (reference state). In figure 5, ${R}_{{\rm{i}}}$ and ${R}_{{\rm{o}}}$ are inner and outer radii of the spherical shell in the stress-free state, respectively.

Zoom In Zoom Out Reset image size

For the spherical shell, all field variables are functions of $R,$ which is the radius of an element in the stress-free state. All energy statements are recast with respect to the reference state through substituting $J={{\lambda }_{0}}^{3}\bar{J},\;{I}_{1}={{\lambda }_{0}}^{2}{\bar{I}}_{1}$ in the free energy statements in which $\bar{J}$ and ${\bar{I}}_{1}$ are corresponding values of $J$ and ${I}_{1}$ with respect to this state. The temperature and pH at the stress-free state are 290 K and 2, respectively. So, the corresponding value of ${\lambda }_{0}$ is obtained as 2.190. Due to the symmetry, the equilibrium equation in the radial direction is:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{P}_{{\rm{r}}}}{{\rm{d}}R}-\displaystyle \frac{2({P}_{{\rm{r}}}-{P}_{\theta })}{R}=0,\end{eqnarray} \tag{ 37 }$$

where P_r and P_θ are the nominal radial and hoop stresses, respectively. These stress components are obtained from the free energy expressions as:

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{{\rm{r}}}\nu }{KT}=N\nu \left({\lambda }_{0}{\lambda }_{{\rm{r}}}-\displaystyle \frac{1}{{\lambda }_{0}{\lambda }_{{\rm{r}}}}\right)+\left(\displaystyle \frac{-\mathrm{1/2}+({\chi }_{0}-{\chi }_{1})}{{{\lambda }_{0}}^{4}{{\lambda }_{\theta }}^{2}{{\lambda }_{{\rm{r}}}}^{2}}\right.\\ &&\quad \quad \quad \quad \;+\;\left.\displaystyle \frac{-\mathrm{1/3}\;+\;2{\chi }_{1}}{{{\lambda }_{0}}^{7}{{\lambda }_{\theta }}^{4}{{\lambda }_{{\rm{r}}}}^{3}}-\displaystyle \frac{1}{{{\lambda }_{0}}^{10}{{\lambda }_{\theta }}^{6}{{\lambda }_{{\rm{r}}}}^{4}}\right)\\ &&\quad \quad \quad \quad \;-\;({c}_{{H}^{+}}+{c}_{+}+{c}_{-}-{\bar{c}}_{{H}^{+}}-{\bar{c}}_{+}-{\bar{c}}_{-}){{\lambda }_{0}}^{2}{{\lambda }_{\theta }}^{2},\\ &&\displaystyle \frac{{P}_{\theta }\nu }{KT}=N\nu \left({\lambda }_{0}{\lambda }_{\theta }-\displaystyle \frac{1}{{\lambda }_{0}{\lambda }_{\theta }}\right)+\left(\displaystyle \frac{-\mathrm{1/2}+({\chi }_{0}-{\chi }_{1})}{{{\lambda }_{0}}^{4}{{\lambda }_{\theta }}^{3}{\lambda }_{{\rm{r}}}}\right.\\ &&\quad \quad \quad \quad \;+\;\left.\displaystyle \frac{-\mathrm{1/3}\;+\;2{\chi }_{1}}{{{\lambda }_{0}}^{7}{{\lambda }_{\theta }}^{5}{{\lambda }_{{\rm{r}}}}^{3}}-\displaystyle \frac{1}{{{\lambda }_{0}}^{10}{{\lambda }_{\theta }}^{7}{{\lambda }_{{\rm{r}}}}^{3}}\right)\\ &&\quad \quad \quad \quad \;-\;({c}_{{H}^{+}}+{c}_{+}+{c}_{-}-{\bar{c}}_{{H}^{+}}-{\bar{c}}_{+}-{\bar{c}}_{-}){{\lambda }_{0}}^{2}{\lambda }_{{\rm{r}}}{\lambda }_{\theta }.\end{eqnarray} \tag{ 38 }$$

Also, the stretch components in terms of the current element radius are:

$$\begin{eqnarray}&&{\lambda }_{{\rm{r}}}=\displaystyle \frac{r(R)}{R},\;\;{\lambda }_{\theta }=\displaystyle \frac{{\rm{d}}}{{\rm{d}}R}r(R),\end{eqnarray} \tag{ 39 }$$

in which $r(R)$ is the radius of the element of the hydrogel in the current state. Substituting (38) and (39) in (37), a differential equation in terms of $r(R)$ and ${c}_{{H}^{+}}$ is obtained, and accompanying (33), a system of two nonlinear equations established. One method of dealing with this problem is elimination of ${c}_{{H}^{+}}$ from the two equations and solving the problem for the displacement field $r(R).$ The relevant boundary conditions for this problem are then expressed as:

$$\begin{eqnarray}&&r={R}_{{\rm{i}}}\;{\rm{at}}\;R={R}_{{\rm{i}}},\;\;\;{P}_{{\rm{r}}}=0\;{\rm{at}}\;R={R}_{{\rm{o}}}.\end{eqnarray} \tag{ 40 }$$

The material parameters are selected as $N\nu =0.04$ and $f=0.10.$ Also, the external solution molarity is $0.001\;M.$ To solve the present boundary value problem, a finite difference technique with Richardson extrapolation is used [44]. The calculated results as depicted in figure 6, are reported for different values of thickness ratio, namely ${R}_{{\rm{o}}}/{R}_{{\rm{i}}}.$ The normalized outer radius of the hydrogel shell is plotted versus the pH and temperature in figures 6(a) and (b), respectively. The normalized outer radius is defined as ratio of the current outer radius of the hydrogel shell $r({R}_{{\rm{o}}}),$ to its value at the reference point ${R}_{{\rm{o}}}.$ As shown in figure 6(a) increasing the pH results in an increase in the normalized outer radius. As the thickness ratio increases, the added material in the outer radius experiences a low level of confinement and as a result, can swell more freely with increase in pH. Therefore, increasing the thickness ratio, the behavior of the shell approaches to the free swelling state as understood from figure 6(a). The effect of temperature increase is also illustrated in figure 6(b) in which the hydrogel experiences shrinkage due to an increase in the temperature. In contrast to figures 6(a), in (b) the normalized outer radius is larger for smaller values of the thickness ratio. Again, the behavior of the shell approaches to the free swelling state as the thickness ratio increases.

Zoom In Zoom Out Reset image size

3.2.2. Numerical implementation

In the last section, the model is successfully used for solving an inhomogeneous swelling problem. Implementing the constitutive models of smart hydrogels as UHYPER in ABAQUS is popular in the literature because of the vast facilities of ABAQUS for FEM simulations [4, 43, 45]. Following other researchers, the proposed model is implemented in an ABAQUS user subroutine, namely UHYPER, to obtain a numerical tool for solving more complicated boundary value problems. To show the capability of the numerical scheme, inhomogeneous swelling of a hydrogel shell on a rigid pillar with boundary conditions identical to the analytical solution, presented in the previous section, is solved and the results are compared in figure 7 for ${R}_{{\rm{o}}}/{R}_{{\rm{i}}}=2.$ The material parameters and the external solution molarity are the same as those of the analytical solution presented in section 3.2.1. As shown in figure 7, the numerical and analytical results are in excellent agreement and the validity of the FEM solution is successfully checked.

Zoom In Zoom Out Reset image size

In addition, for ${R}_{{\rm{o}}}/{R}_{{\rm{i}}}=2$ we solve inhomogeneous swelling of a hydrogel shell on a rigid pillar which may experience a contact with rigid flat walls of a channel as shown in figure 5. The half-width of the channel, $D,$ is selected such that $D/{R}_{{\rm{i}}}=2.25.$ Due to the symmetry, simulations are carried out for a quarter of the real spherical shell and the results are shown in figure 8.

Zoom In Zoom Out Reset image size

In this figure the columns from left to right show the results at temperatures 290 K, 300 K and 310 K, respectively and also, the top-down rows are the results for pH values of 2, 5, 7 and 10, respectively. The material parameters, the stress-free condition and the external solution molarity are identical to the previous numerical problem. The inner radius of the hydrogel shell is constrained to the rigid pillar while the outer radius of the shell experiences a contact with the walls. Thus, this problem contains a complicated coupling between the temperature changes, pH changes, large deformations and a mechanical contact. As shown in figure 8, a temperature decrease or a pH increase may result in contact between the shell and the flat walls. For example when pH = 10, the shell is initially in contact with the wall at T = 290 K. Any increase in the temperature or decrease in pH lead to shrinkage of the shell and as a result they are not in contact anymore as observed at T = 300 K and T = 310 K. For the case that no contact occurs at the low temperature, an increase in the channel opening eventuates from increasing in the temperature or decreasing in pH. These observations are reasonable and the validated numerical scheme derived in this section can be employed to simulate boundary value problems for the swelling of the pH/temperature sensitive hydrogels.

In this work the equilibrium behavior of the pH/temperature sensitive PNIPAM hydrogel was studied and a model was presented with a good agreement with the experimental results available in the literature. The free and biaxial constrained swelling of the pH/temperature sensitive hydrogel was investigated and some parameter studies were performed. One of the main advantageous of the presented model was its continuity especially against the temperature variations which makes it appropriate for implementing in the inhomogeneous swelling problems whether analytically or numerically. Thereafter, the inhomogeneous swelling of a spherical shell on a hard core was analytically studied and the results were presented for various thicknesses of the shell. Also, the model was implemented in ABAQUS via a user subroutine UHYPER. Then, the inhomogeneous swelling of a spherical shell on a hard core, with boundary conditions identical to the analytical one, was solved and the validity of the numerical results was successfully confirmed. Finally, the numerical simulation of the complicated problem of contact between the hydrogel shell and a channel with rigid walls due to the pH or temperature variations was performed with reasonable results. Based on the results, it is concluded that the presented model is suitable for predicting the behavior of pH/temperature hydrogels and has the capability for implementing in the homogeneous and inhomogeneous problems.

The second author is grateful for the research support of the Iran National Science Foundation (INSF).