2.1. Comparing thermo-responsive LENS with a fully nonlinear model

To show that the predictions of LENS modelling compare well with those of commonly used nonlinear modelling of thermo-responsive hydrogels, we consider the swelling and drying of a poly(N-isopropylacrylamide) (PNIPAM) sphere when heated or cooled around its critical temperature $T_C$ . This problem has been treated extensively in the literature owing to its geometric simplicity and tractability (Matsuo & Tanaka Reference Matsuo and Tanaka1988; Tomari & Doi Reference Tomari and Doi1995; Butler & Montenegro-Johnson Reference Butler and Montenegro-Johnson2022). Since no constitutive laws for material parameters are specified in LENS, we can deduce functional forms for $\Pi (\phi )$ , $\mu _s(\phi )$ and $k(\phi )$ given any model of our choice, including the fully nonlinear models used by other authors.

Starting from the most common approach of choosing a Gaussian-chain nonlinear elastic model for the polymer scaffold coupled with Flory–Huggins theory for interaction between polymer and water molecules (Cai & Suo Reference Cai and Suo2011), we derive the form of $\Pi (\phi )$ in Appendix A,

(2.10) \begin{equation} \Pi (\phi ) = \frac {k_B T}{\Omega _f}\left [\Omega ^{-1}\left (\phi -\phi ^{1/3}\right )-\phi -\textrm {ln}{(1-\phi )}-\phi ^2\chi +\phi ^2(1-\phi )\frac {\partial \chi }{\partial \phi }\right ], \end{equation}

where $k_B$ is the Boltzmann constant, $\Omega _f$ is the volume occupied by a single water molecule, $\Omega$ is the volume of polymer molecules relative to water molecules and $\chi$ is the Flory interaction parameter, quantifying the affinity of polymer chains for water molecules. This osmotic pressure function permits us to deduce the equilibrium polymer fraction as a function of temperature, with a sharp change at the critical temperature $T_C$ . A similar approach gives the shear modulus $\mu _s$ , which is found to be independent of polymer fraction. Then, we fit the measured parameters of Hirotsu et al. (Reference Hirotsu, Hirokawa and Tanaka1987) to find the osmotic modulus, shear modulus and permeability for such gels in our formalism, as detailed fully in Appendix A and the supplementary material. This parameter set was chosen to avoid the complicated hysteresis behaviour seen in other such fitting parameters (for example, those measured by Afroze et al. Reference Afroze, Nies and Berghmans2000), which can be modelled using our approach but we do not discuss here. A more detailed discussion of differences between swelling and drying, and spinodal decomposition between different sets of parameters can be found from Butler & Montenegro-Johnson (Reference Butler and Montenegro-Johnson2022). Further details of the governing equations in both cases, and the precise forms of the non-dimensional times and lengths $t_{BMJ}$ and $r_{BMJ}$ can be found in the supplementary material.

Figure 3.

Plots illustrating the swelling of a hydrogel bead after the temperature is lowered from $308\,\mathrm {K}$ to $304\,\mathrm {K}$ . The parameters used are the same as those used by Butler & Montenegro-Johnson (Reference Butler and Montenegro-Johnson2022), with the fully nonlinear results plotted for comparison, and $t_{BMJ}$ is the non-dimensional time used by Butler & Montenegro-Johnson (Reference Butler and Montenegro-Johnson2022). (a) Evolving polymer fraction with the growth of the radius in the fully nonlinear model shown as a red curve. (b) Porosity profiles at $t_{BMJ} = 0.0001$ , $0.0002$ , $0.0005$ , $0.001$ , $0.0025$ , $0.01$ , $0.05$ , $0.1$ , $0.2$ , $0.5$ and $1$ , with darker blue representing later times. Results from the fully nonlinear model are shown as dashed lines.

We first compute the swelling behaviour of a gel that is initially in equilibrium at $T=308\,\mathrm {K}$ before the temperature is rapidly decreased to $304\,\mathrm {K}$ . This leads to swelling from an initial polymer fraction of $\phi _0 \approx 0.64$ to a much lower value $\phi _0 \approx 0.05$ . Figure 3 shows good quantitative and qualitative agreement with the results of Butler & Montenegro-Johnson (Reference Butler and Montenegro-Johnson2022), with marginally slower growth of the radius but the same diffusive transport of water from surroundings into the bulk of the gel. Repeating the same analysis for smooth drying (where there is no formation of a drying front), we raise the temperature from $304\,\mathrm {K}$ to $307.6\,\mathrm {K}$ , with figure 4 showing the good qualitative, but weaker quantitative, agreement in this case. Our model does, however, capture the rapid initial and later-time drying, with a plateau of slower drying present when $2 \leqslant t_{BMJ} \leqslant 5$ .

Figure 4.

Plots illustrating the drying of a hydrogel bead after the temperature is raised from $304\,\mathrm {K}$ to $307.6\,\mathrm {K}$ , with the same parameters as before and the fully nonlinear solution plotted for comparison. (a) Evolving porosity with the shrinkage of the radius in the fully nonlinear model shown as a red curve. (b) Porosity profiles are shown at $t_{BMJ} = 0$ , $1$ , $2$ , $3$ , $4$ , $5$ , $6$ , $7$ and $8$ , with darker blue representing later times. Results from the fully nonlinear model are shown as dashed lines.

There is a more significant discrepancy in the predictions of LENS and the fully-nonlinear model in this case due to the significant polymer fraction gradients present close to $t_{BMJ} =5$ . In Butler & Montenegro-Johnson (Reference Butler and Montenegro-Johnson2022), the criteria for gel deswelling with phase separation are deduced, and in this smooth deswelling problem, we pass close to a region of parameter space where phase separation can occur. The presence of a nearby equilibrium solution gives a critical slow-down behaviour akin to that discussed by Gomez, Moulton & Vella (Reference Gomez, Moulton and Vella2017), manifesting itself as the plateau of slow drying at intermediate times.

2.1.1. Phase separation and negative diffusivities

Often during the deswelling process a sharp drying front forms, travelling radially inwards through the bead, with the exterior rapidly drying to its final state and the interior remaining relatively swollen until the front reaches the centre. This occurs when trajectories in $(T,\,\phi )$ -space pass through the spinodal or coexistence regions. In the spinodal region, spontaneous phase separation can occur, with the formation of regions of dried polymer surrounded by swollen gel or vice versa as the system equilibrates. The coexistence region is a special case of this, where a dried gel and a swollen one can coexist in thermodynamic equilibrium with a simple sharp boundary (such as a drying front) separating the two. In the present study, we consider both of these effects to be forms of spinodal decomposition, with coexistence a weaker ‘local’ form. In either case, there are sharp differences in $\phi$ across very short distances, as seen especially in cases where there is significant hysteretic behaviour in the equilibrium curve (for example, in the gel parameters measured by Afroze et al. Reference Afroze, Nies and Berghmans2000).

Since large gradients in polymer fraction lead to large deviatoric strains, we expect that our model is unlikely to capture the dynamics of these sharp fronts exactly, since it is dependent on the assumption that these strains remain small. Attempting to replicate this behaviour regardless, through raising the temperature from $304\,\mathrm {K}$ to $308\,\mathrm {K}$ , shows that the polymer diffusivity is, in fact, negative for this case in our model. This leads to spinodal decomposition, with $D(\phi )\lt 0$ the criterion for such behaviour to occur. Figure 5 shows the trajectory of swelling and drying problems in $(T,\,\phi )$ -space for one particular choice of parameters, making it clear why swelling (when the temperature is lowered) never leads to negative diffusivities and why some drying can occur (such as that of figure 4) without entering the spinodal region. In the remainder of this paper, we will consider cases of smooth drying where phase separation does not occur: indeed, taking the linear form of $\Pi$ in (2.3) enforces this.

Figure 5.

A plot of the region in $(T,\,\phi )$ -space where the polymer diffusivity is negative (spinodal region), alongside the equilibrium polymer fraction $\phi _0(T)$ , using the gel parameters of Hirotsu et al. (Reference Hirotsu, Hirokawa and Tanaka1987), the smooth swelling problem of figure 3 is plotted in blue, with the temperature lowered and the spinodal region never approached, and the smooth drying of figure 4 is plotted in yellow. Phase separation occurs, for example, when the temperature is raised to $308\,\mathrm {K}$ and the path to equilibrium passes through the spinodal region, as shown in the example green trajectory.

3.1. Model problem

We consider an infinite tube, symmetric around $z=0$ , formed from thermo-responsive gel, occupying the region $a_0 \lt r \lt a_1$ . The lumen of this tube is filled with water and it is surrounded by water. Initially, the gel is in a swollen state with uniform polymer fraction $\phi \equiv \phi _{00}$ and the temperature is constant everywhere, equal to $T_C - \Delta T$ , below the critical threshold for deswelling. When the temperature is brought above the critical value, the gel will deswell, leading to a shrinkage of the tube and the expulsion of water. This water can be expelled radially out of the tube, carried (slowly) through the gel parallel to the axis, or can be transported axially in the lumen of the tube. Though the deswelling response to the temperature change is still governed by the poroelastic time scale, the tube can be manufactured to be sufficiently thin that shrinkage is rapid and bulk water can be transported much more rapidly through the hollow lumen than would otherwise be the case for a solid cylinder (as in the case investigated by Webber et al. (Reference Webber, Etzold and Worster2023)), so that the gel device acts like a small-scale displacement pump, reacting on a much faster time scale than $t_{pore}$ .

Figure 6.

An illustration of a section of the hydrogel tube in $z\geqslant 0$ , occupying the region $b_0 \lt r \lt b_1$ with a hollow lumen inside. The heat pulse that starts at $z=0$ has spread out here leading to a collapse of the tube, which is fully swollen as $z \to \infty$ . At temperatures below the LCST, the tube occupies its original position $a_0 \lt r \lt a_1$ .

The deswelling of tubes formed from hydrogels has been studied in the past, specifically in the context of water-filled tubes exposed to the air, losing water through their walls as they dry out (Curatolo et al. Reference Curatolo, Lisi, Napoli and Nardinocchi2023). Qualitatively, many of the same phenomena as seen in our model situation are seen here: water is driven radially from a fluid-filled pore, through the thin walls of the tube and then out into the surroundings as the gel deswells. However, notably, our gels are surrounded by water and not air, so the tube is unstressed, since we are effectively imposing zero external chemical potential on the outside of the tube by taking $p=0$ here. This implies that we do not expect to see the strong suction effects seen in air drying, where negative inner pressures arise, which have been shown to lead to a circumferential buckling instability (Curatolo, Nardinocchi & Teresi Reference Curatolo, Nardinocchi and Teresi2018). In our case, we can therefore assume that the shape of the tube will remain axisymmetric for all time.

To illustrate the response of the tube to a temperature field that varies in space and time, we impose a temperature $T_C + \Delta T$ at $z=0$ for all times $t\geqslant 0$ . As this heat pulse spreads out $z \to \pm \infty$ in space, there is a collapse of the tube in regions with $T \gt T_C$ , and water is driven out into the surroundings and towards the still-swollen sections of tube. The radius of the inner lumen is described by $r=b_0(z,\,t)$ whilst the outer radius is $r=b_1(z,\,t)$ , such that the collapsed tube occupies the region $b_0 \lt r \lt b_1$ and $b_i(z,\,0)=a_i$ ( $i=0,\,1$ ). Figure 6 illustrates a section of tube some time after the heat pulse has spread out from $z=0$ , showing collapse behind the heat pulse front. Symmetry arguments imply that we can restrict our attention to $z \geqslant 0$ in most cases, with $b_0$ , $b_1$ , $T$ and $\phi$ all even functions of  $z$ .

3.2. Deformation of the tube

In line with LENS modelling, we assume that all deformation is locally isotropic, and that deswelling leads to a displacement field (relative to the initial state) with axial component $\eta$ and radial component $\xi$ given by

(3.2) \begin{equation} \frac {\xi }{r} \approx \frac {\partial \xi }{\partial r} \approx \frac {\partial \eta }{\partial z} \approx 1 - \left (\frac {\phi }{\phi _{00}}\right )^{1/3}. \end{equation}

Making this assumption requires the polymer fraction field to be independent of $r$ at leading order, an assumption that is reasonable to make in the slender limit of a tube with much larger horizontal length scale than diameter. Since $\eta = 0$ at $z=0$ , owing to symmetry around this point, the leading-order displacement field is

(3.3) \begin{equation} \xi = \left [1 - \left (\frac {\phi }{\phi _{00}}\right )^{1/3}\right ]r \textrm { and } \eta = \int _0^z{\left [1 - \left (\frac {\phi }{\phi _{00}}\right )^{1/3}\right ]\,\mathrm {d}u}. \end{equation}

Using this expression for $\xi$ allows us to write

(3.4) \begin{equation} \frac {b_0}{a_0} \approx \frac {b_1}{a_1} \approx \left (\frac {\phi }{\phi _{00}}\right )^{-1/3}, \end{equation}

and so the local thickness of the tube is proportional to $\phi ^{-1/3}$ .

3.3. Response to changes in temperature

To derive the response of the tube to changes in temperature, it is necessary to solve for the evolution of the gel composition $\phi (r,\,z,\,t)$ . Initially, $\phi \equiv \phi _{00}$ everywhere and the geometry of the problem shows that

(3.5) \begin{equation} \left .\frac {\partial \phi }{\partial z}\right |_{z=0} = 0\quad \textrm { and } \quad \left .\frac {\partial \phi }{\partial z}\right |_{z\to \pm \infty } = 0, \end{equation}

arising from the symmetry of the tube and heat pulse around $z=0$ and the assumption that there is no axial flow through the tube itself (proportional to ${\partial \phi }/{\partial z}$ through (2.7)) as $z \to \pm \infty$ . Indeed, we can simplify the problem further using the symmetry around $z=0$ , and instead solve for the composition in $z \geqslant 0$ alone, reflecting our solution to extend to negative $z$ values.

There has recently been much discussion on the boundary conditions to be applied at the interface of a gel with its surroundings (Xu et al. Reference Xu, Zhang, Young, Yue and Feng2022, Reference Xu, Yue and Feng2024). Significantly, it has been shown that the nature of the tangential stress and velocity boundary conditions can have a significant effect on the dynamics of flows within and without a hydrogel. In addition, there are potential frictional effects as fluid flows cross the water–gel boundary. When flows are significant and the external fluid cannot be assumed quiescent, it is important to choose boundary conditions with care, but in the present study, the poroelastic time scale is sufficiently long that viscous stresses are negligible (Webber & Worster Reference Webber and Worster2023) and the dominant flows are radial, with only small tangential components, so the external fluid is treated as a quiescent bath (even though there are flows, for example, along the axis through the lumen).

At the surface $r=b_1(z,\,t)$ , we assume that there is no radial stress exerted by the tube on its surroundings (and vice versa), so $\mathsf {\sigma }_{rr} = 0$ . Further assuming that large-scale flows in the water bath surrounding the tube are small, we take the fluid pressure to be a constant $p\equiv 0$ outside of the tube. Continuity of pervadic pressure then implies that $p=0$ on $r=b_1$ , and this combines with the condition on $\mathsf {\sigma }_{rr}$ to give

(3.6) \begin{equation} \left .\Pi (\phi )\right |_{r=b_1} = 2\mu _s \epsilon _{rr} = 0\quad \textrm { so }\ \phi (b_1(z),\,z,\,t) = \phi _0, \end{equation}

since the deviatoric strain is zero by our assumption of local isotropy.

On the inside of the tube, we cannot a priori make the assumption of uniform zero pervadic pressure since viscous stresses arising from the lumen flows should be balanced by gradients in $p$ . As discussed in Appendix B, there is an order- $(a_1/L)^2$ correction to the interior pressure field, leading to a mixed boundary condition with a contribution from ${\partial \phi _2}/{\partial R}$ . However, the viscous stresses are much larger on the interior of the gel than on the exterior, owing to the low permeability, and it can be shown, using (B7), that the lumen pressure field has little effect on the gel dynamics. Thence, we can use the same boundary condition on the interior of the gel tube as on the exterior, taking

(3.7) \begin{equation} \phi (b_0(z,\,t),\,z,\,t) = \phi _0. \end{equation}

To describe the evolution of polymer fraction in time as the gel expels water, (2.9) becomes

(3.8) \begin{align} \frac {\partial \phi }{\partial t} + \boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {\nabla }\phi &= \frac {1}{r}\frac {\partial }{\partial r}\left [rD(\phi , T)\frac {\partial \phi }{\partial r}\right ] + \frac {\partial }{\partial z}\left [D(\phi , T)\frac {\partial \phi }{\partial z}\right ] \quad \text {with} \notag \\ D(\phi , T) &= \frac {k}{\mu _l}\left [\frac {\Pi _0(T)\phi }{\phi _0(T)} + \frac {4\mu _s}{3}\left (\frac {\phi }{\phi _{00}}\right )^{1/3}\right ]. \end{align}

There is no intrinsic axial length scale arising from the geometry of this problem, since the tube is infinite in length, but we introduce a length scale $L$ representing the characteristic distance over which temperature variations occur. To simplify the analysis, we make a slenderness assumption that the characteristic axial length scale $L$ is much greater than the characteristic radial length scale $a_1$ . Define $\varepsilon = a_1/L$ and assume that the polymer fraction field only has leading-order axial variation, with radial differences in polymer fraction being of the order $\delta \ll 1$ (arising from our assumption of local isotropy),

(3.9) \begin{equation} \phi = \phi _1(z,\,t) + \delta \phi _2(r,\,z,\,t), \end{equation}

where we have made the arbitrary choice that $\phi _2$ is zero on the midline $r= [b_0(z,\,t)+b_1(z,\,t) ]/2$ , with $\phi _1$ being the polymer fraction on the middle of the tube. Substituting this form of $\phi$ into the evolution equation (3.8) and separating variables, we deduce that $\delta =\varepsilon ^2$ . Therefore, we need only make a relatively weak slenderness assumption, since only $\varepsilon ^2$ need be small.

The material flux $\boldsymbol {q}=q_r\boldsymbol {\hat {r}} + q_z\boldsymbol {\hat {z}}$ is solenoidal and thus $q_r / a_1 \sim q_z / L$ , so $q_r \sim \varepsilon q_z$ . Therefore,

(3.10) \begin{equation} q_r \frac {\partial \phi }{\partial r} = \varepsilon ^2 q_r \frac {\partial \phi _2}{\partial r} \sim \frac {\varepsilon ^3}{a_1} q_z\quad \textrm { and }\quad q_z \frac {\partial \phi }{\partial z} \sim \frac {\varepsilon }{a_1} q_z, \end{equation}

allowing us to neglect radial advection using this slenderness assumption. Introducing the non-dimensional variables $R=r/a_1$ and $Z=z/L$ , the same non-dimensional scalings can be made to the radii $b_0$ and $b_1$ , so

(3.11) \begin{equation} B_0 = \frac {b_0}{a_1} = \ell \left (\frac {\phi _1}{\phi _{00}}\right )^{-1/3}\quad \textrm { and } \quad B_1 = \frac {b_1}{a_1} = \left (\frac {\phi _1}{\phi _{00}}\right )^{-1/3}, \end{equation}

where $\ell = a_0/a_1 \lt 1$ . The leading-order balance of (3.8) is thus

(3.12) \begin{equation} L^2 \frac {\partial \phi _1}{\partial t} + L q_z \frac {\partial \phi _1}{\partial Z} = \frac {1}{R}\left [R D(\phi _1,\,T)\frac {\partial \phi _2}{\partial R}\right ] + \frac {\partial }{\partial Z}\left [D(\phi _1,\,T)\frac {\partial \phi _1}{\partial Z}\right ]. \end{equation}

We now separate variables for $\phi _2$ , since the only term that depends on $R$ is the first diffusive term on the right-hand side. Hence,

(3.13) \begin{equation} \phi _2(R,\,Z,\,t) = \frac {f(Z,\,T,\,t)}{4 D(\phi _1,\,T)}R^2 + g(Z,\,T,\,t) \textrm {ln}{R} + h(Z,\,T,\,t). \end{equation}

By definition, $\phi _2 = 0$ on the midline of the tube $R=(B_0+B_1)/2$ and $\phi _2 = [\phi _0(T)-\phi _1]/\varepsilon ^2$ on $R=B_0$ and $R=B_1$ from boundary conditions (3.6) and (3.7), so $\phi _2$ is symmetric around the midline, with

(3.14) \begin{align} \phi _2(R,\,Z,\,t) &= \frac {4[\phi _0(T)-\phi _1]}{\varepsilon ^2(1-\ell )^2}\left (\frac {\phi _1}{\phi _{00}}\right )^{2/3}\left [R - \frac {1+\ell }{2}\left (\frac {\phi _1}{\phi _{00}}\right )^{-1/3}\right ]^2 \quad \text {with} \notag \\ f &= \frac {16D(\phi _1,\,T)[\phi _0(T)-\phi _1]}{\varepsilon ^2(1-\ell )^2}\left (\frac {\phi _1}{\phi _{00}}\right )^{2/3}. \end{align}

This allows us to deduce the radial structure of the polymer fraction given the value of $\phi _1$ , the polymer fraction on the inside of the tube. This is found by solving the evolution equation (3.12), which is now fully determined up to the total axial flux. Using the approach outlined by Webber et al. (Reference Webber, Etzold and Worster2023),

(3.15) \begin{align} L q_z &= \frac {D(\phi _1,\,T)}{\phi _1}\frac {\partial \phi _1}{\partial Z} + L \left (\frac {\phi _1}{\phi _{00}}\right )^{-1/3}\frac {\partial \eta }{\partial t} \notag \\ &= \frac {D(\phi _1,\,T)}{\phi _1}\frac {\partial \phi _1}{\partial Z} - \frac {L^2}{3}\left (\frac {\phi _1}{\phi _{00}}\right )^{-1/3}\int _0^Z{\left (\frac {\phi _1}{\phi _{00}}\right )^{-2/3}\frac {\partial \phi _1}{\partial t}\,\mathrm {d}u}. \end{align}

This, alongside the form of $\phi _2$ from (3.14), can then be substituted into (3.12), which we non-dimensionalise by introducing the variables

(3.16) \begin{equation} \tau = \frac {k\Pi _{00}t}{\mu _l L^2},\quad \Phi _{0,\,1,\,2,\,\infty } = \frac {\phi _{0,\,1,\,2,\,0\infty }}{\phi _{00}},\quad \mathcal {M} = \frac {\mu _s}{\Pi _{00}},\quad \mathcal {D}(\Phi _1,\,T) = \frac {\mu _l}{k \Pi _{00}}D(\phi _1,\,T). \end{equation}

Then,

(3.17) \begin{align} &\frac {\partial \Phi _1}{\partial \tau } + \frac {\mathcal {D}}{\Phi _1}\left (\frac {\partial \Phi _1}{\partial Z}\right )^2 - \frac {\Phi _1^{-1/3}}{3}\frac {\partial \Phi _1}{\partial Z}\int _0^Z{\Phi _1^{-2/3}\frac {\partial \Phi _1}{\partial \tau }\,\mathrm {d}u} = \mathcal {F} + \frac {\partial }{\partial Z}\left [\mathcal {D}\frac {\partial \Phi _1}{\partial Z}\right ] \quad \text {with} \notag \\ &\qquad \mathcal {D} = \frac {\Pi _0(T)}{\Pi _{00}}\frac {\Phi _1}{\Phi _0(T)} + \frac {4\mathcal {M}}{3}\Phi _1^{1/3}\quad \textrm { and }\quad \mathcal {F} = \frac {16 \Phi _1^{2/3}\mathcal {D}[\Phi _0(T)-\Phi _1]}{\varepsilon ^2(1-\ell )^2}. \end{align}

This is to be solved subject to the initial condition $\Phi _1 \equiv 1$ , and subject to boundary conditions ${\partial \Phi _1}/{\partial Z}=0$ at both $Z=0$ and as $Z \to \infty$ . In the present study, a finite-difference scheme akin to that summarised in the supplementary material of Webber et al. (Reference Webber, Etzold and Worster2023) is used to solve (3.17). From this solution, (3.14) gives the radial polymer fraction structure and the shape of the tube is given by

(3.18) \begin{equation} \ell \Phi _1^{-1/3} \leqslant R \leqslant \Phi _1^{-1/3}, \end{equation}

at leading order in the small parameter $\varepsilon ^2$ . The form of the function $\mathcal {F}$ implies that, for our model to be consistent, $\Phi _1$ must everywhere be close to the piecewise-constant equilibrium polymer fraction $\Phi _0$ , or else our scaling arguments for the terms in the advection–diffusion equation will be invalid. We can check this assumption after calculating the solution to verify the validity of our modelling.

3.4. Response to uniform temperature change

Before studying the response of a hollow tube to a propagating heat pulse, we first consider the case where the temperature is everywhere brought up from below the LCST to $T_C + \Delta T$ at $t=0$ . The response of the tube is axially uniform, evolving following a simplified form of (3.17),

(3.19) \begin{equation} \frac {\partial \Phi _1}{\partial \tau } = \frac {16(\Phi _\infty -\Phi _1)}{\varepsilon ^2(1-\ell )^2}\left (\tilde {\Pi }\frac {\Phi _1^{5/3}}{\Phi _\infty } + \frac {4\mathcal {M}}{3}\Phi _1\right ), \end{equation}

where $\tilde {\Pi } = \Pi _{0\infty }/\Pi _{00}$ . We can use this equation to understand how the material parameters $\Phi _\infty$ , $\mathcal {M}$ and $\tilde {\Pi }$ affect the response time to a change in temperature without the added complication of spatial variations. We know that the polymer fraction on the interior of the tube wall will approach $\Phi _\infty$ as time goes on, with the outside polymer fraction instantaneously reaching this value, but the rate at which this steady state is approached may vary. To measure the rate of deswelling, define the deswelling time scale $\tau _{99}$ as the time taken for

(3.20) \begin{equation} \Phi _1 \geqslant \Phi ^*-\frac {\Phi ^*-1}{100}. \end{equation}

Figure 7.

Plots of the one-dimensional deswelling of a tube when the temperature is uniformly changed when $\Phi _\infty =2$ and $\varepsilon = 0.1$ . This shows the variation of the deswelling time scale $\tau _{99}$ (the time taken for $\Phi _1 \geqslant 1.99$ ) and the approach to steady state for a number of tube thicknesses.

Straightforwardly, it is clear that deswelling is more rapid when there is a greater contrast between $\phi _{00}$ and $\phi _{0\infty }$ , since the bracketed term $\Phi _\infty -\Phi _1$ is greater in magnitude. Thus, gels with more dramatic deswelling will approach their steady states faster. Figure 7(a) shows how the time taken to reach $\Phi _\infty$ depends on the stiffness of the gel (encoded by  $\mathcal {M}$ ) and the relative strength of the osmotic pressure at higher temperatures (encoded by  $\tilde {\Pi }$ ). Stiffer gels resist the formation of deviatoric strains, which arise from differences in polymer fraction, so the interior must deswell to catch up with the outside of the tube, leading to a much faster deswelling process as $\mathcal {M}$ increases. Similarly, larger values of $\tilde {\Pi }$ lead to more rapid interstitial flows driven by pervadic pressure gradients, and so the time to deswell decreases as $\tilde {\Pi }$ increases.

Figure 7(b) illustrates the approach of the polymer fraction on the interior of the tube wall, $\Phi _1$ , to the equilibrium value $\Phi _\infty$ , showing how the approach is more rapid for thinner tubes where there is a shorter distance for water to diffuse out.

3.5. Heat transfer in the system

If, instead of a uniform temperature field, we impose a fixed temperature $T=T_C + \Delta T$ at $Z=0$ , we expect this heat pulse to spread out in the axial direction, symmetrically around the origin, with a deswelling front behind of which $T \gt T_C$ and in front of which $T \lt T_C$ . Modelling the transport of heat through the water, gel scaffold and within the pore space water is a potentially complicated task, and a number of different transport processes must be accounted for, as well as the energetic contributions of swelling, deswelling and deformation (Kaviany Reference Kaviany1995). In the present study, we will consider the simplest possible case, acknowledging that more complicated phenomena such as dispersion will also contribute to heat transfer, but can be reasonably neglected on the assumption that flows through the gel are sufficiently slow.

There is much discussion in the literature of thermoelasticity with specific applications to hydrogels (Cai & Suo Reference Cai and Suo2011; Drozdov Reference Drozdov2014; Brunner, Seidlhofer & Ulz Reference Brunner, Seidlhofer and Ulz2024) and here we take a necessarily simpler model, justifying why deformation and swelling do not contribute to temperature evolution at leading order. In Appendix C, we derive a temperature evolution equation for a gel in the LENS formalism in the absence of external heat supply ( $R=0$ ),

(3.21) \begin{equation} \frac {\partial T}{\partial t} + \boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {\nabla } T = \kappa \nabla ^2 T + \frac {k(\phi )}{\rho c \mu _l}\left |\boldsymbol {\nabla } p\right |^2 + \frac {1}{\phi }\left (\frac {\Pi (\phi )}{\rho c} + T\right )\left (\frac {\partial \phi }{\partial t} + \boldsymbol {q}\boldsymbol {\cdot }\boldsymbol {\nabla }\phi \right ), \end{equation}

where $c$ is the specific heat capacity, $\rho$ is the gel density and $\kappa$ is a spatially averaged thermal diffusivity. In the water surrounding the hydrogel, heat transfer is described by the advection–diffusion equation

(3.22) \begin{equation} \frac {\partial T}{\partial t} + \boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla } T = \boldsymbol {\nabla } \boldsymbol {\cdot } \left (\kappa _w \boldsymbol {\nabla } T\right ) = \kappa _w \nabla ^2 T, \end{equation}

where $\kappa _w$ is the thermal diffusivity of water, assumed to be spatially uniform. In our model, we assume that $\kappa \approx \kappa _w$ . Certainly, this is true in regions where $\phi \ll 1$ and the gel remains swollen, since the contributions of diffusivity in the solid are limited here, a statement supported by experiment and molecular dynamics simulation (Xu, Cai & Liu Reference Xu, Cai and Liu2018). In deswollen regions where the polymer fraction is larger, $\kappa$ is likely of the same order of magnitude as $\kappa _w$ , since the thermal diffusivity of polymer chains is of the same magnitude as the thermal diffusivity of water (Freeman, Morgan & Cullen Reference Freeman, Morgan and Cullen1987). However, such regions are a small fraction of the total spatial domain and the fact that the collapsed tube is thin here means that we neglect any variation from $\kappa _w$ in this region.

There are two potential time scales in the heat transfer problem in the gel – the poroelastic time scale $t_{pore}$ of (3.1), and the thermal time scale $t_{therm} = L^2/\kappa$ . Their ratio is

(3.23) \begin{equation} \frac {t_{pore}}{t_{therm}} = \frac {\kappa }{D} = Le, \end{equation}

the Lewis number, representing the ratio of thermal to compositional diffusivities. In the case $Le \gg 1$ , (3.21) simply reduces to the diffusion equation, since all but the first term on the right-hand side depend on the (slow) reconfiguration of the gel scaffold.

We know that flow and deformation of the gel is mediated by the low permeability of the polymer scaffold, with $k \sim 10^{-15}\,\mathrm {m}^2$ , and therefore expect that the transfer of heat by conduction will occur much faster than changes in shape to the gel. The compositional diffusivity $k\Pi _{00}/\mu _l$ typically scales like $10^{-8}\,\mathrm {m}^2\mathrm {s}^{-1}$ , whilst $\kappa \sim 10^{-7}\,\mathrm {m}^2 \mathrm {s}^{-1}$ , so $Le \sim 10$ . In the present study, we restrict our attention to the large- $Le$ limit for simplicity, where heat transfer in both the gel and water can be modelled by

(3.24) \begin{equation} \frac {\partial \theta }{\partial \tau } = \frac{\partial^2 \theta}{\partial Z^2} \quad \textrm { with } \begin{cases} \theta (0,\,\tau ) = 1,\\ \theta \to -1 ,&\text {as $Z\to \infty $} ,\end{cases} \end{equation}

where $\theta = (T-T_C)/\Delta T$ . There are reasonable physical situations where these assumptions do not apply, but we do not consider them here – modelling such cases would require a careful consideration of heat transfer by advection, dispersion and diffusion, as well as incorporating the effect of fast fluid flows into boundary conditions at the gel surface (Xu et al. Reference Xu, Zhang, Young, Yue and Feng2022). Equation (3.24) has a solution in terms of the error function, with

(3.25) \begin{equation} \theta = 2 \textrm {erfc}\left (\frac {Z}{2\sqrt {Le\, \tau }}\right ) - 1, \end{equation}

where $\textrm {erfc}$ is the complementary error function (Abramowitz & Stegun Reference Abramowitz and Stegun1970). To understand the response of the gel to the diffusive heat pulse, we first seek the position of the deswelling front $Z=Z_C$ , where $\theta =0$ . This is found using (3.25), with

(3.26) \begin{equation} \textrm {erfc}{\left (\frac {Z_C}{2\sqrt {Le\, \tau }}\right )} = \frac {1}{2}\quad \textrm { so }\quad Z_C = 2 \textrm {erfc}^{-1}{\left (\frac {1}{2}\right )}{\sqrt {Le\, \tau }} \approx 0.9539 {\sqrt {Le\,\tau }}. \end{equation}

3.6. Response to pulses of heat

Using the model summarised in (3.17), we can compute the mechanisms by which a thermo-responsive gel tube will collapse in response to a temporally and spatially varying temperature field (3.25). Key to the behaviour here is the fact that heat diffuses on a faster time scale than the water can diffuse through the polymer, leading to a smooth front centred on the pulse front $Z_C(\tau )$ . From this point onwards, we will use the parameters in table 1 in all modelling, having discussed the effect of varying $\mathcal {M}$ and $\tilde {\Pi }$ on the one-dimensional deswelling in previous sections. Figure 8 shows the thickness of a tube at different times as heat diffuses and the gel shrinks. Notice that the shrinkage, though rapid, is not instantaneous in time, since the slow diffusion of water out of the walls of the tube sets a delayed response.

Table 1.

Parameter values used in the modelling of drying tubes from § 3.6 onwards, with the effect of changing $\Phi _\infty$ , $\tilde {\Pi }$ and $\mathcal {M}$ discussed in § 3.4.

Figure 8.

Plots of the evolution of a hollow thermo-responsive hydrogel tube with parameters from table 1 and $\ell = 0.25$ . The heat pulse diffuses from left to right, with the gel shrinking behind it.

For the gel to deswell, water must flow from the walls of the tube into the surrounding water, the lumen at the centre of the tube, or through the gel itself parallel to the axis. Clearly, if the walls of the tube are thinner, driving water from the hydrogel is more rapid, since the water has less of a distance to diffuse outwards, and we expect a more rapid response to changes in temperature for larger values of $\ell$ . The more rapid approach to steady state is shown in figure 9(a), where the sharper equilibrium profile is approached more closely around the drying front $Z_C(\tau )$ for thinner tube walls. Assuming that the radial fluxes are locally dominant, (3.17) reduces to the one-dimensional case of § 3.4,

(3.27) \begin{equation} \frac {\partial \Phi _1}{\partial \tau } \approx \frac {16\Phi _1^{2/3}\mathcal {D}}{\varepsilon ^2(1-\ell )^2} \times \begin{cases}\Phi _\infty -\Phi _1, \; &Z\lt Z_C, \\ 1-\Phi _1, \; &Z\gt Z_C,\end{cases} \end{equation}

away from the front at $Z=Z_C$ (where ${\partial \Phi _1}/{\partial Z}$ will be significant). Then, time scales decrease like $(1-\ell )^2$ when $\ell$ is increased. In the opposite limit as $\ell \to 0$ , adjustment happens on the unmodified poroelastic time scale.

Figure 9.

Plots of the interior polymer fraction $\Phi _1$ at $\tau = 10^{-2}$ with the same parameters as in figure 8, showing how the relaxation to the steady state $\Phi =\Phi _0(T)$ around the drying front $Z=Z_C(\tau )$ is much faster for thinner tubes $\ell \to 1$ . These profiles can be approximated by a $\tanh$ function, as in (3.28), with fitting parameter $A(\ell )$ shown in the logarithmic plot on the right.

From figure 8, it is clear that the structure of the solution around $Z=Z_C(\tau )$ appears to propagate like a travelling wave centred on the deswelling front, since the contribution of axial flows through the gel is limited compared with that of radial flows. Therefore, we can consider the quasi-one-dimensional problem in the new coordinate $Z-Z_C$ . The plots in figure 9 suggest that polymer fraction can locally be approximated by a smooth step around $Z=Z_C(\tau )$ , with the steepness a function of thickness $\ell$ . We thus propose that

(3.28) \begin{equation} \Phi _1 \approx \Phi _\infty - \frac {\Phi _\infty -1}{2}\left \lbrace 1 + \tanh {\left [A(\ell )\left (Z-Z_C\right )\right ]}\right \rbrace , \end{equation}

for some scaling factor $A$ , a function of $\ell$ , representing the sharpness of the drying front. Figure 9 shows that $A(\ell ) \sim (1-\ell )^{-1}$ and therefore the thickness of the adjustment region around the front $Z=Z_C(\tau )$ scales like $(1-\ell )$ .

3.6.1. Flow through the walls

Flow in the walls of the tube is driven by diffusive transport of water from more swollen regions to drier regions, with an interstitial fluid velocity

(3.29) \begin{equation} \boldsymbol {u_g} = \frac {D(\phi )}{\phi }\boldsymbol {\nabla }\phi = \frac {k \Pi _{00}}{\mu _l}\left (\frac {1}{L}\frac {\partial \Phi _1}{\partial Z}\boldsymbol {\hat {z}} + \frac {\varepsilon ^2}{a_1}\frac {\partial \Phi _2}{\partial R}\boldsymbol {\hat {r}}\right ) \times \begin{cases} \frac {\tilde {\Pi }}{\Phi _\infty } + \frac {4\mathcal {M}}{3}\Phi _1^{-2/3}, \; &Z\lt Z_C, \\ 1 + \frac {4\mathcal {M}}{3}\Phi _1^{-2/3}, \; &Z\gt Z_C, \end{cases} \end{equation}

at leading order in the aspect ratio. We define a dimensionless radial fluid velocity $U_g$ scaled with $a_1$ divided by the poroelastic time scale and an axial velocity $V_g$ scaled with $L$ divided by the same time scale, so that

(3.30) \begin{equation} \left (V_g,\,U_g\right ) = \left (\frac {\partial \Phi _1}{\partial Z},\,\frac {\partial \Phi _2}{\partial R}\right ) \times \begin{cases} \frac {\tilde {\Pi }}{\Phi _\infty } + \frac {4\mathcal {M}}{3}\Phi _1^{-2/3}, \; &Z\lt Z_C, \\ 1 + \frac {4\mathcal {M}}{3}\Phi _1^{-2/3}, \; &Z\gt Z_C. \end{cases} \end{equation}

Figure 10 illustrates an example flow field through the walls of the gel, with flow from more swollen to less swollen regions. In the dried region behind the temperature front, radial fluxes are outwards as water is driven out of the shrinking gel, with fluid transported axially towards the drier regions to the left. In $Z\gt Z_C$ , however, fluxes are inwards towards the gel. To understand why this is, notice that the gel is more swollen on its interior than exterior when $Z\lt Z_C$ (as the tube fully dries from the outside in) and more swollen on its exterior than interior when $Z\gt Z_C$ (as the tube is fully swollen on $R=B_1$ with some loss of fluid in the interior due to axial fluxes towards the drier tube). Hence, there needs to be water drawn in from the surrounding fluid to replenish these regions.

Figure 10.

A plot at $\tau = 0.02$ of a drying gel tube with the same parameters as in figure 8. The colours represent the polymer fraction field, with arrows in the gel showing the direction and magnitude of the interstitial flow field $\boldsymbol {u_g}$ , as defined in (3.29). The arrows within the lumen show the flow within the tube, with the form of (3.34).

Figure 11.

Plots close to $Z=Z_C(\tau )$ when $\tau = 0.025$ , illustrating dominant radial flows when the gel is thinner ( $\ell =0.75$ ) versus the thicker ( $\ell = 0.25$ ) gel. In all other regards, the parameters are the same as in figure 10. Notice the directional change either side of the drying front.

In general, therefore, the tube draws water inwards ahead of the deswelling front and then expels the water behind this front. This is shown in detail in figure 11, where the dominance of radial fluxes in thinner gel layers is also clear.

3.6.2. Flow in the lumen

There is also a flow of water through the lumen of the tube, resulting from mass conservation and the redistribution of fluid as the tube dries out. Assuming that the flow within the tube lumen can be described by a Stokes flow $\boldsymbol {v}=v_r\boldsymbol {\hat {r}} + v_z\boldsymbol {\hat {z}}$ with viscosity $\mu _l$ ,

(3.31) \begin{equation} \mu _l\nabla ^2 \boldsymbol {v} - \boldsymbol {\nabla } p = \boldsymbol {0}\quad \textrm { and }\quad \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {v} = 0. \end{equation}

Using the slenderness approximation and assuming that pervadic pressure is independent of $r$ , axial derivatives can be neglected in the radial component of the momentum equation and it is found that the radial component of the velocity must be linear in $r$ if it is to be regular at $r=0$ , and so

(3.32) \begin{equation} v_r = C(z) r \quad \textrm { and }\quad v_z = -2\int _0^z{C(z^{\prime })\,\mathrm {d}z^{\prime }}, \end{equation}

assuming that $v_z = 0$ at $z=0$ by symmetry. We then non-dimensionalise to find a radial velocity $U$ scaled with $a_1$ divided by the poroelastic time scale, and an axial velocity $V$ scaled with $L$ divided by the poroelastic time scale. The radial velocities in the gel are given by

(3.33) \begin{equation} U_g = \frac {8\Phi _1^{2/3}\mathcal {D}(\Phi _1)\left [\Phi _0(T)-\Phi _1\right ]}{\varepsilon ^2(1-\ell )^2}\left (R - \frac {1+\ell }{2}\Phi _1^{-1/3}\right ), \end{equation}

and therefore, matching velocities at $R=\ell \Phi _1^{-1/3}$ ,

(3.34) \begin{equation} U = -\frac {4\Phi _1\mathcal {D}\left [\Phi _0(T)-\Phi _1\right ]}{\varepsilon ^2(1-\ell )^2}R\quad \textrm { and }\quad V = \frac {8}{\varepsilon ^2(1-\ell )^2}\int _0^Z{\Phi _1\mathcal {D}\left [\Phi _0(T)-\Phi _1\right ]\,\mathrm {d}Z^{\prime }}. \end{equation}

Behind the temperature front, flow is radially inwards, since the gel dries into its interior, and this drives a forwards flow along the tube by mass conservation. The magnitude of this flow decreases ahead of the front as there is a weak outward radial flow here.

We can gain some insight into the nature of the axial fluid transport $V$ by considering the fitted form of (3.28), from which we can calculate the evolution of velocity in time. Figure 12 shows how fluid, stationary far behind the temperature front, is driven in the same direction as the thermal front, with a maximum axial velocity $V_{ {max}}$ at $Z=Z_C(\tau )$ . The velocity the decays to a constant value in the tube ahead of the front, which is non-zero by mass conservation. Figure 12(b) shows that the height of this axial flow pulse scales like $A(\ell )^{-1} \sim 1-\ell$ . Since the thickness of this adjustment region scales like $A(\ell )$ , the height of the pulse increases with $\ell$ , but its width decreases with $\ell$ , such that the total flux carried by the pulse is constant as $\ell$ is varied.

Figure 12.

Plots showing the approximate axial velocity $V$ (computed using the form of $\Phi _1$ in (3.28)) for the parameters in table 1, showing how fluid travels in a pulse from the left to the right, with the height of the pulse inversely proportional to the fit parameter $A(\ell )$ .

3.7. Summary of results

In this section, we have used the conceptually simple model for a responsive tubular pump to illustrate a number of results attainable using the responsive LENS formalism that can then be applied to the design of more complicated responsive hydrogel devices. First, we illustrated how the relative impermeability of hydrogels allows for the gel dynamics problem to be decoupled from the fluid flow within cavities (in this example, the lumen of the gel tube), so that the induced pumping flows can be straightforwardly deduced as an output of our modelling, as justified in Appendix B. This allowed us to construct an explicit model for the deswelling of responsive gels as the temperature is changed, quantifying exactly how this deswelling is more rapid for tubes that are thinner $\sim (1-\ell )^2$ and showing how response times can be tuned by varying $\ell$ , the thickness parameter.

Furthermore, the assumption of slenderness and the separation of time scales between slower poroelastic deformation and faster transport of heat allows for decoupling between the thermal problem and the poroelastic problem, and so modelling the behaviour of a thermo-responsive tube to a time- and space-varying temperature field is possible in this framework. We have seen how a tube relaxes to a new equilibrium state around a thermal front and can quantify the spatial structure of the adjustment region, which is smoother and less well-defined for thicker tubes than thinner tubes that approach the equilibrium faster. Perhaps most instructively, the LENS model gives clear expressions for the interstitial flow velocities both along the axis of the tube and out of the walls, as well as the velocities within the hollow lumen, permitting predictions of the nature of the axial pumping fluxes to be made when a tube is heated from one point.