2.1. Free energy density of the hydrogel

We consider a hydrogel in a quiescent bath of fluid that has absorbed fluid and deformed. Both of these processes have energy costs associated with them: a chemical contribution from the mixing of fluid molecules and polymer chains, and a physical contribution from the elastic deformation of the polymer. The nominal free energy density of the hydrogel in such a state can then be written as the sum of these two contributions (Flory & Rehner Reference Flory and Rehner1943a,Reference Flory and Rehnerb)

(2.1)\begin{equation} \tilde{W} = \tilde{W}_{{elast}} + \tilde{W}_{{mix}}, \end{equation}

where $\tilde {W}_{{elast}}$ represents the elastic energy density and $\tilde {W}_{{mix}}$ is the mixing energy density. Here, $\tilde {W}$ is the Helmholtz free energy of the mixture per unit volume of the undeformed dry polymer. Throughout this work, we use tildes, $\tilde {\cdot }$, to denote dimensional quantities.

We consider a gel deformation from a reference state, resulting in a local factor increase in volume of the body, $J$, relative to this reference state (at which $J=1$). The principal stretches $\lambda _i$ for $i\in \{1,2,3\}$ are the factor increase in lengths from the reference state in each direction; they are the major and minor axes of an ellipsoid that results from the deformation (locally) of a unit sphere in the reference state. The local volume change $J$ is therefore the product of the stretches, $J=\lambda _1\lambda _2\lambda _3$.

Such a deformation has an energetic cost to stretching (or compressing) the polymer chains. We consider a neo-Hookean model for the elastic energy density:

(2.2)\begin{equation} \tilde{W}_{{elast}} = \frac{\tilde{G}}{2} \left[ \lambda_1^2 + \lambda_2^2 + \lambda_3^2 - 3 - 2 \log J \right], \end{equation}

where $\tilde {G}$ is the shear modulus of the gel. Commonly, the shear modulus is found to follow a linear behaviour with temperature, $\tilde {G}=\tilde {k}_B\tilde {T}/\tilde {\varOmega }_p$, where $\tilde {T}$ is the temperature (in kelvin), $\tilde {k}_B$ is the Boltzmann constant, and $\tilde {\varOmega }_p$ is the volume of polymer per polymer molecule in the reference state. (Note that often this is written as $\tilde {G}=\tilde {N}\tilde {k}_B\tilde {T}$, with $\tilde {N}$ the number density of polymer chains in the reference state, but $\tilde {N}$ and $\tilde {\varOmega }_p$ are related simply through $\tilde {N}=1/\tilde {\varOmega }_p$.)

For the mixing contribution to the free energy density, we use the Flory–Huggins polymer theory (Flory Reference Flory1942; Huggins Reference Huggins1942; Flory & Rehner Reference Flory and Rehner1943b) to get the mixing energy density

(2.3)\begin{equation} \tilde{W}_{{mix}}= \frac{\tilde{k}_B \tilde{T}}{\tilde{\varOmega}_f}\,J \left[ \phi\log\phi + \chi\phi(1-\phi) \right]. \end{equation}

Here, $\phi$ is the porosity (i.e. the local proportion of fluid per unit mixture volume), and $\tilde {\varOmega }_f$ is the volume of fluid per fluid molecule in the unmixed state. Note that in (2.3), there is a factor $J$ to account for the local volume change relative to the reference state, due to deformation of the gel.

The term $\chi$ in (2.3) is the (Flory) mixing parameter, and gives the degree to which the polymer and fluid like to mix; in general, this mixing parameter can depend on both temperature and mixture composition, $\chi = \chi (\tilde {T},\phi )$. It has been suggested that both components are required to reproduce experimental data accurately (Lopez-Leon & Fernandez-Nieves Reference Lopez-Leon and Fernandez-Nieves2007). There are many possible options for this function, but for simplicity, we focus on a relationship that is linear in both the temperature and the polymer fraction:

(2.4)\begin{equation} \left.\begin{gathered} \chi(\tilde{T},\phi) = \chi_0(\tilde{T}) + \chi_1(\tilde{T})\,(1-\phi), \\ \chi_0(\tilde{T})= A_0 + B_0 \tilde{T},\quad \chi_1(\tilde{T}) = A_1 + B_1 \tilde{T}, \end{gathered}\right\} \end{equation}

where $A,B$ are constants, as modelled by Cai & Suo (Reference Cai and Suo2011). We note, however, that many other models for this mixing parameter exist: constant (Hong et al. Reference Hong, Zhao, Zhou and Suo2008; Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016; Hennessy et al. Reference Hennessy, Münch and Wagner2020), dependent only on the temperature (Tomari & Doi Reference Tomari and Doi1995; Chester & Anand Reference Chester and Anand2011), or with more complicated behaviours, such as including terms inversely proportional to temperature (Tanaka Reference Tanaka1978; Quesada-Pérez et al. Reference Quesada-Pérez, Maroto-Centeno, Forcada and Hidalgo-Alvarez2011) or quadratic powers (Drozdov & Christiansen Reference Drozdov and Christiansen2017). We have chosen to use the linear relation (2.4) because it is simple, yet still includes the important local composition dependence. Despite this simplicity, different choices of the constants can give a range of different behaviours, both quantitatively and qualitatively.

To relate the polymer fraction and $J$, we must define the reference state, where $J=1$ everywhere. For simplicity, and in line with many other works on hydrogel swelling (Doi Reference Doi2009; Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016; Hennessy et al. Reference Hennessy, Münch and Wagner2020), we consider the reference state to be the dry state, where the polymer is a solid block with $\phi =0$; we note, however, that there is some debate as to what an appropriate reference state is, and whether such a completely dry state is indeed physical (Quesada-Pérez et al. Reference Quesada-Pérez, Maroto-Centeno, Forcada and Hidalgo-Alvarez2011). Treating the polymer chains as incompressible and the mixture as ideal, the relation between the volume change and solid fraction is then given simply by

(2.5)\begin{equation} J=\frac{1}{1-\phi}. \end{equation}

2.2. Free-swelling equilibria

We now turn to consider the equilibrium swelling of a hydrogel. Under free-swelling conditions, where no external stresses are applied to the hydrogel, and the background fluid in the surrounding environment is static, the stretches in each direction must be equal, $\lambda _i=\lambda$ with $J=\lambda ^3$. The free energy density can therefore be rewritten as a function of $\lambda$ only:

(2.6)\begin{align} \tilde{W}(\lambda) = \frac{3\tilde{G}}{2} \left[ \lambda^2 - 1 - 2 \log \lambda \right] + \frac{\tilde{k}_B \tilde{T}}{\tilde{\varOmega}_f} \left[ (\lambda^3-1)\log\left(1-\frac{1}{\lambda^3}\right) + \chi\left(1-\frac{1}{\lambda^3}\right) \right], \end{align}

with $\chi =\chi _0(\tilde {T}) + \chi _1(\tilde {T})/\lambda ^3$.

Equilibria are found by minimising the free energy density so that $\mathrm {d} \tilde {W}(\lambda )/\mathrm {d} \lambda =0$; the equilibrium stretch $\lambda$ must therefore satisfy the equation

(2.7)\begin{equation} \frac{\lambda^2 - 1}{\lambda^3} + \varOmega \left[ \log\left(1-\frac{1}{\lambda^3}\right) + \frac{1}{\lambda^3} + \frac{\chi_0-\chi_1}{\lambda^6} + \frac{2\chi_1}{\lambda^9} \right] = 0. \end{equation}

Given the mixing parameter $\chi (\tilde {T},\phi )$, the equilibrium solutions $\lambda =\lambda (\tilde {T};\varOmega )$ can then be calculated as a function of temperature with a single parameter, $\varOmega = \tilde {\varOmega }_p/\tilde {\varOmega }_f$. This parameter encodes the relative importance of the mixing and elastic energies, but it is also the ratio of volume per polymer chain compared to a fluid molecule (or, alternatively, a ratio of their densities) and is in general expected to be large. The value of this parameter may vary depending on factors such as the gel composition and curing conditions (since these affect the stiffness of the gel), and can take values from tens (Quesada-Pérez et al. Reference Quesada-Pérez, Maroto-Centeno, Forcada and Hidalgo-Alvarez2011; Drozdov Reference Drozdov2014) to hundreds (Tanaka Reference Tanaka1978; Cai & Suo Reference Cai and Suo2011) and above (Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016).

To complete the calculation of the equilibrium swelling, it therefore remains to define the functions $\chi _0$ and $\chi _1$, i.e. choose values for the constants $A$ and $B$ in (2.4). Even after restricting the equilibrium swelling to certain qualitative behaviour (such as being more swollen at low temperatures, with a swelling transition within a certain temperature range), there remains a wide range of possible choices for these parameters; we consider two possibilities using equilibrium data extracted from the literature.

In the first case, we take the values used previously by Cai & Suo (Reference Cai and Suo2011) in modelling the mechanics of a thermo-responsive hydrogel:

(2.8)\begin{equation} \left.\begin{gathered} \varOmega = 100, \\ A_0 ={-}12.947,\quad A_1 = 17.92, \\ B_0 = 0.04496\,\text{K}^{{-}1},\quad B_1 ={-}0.0569\,\text{K}^{{-}1}. \end{gathered}\right\} \end{equation}

The values of $A$ and $B$ given here were determined by Afroze, Nies & Berghmans (Reference Afroze, Nies and Berghmans2000) from a prepared sample of PNIPAM. In this case, the equilibrium curve is multi-valued, taking a characteristic S-shape as shown in figure 3(a). We will refer to this as the Afroze–Nies–Berghmans (ANB) solution. Small changes in temperature can result in a large (discontinuous) jump in the level of equilibrium swelling. This equilibrium behaviour has been observed experimentally in some thermo-responsive hydrogels (e.g. Tanaka Reference Tanaka1978; Hirokawa & Tanaka Reference Hirokawa and Tanaka1984; Sato Matsuo & Tanaka Reference Sato Matsuo and Tanaka1988). The volume phase transition temperature for shrinking in this example is $\tilde {T}_c\approx 305.8$ K, and we expect a large degree of shrinking when the temperature is increased above this. However, we note that there is a small range of temperatures over which there are multiple possible equilibrium solutions (the system exhibits hysteresis); when decreasing the temperature from an originally deswollen state, the hydrogel will swell significantly only once below the volume phase transition for swelling, $\tilde {T}_c\approx 304.7$ K, which is lower than that for the shrinking transition.

Figure 3. Equilibrium swelling of thermo-responsive hydrogels under free-swelling conditions (solid curves) for (a) the ANB parameters (2.8), and (b) the HHT parameters (2.9). The stretch (or factor increase in lengths compared to the dry state) $\lambda$ is calculated as a function of temperature $\tilde {T}$ using (2.7). The dashed lines show the spinodal curves, where $\partial ^2 \tilde {W}/\partial \lambda _i^2=0$, with the spinodal regions being the darker shaded areas, $\partial ^2 \tilde {W}/\partial \lambda _i^2<0$. The lighter shaded regions are where the isotropically-stretched state can coexist alongside another with a discontinuous normal stretch across the interface. In (b), the fitted curve is compared to (rescaled) data from Hirotsu et al. (Reference Hirotsu, Hirokawa and Tanaka1987), and the inset shows a zoom of the region close to the swelling transition.

We also consider an example that is found by fitting (2.7) to data captured from figure 3 of Hirotsu, Hirokawa & Tanaka (Reference Hirotsu, Hirokawa and Tanaka1987) for another sample of PNIPAM. To do this, we rearrange (2.7) to $\tilde {T}=\tilde {T}(\lambda )$, and fit the appropriate six parameters, which are $\varOmega,A_0,B_0,A_1,B_1$ and the scaling between their gel radius and the stretch (i.e. determine the radius at which $\lambda =1$, since this is unknown in the experiment). This fitting was implemented using the MATLAB least-squares fitting nonlinlsq from the Optimization Toolbox, running the fitting a large number of times to remove any local minima that were found. The best-fit parameters are found to be

(2.9)\begin{equation} \left.\begin{gathered} \varOmega = 720, \\ A_0 ={-}62.22,\quad A_1 ={-}58.28, \\ B_0 = 0.20470\,\text{K}^{{-}1},\quad B_1 = 0.19044\,\text{K}^{{-}1}, \end{gathered}\right\} \end{equation}

with the scaling between the data and the equilibrium stretch found to be 0.4127. These parameters are noticeably different from those in (2.8), with each changing by at least several fold between the two examples, but this is not particularly unexpected considering $\varOmega$ has a wide range of reported values in the literature, depending on the choice of monomer and cross-linker and how the hydrogel is prepared.

The equilibrium curve that is calculated from these parameters is shown in figure 3(b), with the comparison between the fitting and the data also illustrated. There is a sharp transition in the equilibrium degree of swelling close to $\tilde {T}_c\approx 307.6$ K (when both swelling and shrinking). The fitted equilibrium solutions are again multi-valued around the transition, but over a much smaller region (approximately $0.01$ K wide) that would certainly not be detectable experimentally, as can be seen in the inset to figure 3(b). This solution and parameter set will be referred to as the Hirotsu–Hirokawa–Tanaka (HHT) solution, to distinguish it from the ANB case (2.8).

2.3. Coexistence and spinodal curves

In addition to the equilibrium curves in figure 3, we have also plotted two shaded regions for each model hydrogel. To understand these, we first note that the free energy density can be written as a function of the principal stretches only, $\tilde {W}=\tilde {W}(\lambda _1,\lambda _2,\lambda _3)$, since $J=\lambda _1\lambda _2\lambda _3=1/(1-\phi )$. The darker shaded region is bounded by the dashed curve

(2.10)\begin{equation} \frac{\partial^2\tilde{W}}{\partial \lambda_i^2} (\lambda_1=\lambda_2=\lambda_3=\lambda; \tilde{T}) = 0. \end{equation}

Inside the darker shaded region, a uniformly swollen hydrogel is linearly unstable to uniaxial perturbations. We call this region the spinodal region, and the curve the spinodal curve.

We also consider when an isotropically swollen state, with stretch $\lambda$, can be in local equilibrium with an adjoining gel state with a different degree of swelling, with a sharp interface between them. The two stretches in directions tangential to the interface must vary continuously across the interface, $\lambda _t=\lambda$, but the stretch in the normal direction can be discontinuous, $\lambda _n=\lambda _{{new}}\neq \lambda$, and can be either more or less swollen than the isotropic state. At the interface, the following balance condition must be satisfied (Sekimoto & Kawasaki Reference Sekimoto and Kawasaki1989; Sekimoto Reference Sekimoto1993; Tomari & Doi Reference Tomari and Doi1995):

(2.11)\begin{equation} \frac{\partial \tilde{W}}{\partial \lambda_n} (\lambda,\lambda;\tilde{T}) = \frac{\partial \tilde{W}}{\partial \lambda_n} (\lambda_{{new}}, \lambda;\tilde{T}) = \frac{\tilde{W}(\lambda_{{new}},\lambda;\tilde{T})-\tilde{W} (\lambda,\lambda;\tilde{T})}{\lambda_{{new}}-\lambda}, \end{equation}

where, for simplicity, we here write the arguments of the free energy density as $\tilde {W}=\tilde {W}(\lambda _n,\lambda _t;\tilde {T})$. As we will see later when defining the stress, the first equality here is enforcing a mechanical stress balance across the interface so that the normal stress is continuous. The second equality ensures that the two sides are in chemical equilibrium, so the gradient $\partial \tilde {W}/\partial \lambda _n$ is well-defined at the interface. At a given temperature $\tilde {T}$, (2.11) therefore gives two equations in two unknowns, $\lambda$ and $\lambda _{{new}}$, which can be solved and plotted in $(\tilde {T},\lambda )$-space; we instead choose values of $\lambda$ and solve for both $\tilde {T}$ and $\lambda _{{new}}$. Since the free energy density is symmetric in the principal stretches (i.e. $\lambda _1,\lambda _2,\lambda _3$), we can solve (2.11) by taking derivatives with respect to any principal stretch as the normal direction, so that, for example, $\lambda _n=\lambda _1$ without loss of generality, and plot the solution as the dotted line in figure 3, called the coexistence curve. In the lighter shaded region, bounded by this curve, the described coexistence between the neighbouring isotropic and anisotropic states can still occur mechanically, but the latter is energetically favourable and will be expected to invade the domain.

Considering a hydrogel that is initially isotropically swollen with stretch $\lambda$ at temperature $\tilde {T}$, we may expect to observe different dynamic behaviours depending upon which region of $(\tilde {T},\lambda )$-space the hydrogel starts in. In the spinodal region, the initial state is unstable, so we might see spontaneous swelling or shrinkage within the interior of the hydrogel due to the growth of any small perturbations – spinodal decomposition. If, instead, the initial state is within the coexistence region, then there is an energetically favourable alternative state. However, the gel must be stimulated at some specific location to overcome the energy barrier required to generate this alternative state, which can then invade the domain – this is nucleation and growth. Therefore, when starting in these two regions, we expect to see phase separation in the hydrogel. The behaviours of gels in the non-shaded regions of figure 3 are slightly more complicated to predict, since as a gel swells or shrinks, it is generally no longer isotropically swollen, thus phase separation could occur after a delay, or not at all. We will turn to investigate some of these behaviours in a spherical thermo-responsive hydrogel, but first we introduce a model for the dynamics of such a gel.

3.1. Strains

For a sphere undergoing a spherically symmetric deformation with a radial displacement $\tilde {u}(\tilde {r},\tilde {t})$, the principal stretches are (Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016)

(3.1a,b)\begin{equation} \lambda_{r} = \left( 1 - \frac{\partial \tilde{u}}{\partial \tilde{r}} \right)^{{-}1},\quad \lambda_\theta = \lambda_\varphi = \left( 1 - \frac{\tilde{u}}{\tilde{r}} \right)^{{-}1}, \end{equation}

where $\tilde {r}$ is the radial coordinate from the centre of the sphere, $\theta$ and $\varphi$ are the spherical angle coordinates, and $\tilde {t}$ is time. This can be seen by considering a small piece of the gel that is initially at radial distance $\tilde {r}_0$ in the undeformed reference state; once deformed, it is stretched locally by an amount $\partial \tilde {r}/\partial \tilde {r}_0$ in the radial direction and $\tilde {r}/\tilde {r}_0$ in the angular directions, and we have that $\tilde {u}=\tilde {r}-\tilde {r}_0$.

We can then relate the stretches to the porosity through the local volume change $J$, since we have $J = \lambda _r \lambda _\theta ^2$, but also (2.5) holds (since the reference volume is taken to be the dry state with $\phi =0$). Combining these relations for $J$, we find a first-order differential equation in $\tilde {r}$ for the displacement $\tilde {u}$. This can be integrated once to determine the displacement in terms of the porosity:

(3.2)\begin{equation} \tilde{u}(\tilde{r},\tilde{t}) = \tilde{r} - \left(\tilde{r}^3-3 \int_{0}^{\tilde{r}} x^2\, \phi(x,\tilde{t})\,\mathrm{d} x \right)^{1/3}, \end{equation}

with $\tilde {u}(0,\tilde {t})=0$ at the centre of the sphere, and $\tilde {u}(\tilde {a},\tilde {t})=\tilde {a}-\tilde {a}_d$ at the sphere edge.

Therefore, given the porosity $\phi (\tilde {r},\tilde {t})$ in the gel sphere, we can calculate the radial displacements $\tilde {u}(\tilde {r},\tilde {t})$ and hence the stretches $\lambda _i(\tilde {r},\tilde {t})$.

3.2. Stresses

Having determined the strains in the sphere, it remains to calculate the stresses. The mechanical stresses $\tilde {\sigma }$ acting on the gel can be written in terms of gradients of the free energy density $\tilde {W}$, given in (2.1)–(2.3). We decompose the stress components into contributions from the elastic deformation and the chemical mixing:

(3.3)\begin{equation} \tilde{\sigma}_i = \tilde{\sigma}_i' - \tilde{p}, \end{equation}

where $\tilde {\sigma }_i$ are the principal Cauchy stresses within the mixture. The Terzaghi effective stress (which is the elastic contribution) is defined by

(3.4)\begin{equation} \tilde{\sigma}_i' = \frac{\lambda_i}{J}\,\frac{\partial \tilde{W}_{{elast}}}{\partial \lambda_i} = G\, \frac{\lambda_i^2-1}{J}. \end{equation}

The pressure is split into an osmotic pressure $\tilde {\varPi }$ due to the interaction of the fluid and polymer, and the chemical potential $\tilde {\mu }$:

(3.5)\begin{equation} \tilde{p} = \tilde{\varPi} + \frac{\tilde{\mu}}{\tilde{\varOmega}_f}, \end{equation}

where the osmotic pressure is given by

(3.6)\begin{equation} \tilde{\varPi} ={-} \frac{\mathrm{d} \tilde{W}_{{mix}}}{\mathrm{d} J} ={-}\frac{\tilde{k}_B \tilde{T}}{\tilde{\varOmega}_f} \left[\log\left(1-\frac{1}{J}\right) + \frac{1}{J} + \frac{\chi_0-\chi_1}{J^2} + \frac{2\chi_1}{J^3} \right]. \end{equation}

The justification for the elastic stress and osmotic pressure being written in terms of derivatives of the free energy density can be found by considering the work done when changing the hydrogel mixture by a small volume locally (see (4) and (5) of Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016). More generally, the stress tensor can be written in terms of derivatives of the free energy density with respect to the deformation tensor (e.g. Doi Reference Doi2009; Tadmor, Miller & Elliott Reference Tadmor, Miller and Elliott2012), but we stick to the simpler spherically symmetric model outlined above.

The term $\tilde {\mu }/\tilde {\varOmega }_f$ in (3.5) arises to account for the work required to move fluid into the mixture. An alternative view of this term is that it is the pervadic pressure of the fluid in the mixture – a connected fluid-only bath in local equilibrium with the fluid in the mixture would be at this pressure (Peppin, Elliott & Worster Reference Peppin, Elliott and Worster2005; Etzold et al. Reference Etzold, Linden and Worster2021).

Note that the precise value of the chemical potential does not affect the equilibrium calculations, since it must be uniform and equal to the ambient value in the surrounding fluid. As we will see, only gradients in the chemical potential alter the dynamics. In reality, the value of the ambient chemical potential should differ at each temperature, but this will turn out to be unimportant for our simulations as we consider only dynamics at a fixed temperature. More care may be needed in more complex scenarios, but note that adding a uniform amount to the chemical potential results in only a uniform increase in the stress, with otherwise no effect on our results.

3.3. Poro-elastic flow

To resolve the dynamic behaviour of the swelling and shrinking hydrogel sphere, we need to understand the fluid flow within the gel pore space. The fluid velocity $\tilde {v}_f$, relative to the motion of the gel $\tilde {v}_p$, is driven by the chemical potential gradient (Darcy flow):

(3.7)\begin{equation} \phi(\tilde{v}_f-\tilde{v}_p) ={-} \frac{\tilde{k}(\phi)}{\tilde{\eta}}\,\frac{\partial }{\partial \tilde{r}} \left(\frac{\tilde{\mu}}{\tilde{\varOmega}_f} \right), \end{equation}

where $\tilde {\eta }$ is the fluid viscosity and $\tilde {k}$ is the permeability of the hydrogel. Note that the Darcy flow here is driven only by the chemical potential because, as mentioned previously, it is equivalent to the pressure of the fluid within the mixture by itself (Peppin et al. Reference Peppin, Elliott and Worster2005). This Darcy flow is what makes the dynamic modelling different from that of Tomari & Doi (Reference Tomari and Doi1995), since here we account explicitly for the interstitial fluid flow relative to the solid matrix, rather than accounting for it through an energy dissipation via friction.

We take the permeability to be isotropic and a function of the porosity $\tilde {k}(\phi )$. For simplicity, we will present results with a constant permeability $\tilde {k}_0$, but keep a general form in our equations. We discuss an alternative permeability model used by Bertrand et al. (Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016) in Appendix C.

As the hydrogel swells or shrinks, any deformation of the solid structures must be accompanied by a corresponding displacement of fluid; in this manner, conservation of solid volume dictates that the porosity must obey

(3.8)\begin{equation} \frac{\partial }{\partial \tilde{t}}(1-\phi) + \frac{1}{\tilde{r}^2}\, \frac{\partial }{\partial \tilde{r}} [\tilde{r}^2 (1-\phi)\tilde{v}_p] = 0, \end{equation}

and a similar equation holds for the fluid velocity $\tilde {v}_f$, with $\phi$ replacing $1-\phi$. Combining these two relations gives a simple relation of no net flux in any cross-section:

(3.9)\begin{equation} \phi \tilde{v}_f + (1-\phi) \tilde{v}_p = 0. \end{equation}

Combining (3.7)–(3.9), we find a partial differential equation for the evolution of the porosity in terms of the gradient of the chemical potential:

(3.10)\begin{equation} \frac{\partial \phi}{\partial \tilde{t}} = \frac{1}{\tilde{r}^2}\,\frac{\partial }{\partial \tilde{r}} \left[\tilde{r}^2(1-\phi)\,\frac{\tilde{k}(\phi)}{\tilde{\eta} \tilde{\varOmega}_f}\, \frac{\partial \tilde{\mu}}{\partial \tilde{r}} \right]. \end{equation}

Note that this is a Reynolds equation of the form $\partial \phi /\partial \tilde {t} + \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {Q} = 0$, with radial flux $Q = -[(1-\phi )\,\tilde {k}(\phi ) /(\tilde {\eta }\tilde {\varOmega }_f)] \times \partial \tilde {\mu }/\partial \tilde {r}$. Conservation of fluid volume then means that the sphere radius changes due to the outward radial flux at the sphere edge:

(3.11)\begin{equation} \frac{\mathrm{d} \tilde{a}}{\mathrm{d} \tilde{t}} = \frac{\tilde{k}(\phi)}{\tilde{\eta}\tilde{\varOmega}_f}\, \frac{\partial \tilde{\mu}}{\partial \tilde{r}}\quad \text{at }\tilde{r}=\tilde{a}. \end{equation}

A stress balance in the solid enforces $\boldsymbol {\nabla }\boldsymbol {\cdot }\tilde {\sigma }=0$; considering the radial component of this, we can determine the chemical potential gradient in terms of the Terzaghi stresses and osmotic pressure via

(3.12)\begin{equation} \tilde{\varOmega}_f^{{-}1}\,\frac{\partial \tilde{\mu}}{\partial \tilde{r}} = \frac{\partial \tilde{\sigma}_r'}{\partial \tilde{r}} + 2\, \frac{\tilde{\sigma}_r'-\tilde{\sigma}_\theta'}{\tilde{r}} - \frac{\partial \tilde{\varPi}}{\partial \tilde{r}}. \end{equation}

3.4. Non-dimensionalisation

We rescale the governing equations using the dry sphere radius $\tilde {a}_d$, the permeability scale $\tilde {k}_0$, the stress scale $\tilde {k}_B \tilde {T}_{{end}}/\tilde {\varOmega }_p$ (for a final temperature $\tilde {T}_{{end}}$), the chemical potential scale $\tilde {k}_B \tilde {T}_{{end}}\tilde {\varOmega }_f/\tilde {\varOmega }_p$, and the time scale

(3.13)\begin{equation} \tilde{\tau} = \frac{\tilde{\eta} \tilde{\varOmega}_p \tilde{a}_d^2}{\tilde{k}_0\tilde{k}_B \tilde{T}_{{end}}}. \end{equation}

Note that the time scale here is proportional to the square of the size of the gel sphere, $\tilde {\tau }\propto \tilde {a}_d^2$, just as was found for the swelling of gels by Tanaka & Fillmore (Reference Tanaka and Fillmore1979); here, the effective diffusivity is $\tilde {D} \sim (\tilde {k}_0\tilde {k}_B \tilde {T}_{{end}}) / (\tilde {\eta } \tilde {\varOmega }_p)$. In general, we expect this time scale to be much longer than that for thermal diffusion: Tanaka & Fillmore (Reference Tanaka and Fillmore1979) reported a diffusivity of $3\times 10^{-7}\,\textrm {cm}^2\,\textrm {s}^{-1}$, which is $10^4$ times smaller than a typical thermal diffusivity for water, suggesting that typically, swelling takes $100$ times longer than the diffusion of heat. We obtain a similar effective diffusivity for our swelling gel model when $\tilde {\eta }=10^{-3}\,\textrm {Pa}\,\textrm {s}$, $\tilde {\varOmega }_p=3\times 10^{-25}\,\textrm {m}^3$ and $\tilde {k}_0=8\times 10^{-20}\,\textrm {m}^3$ (Bertrand et al. Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016). We therefore model temperature changes as being effectively instantaneous and uniform over the whole system.

We also note that, in general, the fluid viscosity $\tilde {\eta }$ may be temperature-dependent. Since this viscosity appears in the dimensionless model only through the time scale $\tilde {\tau }$, defined in (3.13), and all simulations are to be conducted at fixed temperature, this does not affect any of our dimensionless results. However, there will be a quantitative effect that must be accounted for when considering the real dimensional dynamic swelling and shrinking times.

From this point onwards, we write all equations and results in dimensionless terms, removing the tildes used in earlier equations that denote the dimensional forms of variables. The only remaining dimensional quantity is the temperature $\tilde {T}$, since its dimensional value matters, although in all of our equations it will be non-dimensionalised by multiplying by $B_0$ or $B_1$ in $\chi$.

The dimensionless partial differential equation for the evolution of the porosity is

(3.14)\begin{equation} \frac{\partial \phi}{\partial t} = \frac{1}{r^2}\,\frac{\partial }{\partial r} \left[r^2(1-\phi)\,k(\phi)\,\frac{\partial \mu}{\partial r} \right], \end{equation}

for a dimensionless permeability $k(\phi )$ (which we take to be $k=1$). The chemical potential gradient is given in terms of the stresses:

(3.15)\begin{equation} \frac{\partial \mu}{\partial r} = \frac{\partial \sigma'_r}{\partial r} + 2\, \frac{\sigma'_r-\sigma'_\theta}{r} - \frac{\partial \varPi}{\partial r}. \end{equation}

These stresses are defined in terms of the strains through

(3.16)\begin{equation} \sigma'_i = \frac{\lambda_i^2-1}{J}, \end{equation}

and similarly for the osmotic pressure

(3.17)\begin{equation} \varPi ={-} \varOmega \left[ \log\left(1-\frac{1}{J}\right) + \frac{1}{J} + \frac{\chi_0-\chi_1}{J^2} + \frac{2\chi_1}{J^3} \right], \end{equation}

where $\varOmega =\tilde {\varOmega }_p/\tilde {\varOmega }_f$ is the parameter mentioned previously that relates the relative importance of the mixing and elastic energies.

The strains can be calculated directly from the porosity $\phi (r,t)$ by first determining the radial displacement

(3.18)\begin{equation} u = r - \left( r^3 - 3 \int_{0}^{r} x^2 \phi\,\mathrm{d} x \right)^{1/3}. \end{equation}

The stretches and volume change are then

(3.19a–c)\begin{equation} J = \frac{1}{1-\phi}, \quad \lambda_\theta = \frac{1}{1-u/r}, \quad \lambda_r = \frac{J}{\lambda_\theta^2}. \end{equation}

The sphere edge moves according to the chemical potential gradient there:

(3.20)\begin{equation} \frac{\mathrm{d} a}{\mathrm{d} t} = k(\phi)\,\frac{\partial \mu}{\partial r} \quad \text{at }r=a. \end{equation}

Note that at the sphere edge, the displacement is $u=a-1$ and $\lambda _\theta =a$. We also enforce chemical equilibrium, $\mu =\mu _0$, and a radial stress balance, $\sigma _r=-\mu _0$, at the sphere boundary, where $\mu _0$ is the chemical potential of the surrounding background fluid. Note that the radial stress just inside the sphere is not zero because it must balance the external fluid pressure (cf. pervadic pressure). Combining these two conditions enforces $\sigma '_r=\varPi$, which gives a nonlinear equation for $J$ at the sphere edge:

(3.21)\begin{equation} \frac{J}{a^4} - \frac{1}{J} + \varOmega \left[\log\left(1-\frac{1}{J}\right) + \frac{1}{J} + \frac{\chi_0-\chi_1}{J^2} + \frac{2 \chi_1}{J^3} \right] = 0. \end{equation}

Note that these boundary conditions mean that the hydrogel dynamics has no dependence on the external chemical potential $\mu _0$. This differs from the model of Bertrand et al. (Reference Bertrand, Peixinho, Mukhopadhyay and MacMinn2016), where there is a stress discontinuity at the sphere edge, but does align with other models (e.g. Hennessy et al. Reference Hennessy, Münch and Wagner2020).

Further, we are neglecting any mechanical feedback from the subsequent external fluid flow. We expect the outside flow to be unimportant to the swelling dynamics because the polymer's incompressibility means that any fluid flux through the boundary must be compensated for exactly by the boundary motion itself (except for, perhaps, some inhomogeneity at the pore scale). Some recent studies have shown that responsive hydrogels are able to generate substantial bulk fluid motions, such as cyclic swelling of bilayer ribbons causing net swimming (Tanasijević et al. Reference Tanasijević, Jung, Koens, Mourran and Lauga2022), but in this reported case, this was possible because of some significant asymmetry in the material composition of the particle due to impermeable sections of the boundary, which caused shape change and asymmetric fluid flow across the boundary.

In summary, our model comprises a neo-Hookean free energy density for deformation (2.2), the Flory–Huggins free energy density of mixing for a non-ionic gel (2.3), a particular form for the temperature and composition dependence of the mixing parameter (2.4), and poro-elastic dynamics with fluid motion given by Darcy flow (3.7). This combination results in the dimensionless governing equations for the evolution of the swelling or shrinking thermo-responsive gel, given in (3.14)–(3.21).

3.5. Numerical scheme

We solve these equations numerically using a spatial grid in the range $N=200{-}400$ points and time step ${\rm \Delta} t = 10^{-7}{-}10^{-8}$, with results recorded every $t=10^{-4}$. Such a small time step was used to maintain numerical stability with high spatial resolution for this simple numerical scheme. Details of the scheme are given in Appendix A.

4.3. Comparison of results

In both model systems, we have seen a range of different dynamic behaviours. Both cases exhibited swelling that progressed smoothly inwards, and shrinking that resulted in a sharp travelling front between a swollen core and shrunken shell. In each model gel, we were also able to find an example of the opposite behaviour.

The occurrence of front formation during shrinking seems to align with whether the initial homogeneous state is within the bounds of the coexistence regions that are illustrated in figure 3. When starting there, the gel initially phase separates close to the boundary, forming two regions with different degrees of swelling; the newly created shrunken state then spreads inwards, replacing the swollen core.

Similar behaviour for shrinking gels was also found by Tomari & Doi (Reference Tomari and Doi1995). They showed that (for their model hydrogel) there are equilibria that appear linearly stable but exist within the coexistence region and so are thermodynamically unstable, and calculated the intersection between the coexistence and equilibrium curves. Their dynamic simulations revealed an inwards-propagating front, similar to ours, but with an an overshoot in the amount of shrinking behind (outside) the front that we do not observe. It is not clear whether this difference is a numerical artefact or related to the change in the dynamic model to incorporate Darcy flow.

We also found that phase separation, although apparently common in the phase space of our shrinking hydrogel examples, does not always occur when heating a thermo-responsive hydrogel. For small temperature changes above the transition, one of our model hydrogels (see figure 10) shrank smoothly with no front formation.

Swelling of the hydrogel generally occurs smoothly towards the ultimate equilibrium state, since the dynamics is often not affected by the phase separation mechanisms seen in the shrinking dynamics; the coexistence and spinodal regions are mostly at higher temperatures. However, we did find that phase separation could occur when swelling in some cases, with the formation of a small but sharp front in one of our examples (figure 7). Here, the initial condition was not in the coexistence region, so no front immediately formed, but as the solution evolved, some (now anisotropic) part of the hydrogel passed its coexistence limit and a small inwards-travelling front formed between a swollen shell and a shrunken core. (Recall that the coexistence curve given by (2.11) was calculated for an isotropic hydrogel; a similar equation holds for hydrogels in an anisotropic state but with normal and tangential stretches that differ from one another.) Similar behaviour may explain the two-stage dynamics observed in the experiments of Sato Matsuo & Tanaka (Reference Sato Matsuo and Tanaka1988) for a swelling thermo-responsive hydrogel just above the transition temperature.

Although the shaded coexistence and spinodal regions illustrated in figure 3 are valid only for an isotropically swollen hydrogel, and as the hydrogel swells or shrinks it will in general not remain isotropically swollen, they are still useful in determining when this delayed front formation may occur. This is because, in this enforced spherical symmetry, the centre of the sphere must always remain isotropically stretched – the radial and angular stretches converge near the centre, $\lambda _r-\lambda _\theta \to 0$ as $r\to 0$, since the stretches are given by (3.19a–c) and $u/r\to \partial u/\partial r$ as $r\to 0$ (consider a Taylor expansion close to $r=0$ with $u(r=0)=0$). As such, any trajectory (see figure 4) at fixed temperature that begins outside the coexistence region and passes through the coexistence curve will be expected to undergo (delayed) phase separation. For example, shrinking the ANB hydrogel from 302 to 306 K resulted in delayed front formation (see Appendix B).

Beyond this, we have not presented any results that started in the spinodal region. Simulations in this case exhibit spontaneous phase separation, with shrunken regions appearing in the interior of the hydrogel (see Appendix B). Similar dynamics has been observed in one-dimensional hydrogels (Hennessy et al. Reference Hennessy, Münch and Wagner2020). However, we do not explore this more fully because the spherically symmetric constraint of our model causes unrealistic results: the internal phase separation observed is not likely to occur at a given radius all at once to form shrunken concentric spherical shells, but instead occurs in small localised regions across the whole hydrogel that break the spherical symmetry.

We have therefore characterised and classified a range of different dynamic behaviours at fixed temperature using the coexistence and spinodal regions shown in figure 3. The importance of the coexistence and spinodal conditions has been highlighted previously in thermo-responsive hydrogels (e.g. Sekimoto Reference Sekimoto1993), and our results are an extension of the predictions of Tomari & Doi (Reference Tomari and Doi1995), who found phase separation behaviour inside the coexistence region that caused a small range of apparently stable equilibrium states close to the transition to be effectively unstable in practice. We have demonstrated that this dynamic behaviour extends beyond the transition region, and found where it occurs.

The results from our two model hydrogels suggest that phase separation is much common for shrinking hydrogels than those that are swelling. In many scenarios, such as during the expulsion phase of a drug-laden hydrogel, it is particularly important to understand the shrinking behaviour and dynamics of these thermo-responsive hydrogels, and we have seen that front propagation is a key feature of shrinking gels. We therefore turn to investigate the dynamics of the front in more detail now.

5.1. Step function approximation

Solutions that begin in the coexistence region form a sharp front on short time scales. Once the front has formed, profiles of the porosity suggest that it is very close to being uniform in space on either side of the front (see figures 6 and 9). In figure 11(a), we plot the porosity as a function of radius at an arbitrarily chosen time, $t=0.2$, from the simulation for figure 6. We see that this solution looks like a step function, and similar profiles are seen at other times and in the other shrinking front solutions.

Figure 11. (a) Porosity for the ANB shrinking front of figure 6 at time $t={0.2}$, compared to a step function approximation that has levels matched to the initial equilibrium value at the centre and so that the leading-order osmotic pressure is zero towards the edge. (b) The perturbation in porosity from the step function in the numerical solution (solid curves) compared to the calculated values of ${{\phi }^{(1)}_\pm }$ from solving (5.14) (dashed), also at time $t=0.2$. The position of the front is denoted by a vertical dotted line.

To justify the formation of this step function porosity profile, we note that in our system, the parameter relating the relative importance of the mixing and elastic contributions is large: $\varOmega \gg 1$. This means that in general, the osmotic pressure, given by (3.17), dominates the total stress in the mixture. However, we note that the osmotic pressure is simply a function of the porosity multiplied by the large factor $\varOmega$, i.e. we can write $\varPi =\varOmega \,P(\phi )$ for a function $P$. Therefore, any $O(1)$ gradients in the porosity will be expected to result in $O(\varOmega )$ gradients in the osmotic pressure, and hence also in the total radial stress and chemical potential gradients due to the stress balance, (3.15). These large stresses will be relaxed on a fast time scale of $O(\varOmega ^{-1})$, so on an $O(1)$ time scale, we expect to see solutions with constant osmotic pressure and porosity.

Based on this observation, we look for solutions to the poro-elastic model with a front at $r=R_f$ using an expansion in terms of the small parameter $\varOmega ^{-1}\ll 1$:

(5.1)\begin{equation} \phi(r,t) = \begin{cases} {\phi}^{(0)}_-{+} \varOmega^{{-}1}\,{\phi}^{(1)}_-(r,t) + O(\varOmega^{{-}2}) & \text{for }0\leq r< R_f,\\ {\phi}^{(0)}_+{+} \varOmega^{{-}1}\,{\phi}^{(1)}_+(r,t) + O(\varOmega^{{-}2}) & \text{for }R_f< r\leq a, \end{cases} \end{equation}

where the $-$ terms refer to the solution inside the front, and the $+$ terms to solutions outside the front.

For the leading-order porosity, we take the inner value as the initial equilibrium, ${\phi }^{(0)}_-=\phi (r=0,t=0)$, as we assume that the inner region has not been able to expel any fluid during the initial stages of front formation. For the porosity in the outer shell, we expect a uniform value that satisfies the no-radial-stress boundary condition at the edge. Since the osmotic pressure is dominant, this would mean that the porosity is such that the leading-order osmotic pressure vanishes in the entire region, $P({\phi }^{(0)}_+)=0$. We therefore have a shrunken shell where the leading-order osmotic pressure vanishes. A step function with these two levels is plotted in figure 11(a) against the numerical solution, and we see good agreement. The solid line in figure 11(b) shows that the difference between the numerical result and this step function approximation is only a few per cent.

We now calculate the leading-order contribution to the strains and stresses in the hydrogel. In each region, the local volume change (relative to the dry reference state) is found from expanding (2.5) to give

(5.2)\begin{equation} J = \frac{1}{1-{\phi}^{(0)}_\pm} + O(\varOmega^{{-}1}), \end{equation}

and the radial displacement, calculated using (3.18), is

(5.3)\begin{equation} u =\begin{cases} \left[ 1-\left(1-{\phi}^{(0)}_-\right)^{1/3} \right]r + O({\varOmega^{{-}1}}) & \text{for }0\leq r< R_f,\\ \left[ 1-f \right]r + O({\varOmega^{{-}1}}) & \text{for }R_f< r\leq a. \end{cases} \end{equation}

Here, we have defined the function

(5.4)\begin{equation} f(r) = \left[\left(1-{\phi}^{(0)}_+\right)+ \left({\phi}^{(0)}_+{-}{\phi}^{(0)}_-\right) \left( \frac{R_f}{r} \right)^3 \right]^{1/3} \end{equation}

for $R_f\leq r\leq a$, which takes the values $f(R_f)=(1-{\phi }^{(0)}_-)^{1/3}$ and $f(a)=1/a$. Note that $(4{\rm \pi} r^3/3)\times f(r)^3$ is equal to the total volume of solid material inside $r$ (when $r>R_f$), so $f^3$ could be viewed as the average solid fraction within a sphere of radius $r$ once the front has formed.

We then use the displacements to calculate the stretches $\lambda$ from (3.19a–c), which are

(5.5)\begin{align} \lambda_\theta = \begin{cases} \left(1-{\phi}^{(0)}_-\right)^{{-}1/3} + O({\varOmega^{{-}1}}) & \text{for }0\leq r< R_f,\\ f^{{-}1} + O({\varOmega^{{-}1}}) & \text{for }R_f< r\leq a, \end{cases} \end{align}

(5.6)\begin{align} \lambda_r = \begin{cases} \left(1-{\phi}^{(0)}_-\right)^{{-}1/3} + O({\varOmega^{{-}1}}) & \text{for }0\leq r< R_f,\\ \dfrac{f^2}{1-{\phi}^{(0)}_+} + O({\varOmega^{{-}1}}) & \text{for }R_f< r\leq a. \end{cases} \end{align}

The Terzaghi stresses are calculated from the strains using (3.16), giving

(5.7)\begin{align} \sigma'_\theta = \begin{cases} \left(1-{\phi}^{(0)}_-\right)^{1/3} - (1-{\phi}^{(0)}_-) + O({\varOmega^{{-}1}}) & \text{for }0\leq r< R_f,\\[6pt] \left(1-{\phi}^{(0)}_+\right)\left(f^{{-}2}-1\right) + O({\varOmega^{{-}1}}) & \text{for }R_f< r\leq a, \end{cases} \end{align}

(5.8)\begin{align} \sigma'_r = \begin{cases} \left(1-{\phi}^{(0)}_-\right)^{1/3} - (1-{\phi}^{(0)}_-) + O({\varOmega^{{-}1}}) & \text{for }0\leq r< R_f,\\ \dfrac{f^4-\left(1-{\phi}^{(0)}_+\right)^2}{1-{\phi}^{(0)}_+} + O({\varOmega^{{-}1}}) & \text{for }R_f< r\leq a. \end{cases} \end{align}

The (dominant) osmotic pressure, given in (3.17), can be expanded in powers of $\varOmega ^{-1}$ to get

(5.9)\begin{equation} \varPi = {\varOmega\,P\left({\phi}^{(0)}_\pm\right) + P'\left({\phi}^{(0)}_\pm\right) {\phi}^{(1)}_\pm{+} O(\varOmega^{{-}1})}, \end{equation}

where we recall that ${\phi }^{(0)}_+$ satisfies $P({\phi }^{(0)}_+)=0$, and note that the function $P$ and its first derivative are given by

(5.10)$$\begin{gather} P(\phi) ={-} \log\phi - (1-\phi) - (\chi_0-\chi_1)(1-\phi)^2 - 2\chi_1(1-\phi)^3 , \end{gather}$$

(5.11)$$\begin{gather}P'(\phi) ={-}\frac{1}{\phi} + 1 + 2(\chi_0-\chi_1)(1-\phi) + 6\chi_1(1-\phi)^2 . \end{gather}$$

The leading-order contribution to the chemical potential is then calculated using (3.15) – noting that $f'(r)=-({\phi }^{(0)}_+-{\phi }^{(0)}_-)R_f^3/r^4f^2$ and simplifying – to obtain

(5.12)\begin{equation} \frac{\partial \mu}{\partial r} =\begin{cases} - P'\left({\phi}^{(0)}_-\right) \dfrac{\partial {\phi}^{(1)}_-}{\partial r} + O(\varOmega^{{-}1}) & \text{for }0\leq r< R_f,\\[6pt] \dfrac{-2\left({\phi}^{(0)}_+{-}{\phi}^{(0)}_-\right)^2 R_f^6}{\left(1-{\phi}^{(0)}_+\right)f^2 r^7} - P'\left({\phi}^{(0)}_+\right) \dfrac{\partial {\phi}^{(1)}_+}{\partial r} + O(\varOmega^{{-}1}) & \text{for }R_f< r\leq a. \end{cases} \end{equation}

To sustain a profile that has uniform porosity on either side of a travelling front, we expect the total fluid flux to be uniform in space at any given time, so that the fluid pushed outwards by the front moving inwards (causing the gel to contract locally as the front passes by) is equal to the amount of fluid expelled from the boundary. For this to be the case, $r^2\,\partial \mu /\partial r$ must be uniform in each region. Since $\partial \mu /\partial r=0$ at the centre of the sphere, the section of swollen core in our solution must therefore have no flux, and the chemical potential should be of the form

(5.13)\begin{equation} \frac{\partial \mu}{\partial r} = \begin{cases} 0 & \text{for }0\leq r< R_f,\\ A/r^2 & \text{for }R_f< r\leq a, \end{cases} \end{equation}

for some $A$ that we wish to determine – its value quantifies the total (inward) fluid flux from the gel sphere. Note that $A$ is independent of $r$ but may vary in time, $A=A(t)$.

Equating (5.12) and (5.13), we find a differential equation for the porosity perturbation in the shrunken shell of the gel, ${{\phi }^{(1)}_+}$, which is valid for $R_f< r< a$:

(5.14)\begin{equation} {P'\left({\phi}^{(0)}_+\right) \frac{\partial {\phi}^{(1)}_+}{\partial r}} = \frac{-2\left({\phi}^{(0)}_+ - {\phi}^{(0)}_-\right)^2 R_f^6}{\left(1-{\phi}^{(0)}_+\right)f(r)^2 r^7} - \frac{A}{r^2}. \end{equation}

At each point in time, this is a first order differential equation in $r$ for ${{\phi }^{(1)}_+}(r,t)$ with one unknown constant, so we require two boundary conditions.

The first of these is simply the radial stress balance at the outer boundary: $\sigma _r=-\mu _0$ (i.e. $\sigma '_r=\varPi$) at $r=a$. This gives one condition on ${{\phi }^{(1)}_+}$ at the sphere edge, namely

(5.15)\begin{equation} {P'\left({\phi}^{(0)}_+\right) {\phi}^{(1)}_+} = \frac{1}{a^4\left(1-{\phi}^{(0)}_+\right)} - \left(1-{\phi}^{(0)}_+\right) \quad \text{at } r=a. \end{equation}

The second condition imposes that the radial stress is continuous across $r=R_f$: $[\sigma _r]_-^+=0$. Enforcing that the chemical potential is continuous across the front (i.e. $[\mu ]_-^+=0$), we then find a condition for the porosity perturbation at the front:

(5.16)\begin{align} P'\left({\phi}^{(0)}_+\right) {\phi}^{(1)}_+ &= \left(1-{\phi}^{(0)}_-\right)^{4/3} \left[\frac{1}{1-{\phi}^{(0)}_+} - \frac{1}{1-{\phi}^{(0)}_-} \right] \nonumber\\ &\quad + \left({\phi}^{(0)}_+ - {\phi}^{(0)}_-\right) + P\left({\phi}^{(0)}_- \right) \quad \text{at } r=R_f. \end{align}

Given a sphere radius $a$, front position $R_f$, and uniform porosities ${\phi }^{(0)}_\pm$, at any point in time, we can therefore solve the boundary value problem given by the differential equation (5.14) and boundary conditions (5.15) and (5.16), to determine the porosity perturbation ${{\phi }^{(1)}_+}$ and (perhaps more importantly) the value of $A$; hence we can calculate the radial fluid flux. Note that the differential equation (5.14) could be integrated to determine ${{\phi }^{(1)}_+}$ analytically and give a transcendental equation for $A$, but we choose instead to solve (5.14) numerically using MATLAB's in-built boundary value problem solver bvp4c.

5.2. Comparison to the numerical solution

Taking the values for the porosity and sphere radius used in figure 11(a) at time ${t=0.2}$, we now plot the leading-order strains and stresses, calculated from (5.5)–(5.8), in figures 12(a,b). Again, we find good agreement between the step function approximation and the numerical solutions.

Figure 12. Comparison of the step function approximation and the numerical solutions for (a) the stretches, (b) the Terzaghi (elastic) stresses, and (c) the chemical potential gradient after a front has formed. Here, solid curves show the full numerical results for the multi-valued shrinking front, as presented in figure 6, at a dimensionless time $t={0.2}$. The approximate solutions – calculated after measuring the front position and the sphere radius – are shown as dashed curves. In (a,b), the orange (lower) solid curves denote the radial stretches and Terzaghi stresses, whereas the blue (upper) solid curves represent the angular stretches and Terzaghi stresses.

We can also consider calculating the chemical potential gradient $\partial \mu /\partial r$ given in (5.13) by solving the boundary value problem (5.14) for $A$ and ${{\phi }^{(1)}_+}$. Using the same parameters as in figures 11(a) and 12(a,b), we determine $A={-3.25}$, and obtain the chemical potential gradient shown in figure 12(c). We see here that the corrected solution for the chemical potential with a constant fluid flux approximates the numerical solution well.

In solving (5.14), we also determine ${{\phi }^{(1)}_+}$, and in figure 11(b) we plot the calculated ${{\phi }^{(1)}_\pm }$ alongside the values of $\phi (r)-{\phi }^{(0)}_+$ obtained from the full numerical solution. We find good agreement between the two, suggesting that our analytic solution is indeed capturing the key features of this shrinking gel. (There is a visible discrepancy just ahead of the front in figures 11(b) and 12(c), but we attribute this to errors arising from calculating gradients numerically close to a sharp front.)

5.3. Dynamics of the front

We have seen that at a given time, we can well-reproduce the stresses, strains and chemical potential gradient after observing the porosity at the centre, the sphere radius and the front position. We should therefore be able to evolve numerically the front solution forwards in time, since we can calculate the fluid flux at each time point after solving the boundary value problem (5.14)–(5.16). The sphere radius evolves according to (3.20), and conservation of solid in the sphere gives the evolution of the front by

(5.17)\begin{equation} \dot{R_f} ={-} \frac{\left(1-{\phi}^{(0)}_+\right)a^2 \dot{a}}{\left({\phi}^{(0)}_+ - {\phi}^{(0)}_-\right)R_f^2} ={-} \frac{\left(1-{\phi}^{(0)}_+\right)k \left({\phi}^{(0)}_+\right)}{\left({\phi}^{(0)}_+ - {\phi}^{(0)}_-\right)}\, \frac{A(t)}{R_f(t)^2}. \end{equation}

In figure 13(a), we show the full numerical results for the evolution of the front position and sphere radius against the evolved step function solution, using the same parameters as in figure 6. The evolution of this solution was initiated after observing the sphere radius and front position from the numerical solution at time $t={0.001}$; the equations are valid only once the front has formed, so it is not possible to initiate the simulation from $t=0$ (we require $R_f< a$ to solve the boundary value problem (5.14)). Similarly, figure 13(b) shows the evolution of the analytic solutions compared to the results of figure 9, started at time $t=0.01$.

Figure 13. Dynamics of a shrinking hydrogel once a front has formed, comparing the full numerical solution (solid curves) with the step function approximation (dashed curves) for (a) the ANB shrinking front of figure 6, and (b) the HHT shrinking front of figure 9. For the numerical solution, the hydrogel sphere radius is shown in blue, whilst the front position is plotted in orange. In (a), the step function solution is initiated from the numerics at a time $t=0.001$; in (b), it is at $t=0.01$. The dotted lines illustrate the radius and predicted time at which $80\,\%$ of the fluid has been expelled, as described in the drug delivery application of § 5.4.

The evolution of the front and sphere edge appears to be reasonably well approximated by the step function solution. The shapes of the curves are matched well, suggesting that this observed core–shell behaviour in the porosity is indeed dominating the dynamics. However, we note that the step function approximation slightly miscalculates the time scale for the front to reach the centre of the sphere, particularly in figure 13(b); calculating the next order correction to the solution using the solution for ${{\phi }^{(1)}_+}$ may improve the accuracy of this approximation and reduce the discrepancy.

5.4. Application: predicting multi-dose strategies for targeted drug delivery

A potential use for this simplified step function solution is in determining dosage strategies for targeted drug delivery. In this scenario, we consider a swollen hydrogel sphere that contains a given quantity of a drug within its interstitial pore space. Upon heating, the hydrogel contracts and expels its previously absorbed fluid, thereby releasing its load of drugs. As we have seen, in many situations we can expect this temperature change to result in front formation and propagation that dominates the dynamics. If we know the equilibrium behaviour of our hydrogel, and therefore the parameters $\varOmega$, $A_i$ and $B_i$, then we can then use our simplified step function front solution to quickly evolve our model system to calculate the (dimensionless) time required for expulsion of a dosage. If any of the inputs change (e.g. starting or final temperature, or gel sphere size), this model can be updated and new results found in realistic time scales for clinical practice.

This could be extended to develop multi-dose strategies: a drug-laden hydrogel sphere could be actuated on more than one occasion to release specific doses at given times, as illustrated in figure 14. Upon each actuation, we want to release a specified amount of the drug. The question is then: how long do we have to apply the temperature stimulus to achieve the required dosage?

Figure 14. Illustration of a multi-dose drug strategy. A swollen drug-laden hydrogel is heated, causing it to shrink and expel its load. The heat source is removed before all of the drug is released, allowing for a second dose to be administered at a later time. To release a given dosage at each stage, the actuation time of the heating must be controlled carefully.

For an incompressible solid (as was assumed for the polymer in our model), the volume of polymer within the boundary of the sphere remains constant, and any change in total volume must be due to expulsion or absorption of fluid into the pore space. The expelled volume of fluid when the sphere changes from an initial radius $a(0)$ to the current radius $a(t)$, regardless of the spatial arrangement of fluid and solid within the hydrogel sphere in each state, can then be calculated simply as

(5.18)\begin{equation} {\rm \Delta} V = \frac{4{\rm \pi}}{3} \left[a(0)^3-a(t)^3\right], \end{equation}

where the dimensionless expelled volume here is ${\rm \Delta} V = {\rm \Delta} \tilde {V}/\tilde {a}_d^3$, for a dimensional expelled volume ${\rm \Delta} \tilde {V}$.

We consider a situation where we are given a hydrogel with known properties (i.e. we know its equilibrium parameters, permeability, etc.) and are changing between two set temperatures, $\tilde {T}_{{start}}$ and $\tilde {T}_{{end}}$, and we wish to expel a dosage volume $\tilde {v}$ that is less than the total fluid volume within the swollen hydrogel. Using the following procedure, we can calculate the required actuation time:

1. (i) Find the initial equilibrium stretch $\lambda$ by solving the energy minimisation equation (2.7) at the initial temperature $\tilde {T}_{{start}}$. This is the initial dimensionless sphere radius $a(0)=\lambda$.

2. (ii) Evolve the analytic step function solution forward in time until the front reaches the sphere centre.

3. (iii) Calculate the appropriate final radius using (5.18):

   (5.19)\begin{equation} a(t) = \left[ a(0)^3 - \frac{\tilde{v}}{4{\rm \pi}\tilde{a}_d^3/3} \right]^{1/3}. \end{equation}

4. (iv) From the equivalent to figure 13, read off the dimensionless time at which this dimensionless sphere radius is attained.

The dimensional time can then be calculated simply by multiplying by the time scale $\tilde {\tau }$ defined in (3.13).

For example, if we were considering the ANB hydrogel with a temperature change from $\tilde {T}_{{start}}=302$ to $\tilde {T}_{{end}}=308$, then we would obtain figure 13(a). To expel $80\,\%$ of the stored fluid volume, which is equivalent to a dosage volume $\tilde {v}=(4{\rm \pi} \tilde {a}_d^3/3)\times (0.2\,a(0)^3+0.8)$, we find that we must shrink from an initial radius $a=2.29$ to a final radius $a=1.47$ (larger than the final equilibrium radius $a=1.11$). From figure 13(a), we then see that we require a dimensionless time $t=0.092$ from our approximate front solution; our numerics suggest that this should instead be $t=0.089$, which results in a difference in the expelled fluid volume of approximately $1\,\%$. For the same volumes with the HHT hydrogel, the shrinking must occur from an initial radius $a=2.67$ to a radius $a=1.67$; the approximate solution suggests a time $t=0.123$ compared to $t=0.137$ for the numerics. While this error in the timing appears significant, the resulting dosage volume is out by less than $5\,\%$.