2.1. The Cauchy strain tensor

Our model, as mentioned previously, is linear-elastic relative to a reference state of isotropic swelling, with all nonlinearities encapsulated by the swelling and drying state of the material. There are a number of strain measures used in finite-strain elasticity, incorporating nonlinear effects and generalising the Cauchy strain tensor of linear elasticity to capture the influence of potentially large displacements. Many theories simply use the deformation gradient tensor $\boldsymbol{\mathsf{F}}$ or the Cauchy–Green deformation tensors $\boldsymbol{\mathsf{F}}\boldsymbol{\mathsf{F}}^{{T}}$ and $\boldsymbol{\mathsf{F}}^{{T}}\boldsymbol{\mathsf{F}}$. The choice of strain measure is largely arbitrary, and is usually made to simplify the constitutive relation for the stress. Our model requires a linear relationship between stresses arising from deviatoric strains and these deviatoric strains themselves, and we therefore use the Cauchy strain tensor of linear elasticity. We show in Appendix B that, even starting from a fully hyperelastic model, once the assumption of small deviatoric strain is applied, the stress can be written in terms of the deviatoric Cauchy strain in addition to a nonlinear isotropic part.

When defining our strain tensor it is also important to consider carefully whether the quantities under consideration are Lagrangian or Eulerian. In linear elasticity, this is less important because the small deformations are the same to leading order in both descriptions of the problem whilst here, with potentially large total strains, this is no longer the case. We start by defining the displacement $\boldsymbol {\xi } = \boldsymbol {x} - \boldsymbol {X}$. The Cauchy strain tensor is then defined by

(2.7)\begin{equation} \boldsymbol{\mathsf{e}} = \tfrac{1}{2}\left[\boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\xi}+\left(\boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\xi}\right)^{\rm T}\right], \end{equation}

where $\boldsymbol {\nabla }_{\boldsymbol {x}}$ denotes a gradient taken with respect to the deformed, Eulerian, coordinates. Even though, for large deformations, it may appear that a Lagrangian description of the gel would be more tractable when dealing with swelling problems, it is noted by MacMinn et al. (Reference MacMinn, Dufresne and Wettlaufer2016) that an Eulerian description often tends to be preferable for poroelastic problems, since the equations governing the fluid flow tend to be stated in an Eulerian framework. Coupling is easier if one consistently uses the same description. We can calculate the value of $\boldsymbol {\nabla }_{\boldsymbol {x}}\boldsymbol {\xi }$ in terms of $\boldsymbol{\mathsf{F}}$ using the chain rule because

(2.8)\begin{align} \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{\xi} = \boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{\xi}\boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{X} &= \left[\boldsymbol{\nabla}_{\boldsymbol{X}}\boldsymbol{x} - \boldsymbol{\mathsf{I}}\right]\boldsymbol{\cdot}\boldsymbol{\nabla}_{\boldsymbol{x}}\boldsymbol{X} \nonumber\\ &= \boldsymbol{\mathsf{F}}\boldsymbol{\mathsf{F}}^{-1} - \boldsymbol{\mathsf{F}}^{-1} \nonumber\\ &= \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{F}}^{-1}. \end{align}

Again working only to first order in the small deviatoric strain, we use the expansion of the inverse of a matrix around the identity (Petersen & Pedersen Reference Petersen and Pedersen2012),

(2.9)\begin{equation} \left(\boldsymbol{\mathsf{I}} + \epsilon \boldsymbol{\mathsf{A}}\right)^{-1} = \boldsymbol{\mathsf{I}} - \epsilon \boldsymbol{\mathsf{A}} + O(\epsilon^2), \end{equation}

to show that

(2.10)\begin{equation} \boldsymbol{\mathsf{F}}^{-1} = \frac{1}{\mathcal{F}}\boldsymbol{\mathsf{I}} - \frac{1}{\mathcal{F}^2}\boldsymbol{\mathsf{f}} + O(\boldsymbol{\mathsf{f}}^2). \end{equation}

This gives an expression for the Cauchy strain tensor, split into an isotropic strain and a deviatoric strain ${\boldsymbol {\epsilon }}$ with $\operatorname {tr}{\boldsymbol {\epsilon }} = 0$,

(2.11)\begin{equation} \boldsymbol{\mathsf{e}} = \left[1-\left(\frac{\phi}{\phi_0}\right)^{1/n}\right]\boldsymbol{\mathsf{I}} + {\boldsymbol{\epsilon}} \quad \text{where}\ {\boldsymbol{\epsilon}} = \frac{1}{2}\left(\frac{\phi}{\phi_0}\right)^{2/n}\left(\boldsymbol{\mathsf{f}}+\boldsymbol{\mathsf{f}}^T\right). \end{equation}

3.1. Pressure

We start by considering the bulk pressure of an $N$-component colloidal mixture, which is defined thermodynamically by the relation

(3.1)\begin{equation} P =- \left(\frac{\partial U}{\partial V}\right)_{S,M_{\left\lbrace 1,2,\dots,N\right\rbrace}}, \end{equation}

where $U$, $V$ and $S$ are the internal energy, volume and entropy of the mixture, respectively, and $M_i$ ($i=1,\dots, N$) is the mass of the $i$th component. In the case of a hydrogel, there are only two such components, the polymer and the water. However, and more usefully for our purposes, the bulk pressure can also be defined mechanically, as the total isotropic force per unit area exerted on the mixture. Note that no reference is made as to the source of this pressure in terms of the microstructural properties of the hydrogel. Recalling the definition of the Cauchy stress tensor, we write

(3.2)\begin{equation} {\boldsymbol{\sigma}} =-P\boldsymbol{\mathsf{I}} + \boldsymbol{\sigma}_{{dev}} \quad {\rm with}\ P =-\frac{1}{n}\operatorname{tr}{\boldsymbol{\sigma}}, \end{equation}

where $\boldsymbol {\sigma }_{{dev}}$ is the stress deviator tensor, with zero trace.

This bulk pressure can be further separated into a pervadic pressure $p$ and a generalised osmotic pressure $\varPi$, which we henceforth refer to simply as the osmotic pressure. This separation is discussed by Peppin, Elliott & Worster (Reference Peppin, Elliott and Worster2005), where it is seen that the pervadic pressure can be linked to the chemical potential $\mu$ often used in discussions of colloidal mixtures. The pervadic pressure is defined as the pressure measured by a transducer separated from the gel by a semipermeable membrane that allows only water to pass through. This pressure can be identified with the pore pressure of the liquid component of a gel in some existing poroelastic models (Hewitt et al. Reference Hewitt, Nijjer, Worster and Neufeld2016), and it is gradients in $p$, which is also referred to as the Darcy pressure (Peppin Reference Peppin2009), that drive flow through the matrix.

The osmotic pressure $\varPi = P-p$ is the difference between the bulk and pervadic pressures. It has contributions on the micro scale both from elastic stresses in the polymer scaffold and from the intermolecular interactions between the water and polymer molecules. However, we do not consider these explicitly in our continuum model, merely concerning ourselves with the resultant macroscopic effect from a combination of these many physical factors. We expect the osmotic pressure to be positive when $\phi > \phi _0$ and negative for $\phi < \phi _0$, by the definition of the equilibrium polymer fraction. We also expect this pressure to be a function of polymer fraction alone; not only does the osmotic pressure depend solely on the state of swelling or drying, it is also reasonable to assume that isotropic stresses should depend only on isotropic strains, which can be written as a function of polymer fraction $\phi$ using (2.11). Combining this with all that has been discussed previously, the Cauchy stress tensor (3.2) can be written in the form

(3.3)\begin{equation} {\boldsymbol{\sigma}} =-\left[p + \varPi(\phi)\right]\boldsymbol{\mathsf{I}} + \boldsymbol{\sigma}_{{dev}}. \end{equation}

3.1.1. The osmotic modulus

Authors including Doi (Reference Doi2009) and Etzold et al. (Reference Etzold, Linden and Worster2021) introduce an osmotic modulus $K$, defined as an analogue of the bulk modulus $\kappa$ in linear elasticity (Landau & Lifshitz Reference Landau and Lifshitz1986; Chung Reference Chung2007). Under linear elasticity, in which all strains are considered small, $\kappa$ is defined by $-V \,{\rm d}P/{\rm d}V$, where $V$ is the volume of the elastic material, and thus an incompressible material has an infinite value of $\kappa$. The osmotic modulus is defined in an analogous manner with respect to the osmotic pressure, as opposed to the bulk pressure, so that

(3.4)\begin{equation} K =- V \frac{{\rm d}\varPi}{{\rm d}V}. \end{equation}

This reflects the fact that, for an incompressible elastic hydrogel such as those that we are modelling, volume changes only result from osmotic effects (swelling and drying) and not from bulk elastic ones. In this case, because the volume of a gel is constant and proportional to $1/\phi$, we define the osmotic modulus $K(\phi )$ by the expression

(3.5)\begin{equation} K(\phi) = \phi \frac{\partial \varPi}{\partial \phi}. \end{equation}

In the fully linear limit (Doi Reference Doi2009; Etzold et al. Reference Etzold, Linden and Worster2021), $K/\phi$ is taken as constant, linearising around the fully swollen state $\phi = \phi _0$ such that the osmotic pressure is linear in $\phi$,

(3.6)\begin{equation} \varPi(\phi) = K_0\frac{\phi-\phi_0}{\phi_0}, \end{equation}

with osmotic modulus $K(\phi ) = K_0 \phi / \phi _0$. More generally than this, using our definition above, it is possible to linearise around any polymer fraction $\phi = \phi ^*$ and find that, close to this value,

(3.7)\begin{equation} \varPi(\phi) - \varPi(\phi^*) = K(\phi^*)\frac{\phi-\phi^*}{\phi^*}. \end{equation}

The difference between these two linearising approximations to $\varPi$ is illustrated in figure 2, showing how linearising around a given value of $\phi$ gives a much better approximation in the neighbourhood of $\phi ^*$ than assuming an entirely linear form for $\varPi (\phi )$.

Figure 2.

An illustration of a representative nonlinear osmotic pressure $\varPi (\phi )$ in black, taken to be of the form $(\phi /\phi _0)\ln (\phi /\phi _0)$ (Appendix B) alongside the linear approximation used by Doi (Reference Doi2009) in red. Our model allows for linearisation around any given polymer fraction $\phi ^*$, an example of which is shown in blue. The osmotic modulus $K(\phi ) = \phi \varPi '(\phi )$.

3.2. Deviatoric stresses

Recall that we treat the gel, swollen to any polymer fraction, as instantaneously incompressible and linear-elastic in nature. Further, we expect that isotropic strains lead to isotropic stress and that deviatoric strains lead only to deviatoric stresses, which suggests that $\boldsymbol {\sigma }_{{dev}}$ should only depend on ${\boldsymbol {\epsilon }}$. Assuming linearity and local isotropy of the material, the founding assumptions of most linear-elastic theories, we write ${\boldsymbol {\sigma }}_{{dev}} = \boldsymbol{\mathsf{C}} \boldsymbol {:} {\boldsymbol {\epsilon }}$, where $\boldsymbol{\mathsf{C}}$ is a fourth-rank isotropic tensor and $[\boldsymbol{\mathsf{A}}\boldsymbol {:}\boldsymbol{\mathsf{b}}]_{ij} = {\mathsf{A}}_{ijkl}{\mathsf{b}}_{kl}$. By the traceless nature of ${\boldsymbol {\epsilon }}$, this reduces to

(3.8)\begin{equation} {\boldsymbol{\sigma}}_{{dev}} = 2\mu_s {\boldsymbol{\epsilon}}, \end{equation}

where the constant $\mu _s$ is chosen to agree with the definition of the shear modulus in linear elasticity (Landau & Lifshitz Reference Landau and Lifshitz1986; Chung Reference Chung2007).

However, because we expect the material properties of the gel to depend on the water content, or equivalently the polymer fraction $\phi$, we allow $\mu _s$ to be a function of this polymer fraction, resulting in the constitutive relation for ${\boldsymbol {\sigma }}$,

(3.9)\begin{equation} {\boldsymbol{\sigma}} =-\left[p + \varPi(\phi)\right]\boldsymbol{\mathsf{I}} + 2\mu_s(\phi){\boldsymbol{\epsilon}}. \end{equation}

Equation (3.9) has the form of the fully linear model presented by Doi (Reference Doi2009), but this linear-elastic-nonlinear-swelling relation captures linearity in the small deviatoric strains ${\boldsymbol {\epsilon }}$ whilst also allowing for arbitrarily large isotropic swelling strains. This is achieved by allowing both $\varPi$ and $\mu _s$, the two material parameters, to depend nonlinearly on polymer fraction $\phi$. It is possible to relate these two material parameters to existing nonlinear theories including the aforementioned Gaussian-chain/Flory–Huggins approach to finite-strain poroelasticity (see Appendices A and B for examples of this). The utility of this constitutive model is that it is relatively tractable and it describes the behaviour of the gel solely in terms of these two macroscopically measurable parameters.

Given this equation, we can use the Cauchy momentum equation to balance the forces on a gel element and describe its dynamics. For a body force $\boldsymbol {f_b}$, which may be a function of space,

(3.10)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\sigma}} + \boldsymbol{f_b} = \boldsymbol{0} \end{equation}

is the equation governing momentum balance if the acceleration of gel elements is neglected. In the majority of problems considered henceforth, the body force will be taken to be zero for simplicity.

3.3. Comparison with linear poroelasticity

As an aside, we show here how our formulation relates to classical linear poroelasticity, building on the effective-stress framework (Terzaghi Reference Terzaghi1925; Biot Reference Biot1941), and used in a number of existing investigations into the behaviour of two-phase systems, for example Hewitt, Neufeld & Balmforth (Reference Hewitt, Neufeld and Balmforth2015). In a two-phase system, both the solid and the liquid components contribute to the overall stress of the system, and Terzaghi (Reference Terzaghi1925) assumed the total [Cauchy] stress tensor ${\boldsymbol {\sigma }}$ to be a volume-fraction weighted average of the stresses due to the solid phase and the liquid phase, such that

(3.11)\begin{equation} {\boldsymbol{\sigma}} = \phi {\boldsymbol{\sigma}}^{(s)} + (1-\phi){\boldsymbol{\sigma}}^{(l)}. \end{equation}

As is familiar in fluid mechanics, the liquid stress

(3.12)\begin{equation} {\boldsymbol{\sigma}}^{(l)} =-p\boldsymbol{\mathsf{I}} + {\boldsymbol{\tau}}, \end{equation}

where ${\boldsymbol {\tau }}$ is the deviatoric fluid stress (equal to $2\mu _l {\boldsymbol {\varepsilon }}$ in the case of a Newtonian rheology, where $\mu _l$ is the dynamic viscosity and ${\boldsymbol {\varepsilon }}$ is the rate-of-strain tensor). Terzaghi quantified the excess in stress due to the elasticity of the matrix over and above the isotropic pore pressure by writing

(3.13)\begin{equation} {\boldsymbol{\sigma}} = {\boldsymbol{\sigma}}^{(e)} - p \boldsymbol{\mathsf{I}} + (1-\phi){\boldsymbol{\tau}}, \end{equation}

where ${\boldsymbol {\sigma }}^{(e)} = \phi ({\boldsymbol {\sigma }}^{(s)} + p\boldsymbol{\mathsf{I}})$ is the effective stress tensor (Wang Reference Wang2000). However, in hydrogels, we neglect the effect of viscous stress on the total stress tensor, dropping the final term in this relation to give

(3.14)\begin{equation} {\boldsymbol{\sigma}} = {\boldsymbol{\sigma}}^{(e)} - p \boldsymbol{\mathsf{I}}. \end{equation}

This is not without precedent: Hewitt et al. (Reference Hewitt, Nijjer, Worster and Neufeld2016) neglected viscous stresses, expecting the elastic response to be more important over the timescales we aim to model, and we provide a post-hoc scaling justification for this assumption in Appendix C. By comparing (3.14) with (3.9), we see that the pervadic pressure is equal to the pore pressure as understood by Terzaghi and Biot, and that our (generalised) osmotic pressure is the isotropic part of Terzaghi's effective stress ${\boldsymbol {\sigma }}^{(e)}$.

In linear poroelastic models, a linear-elastic constitutive relation is specified for ${\boldsymbol {\sigma }}^{(e)}$, (Detournay & Cheng Reference Detournay and Cheng1993) which separates Terzaghi's effective stress into an isotropic part related to the bulk modulus $\kappa$ and a deviatoric part related to the shear modulus $\mu _s$. Through this, we can draw analogues between a bulk modulus for the matrix and the osmotic modulus, and also note that the shear modulus in our formulation plays exactly the same role as the shear modulus in linear poroelasticity.

4.1. Rewriting the transport equation in terms of polymer fraction alone

In one dimension, ${\boldsymbol {\epsilon }} = \boldsymbol{\mathsf{0}}$ and (4.7) is expressed solely in terms of $\phi$, and alone provides a general equation to describe nonlinear swelling. In higher dimensions, it is desirable to eliminate deviatoric strain from this transport equation so that an explicit knowledge of the displacements $\boldsymbol {\xi }$ is not needed to explain how the polymer fraction evolves. We can do this to leading order in the small deviatoric strains ${\boldsymbol {\epsilon }}$, starting by combining (2.7) and (2.11) to give

(4.8)\begin{equation} \frac{1}{2}\left[\boldsymbol{\nabla} \boldsymbol{\xi} + \left(\boldsymbol{\nabla} \boldsymbol{\xi}\right)^{\mathrm{T}}\right] = \left[1-\left(\frac{\phi}{\phi_0}\right)^{1/n}\right]\boldsymbol{\mathsf{I}} + {\boldsymbol{\epsilon}}. \end{equation}

The trace of this equation shows that

(4.9)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\xi} = n\left[1-\left(\frac{\phi}{\phi_0}\right)^{1/n}\right], \end{equation}

and its divergence gives

(4.10)\begin{align} \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\epsilon}} &= \frac{1}{2}\left[\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\xi} + \boldsymbol{\nabla}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\xi}\right)\right] - \boldsymbol{\nabla}\left[1-\left(\frac{\phi}{\phi_0}\right)^{1/n}\right], \nonumber\\ &=\frac{1}{2}\nabla^2\boldsymbol{\xi}+ \left(1-\frac{n}{2}\right)\boldsymbol{\nabla}\left(\frac{\phi}{\phi_0}\right)^{1/n}. \end{align}

The fact that deformation is, to leading order, isotropic in our model, combined with (4.8), indicates that

(4.11)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\xi} = \left[1-\left(\frac{\phi}{\phi_0}\right)^{1/n}\right]\boldsymbol{\mathsf{I}} + O(\varepsilon) \quad{\rm thus}\ \nabla^2 \boldsymbol{\xi} =-\boldsymbol{\nabla} \left(\frac{\phi}{\phi_0}\right)^{1/n} + O(\varepsilon/L), \end{equation}

where $\varepsilon = \max _{i,\,j}|\epsilon _{ij}|$ and $L$ is a length scale for the problem. Therefore, combining this result with (4.10) shows that

(4.12)\begin{equation} \boldsymbol{\nabla} \left(\frac{\phi}{\phi_0}\right)^{1/n} = O(\varepsilon/L) \quad{\rm so}\ \boldsymbol{\nabla} \phi = O(\varepsilon/L), \end{equation}

or that gradients in $\phi$ are of the same order as gradients in the deviatoric strain and are therefore in this sense ‘small’ when the deviatoric strain is small. Hence, for any given functions $g$ and $h$ of $\phi$,

(4.13a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ g(\phi) {\boldsymbol{\epsilon}}\right] = g(\phi) \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\epsilon}} \quad{\rm and}\quad \boldsymbol{\nabla} \boldsymbol{\cdot} \left[h(\phi) \boldsymbol{\nabla} \boldsymbol{\cdot}{\boldsymbol{\epsilon}}\right] = h(\phi) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\epsilon}} \end{equation}

to leading order in the small parameter $\varepsilon$. For example,

(4.14)\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} \left\lbrace \phi k(\phi) \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\mu_s(\phi){\boldsymbol{\epsilon}}\right]\right\rbrace &= \phi k(\phi) \mu_s(\phi) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{\epsilon}} + O(\varepsilon^2/L)\nonumber\\ &= \phi k(\phi) \mu_s(\phi) \left(1-n\right) \nabla^2 \left(\frac{\phi}{\phi_0}\right)^{1/n}+ O(\varepsilon^2/L), \end{align}

taking the divergence of (4.10). This can be rewritten as a divergence, also using (4.13a,b), because

(4.15)\begin{equation} \phi k(\phi) \mu_s(\phi) \left(1-n\right) \nabla^2 \left(\frac{\phi}{\phi_0}\right)^{1/n} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[(1-n)\phi k(\phi)\mu_s(\phi) \boldsymbol{\nabla} \left(\frac{\phi}{\phi_0}\right)^{1/n}\right] + O(\varepsilon^2/L). \end{equation}

Using this result, neglecting terms of second order and above in the deviatoric strain, the governing equation for polymer transport can be written, in conservation form, as

(4.16)\begin{equation} \frac{{\rm D}_{\boldsymbol{q}}\phi}{{\rm D}t} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\frac{k(\phi)}{\mu_l}\left\lbrace K(\phi) + \frac{2(n-1)}{n} \mu_s(\phi)\left(\frac{\phi}{\phi_0}\right)^{1/n}\right\rbrace\boldsymbol{\nabla} \phi\right]. \end{equation}

This equation shows that the polymer fraction satisfies a nonlinear advection–diffusion equation with a diffusivity made up of both generalised osmotic effects and bulk elastic effects. Note that this equation makes no assumptions regarding the properties of the hydrogel being investigated other than small deviatoric strains; the constitutive relations for the macroscopic material parameters $K(\phi )$, $\mu _s(\phi )$ and $k(\phi )$ are to be determined.

4.2. Comparison with existing models

We have already seen how the transport equation for polymer of (4.7) compares with the model used by Etzold et al. (Reference Etzold, Linden and Worster2021) in the one-dimensional ($n=1$) case, with an absence of deviatoric strains. In this limit, (4.16) becomes

(4.17)\begin{equation} \frac{{\rm D}_{\boldsymbol{q}}\phi}{{\rm D}t} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\frac{k(\phi) K(\phi)}{\mu_l}\boldsymbol{\nabla} \phi\right], \end{equation}

with the nonlinear diffusivity $D(\phi ) = k(\phi )K(\phi )/\mu _l$. It is equivalent to that for a colloid with solid particle fraction $\phi$ (Doi Reference Doi2009; Hewitt et al. Reference Hewitt, Nijjer, Worster and Neufeld2016; Worster, Peppin & Wettlaufer Reference Worster, Peppin and Wettlaufer2021), because $K(\phi ) = \phi \partial \varPi /\partial \phi$, and is also referred to as the poroelastic diffusivity.

Now consider the case where we would expect linear elasticity to hold both in deviatoric and isotropic strains, and linearise around the equilibrium state $\phi = \phi _0$. If material parameters are assumed constant, (4.16) becomes

(4.18)\begin{equation} \frac{\partial\phi}{\partial t} + \boldsymbol{q}\boldsymbol{\cdot}\boldsymbol{\nabla} \phi = \frac{k}{\mu_l}\left[K + 2\left(1-1/n\right)\mu_s\right]\nabla^2 \phi. \end{equation}

Doi (Reference Doi2009), for example, considers the special case in which $n=3$ and $\boldsymbol {q} = \boldsymbol {0}$ and finds this same result: a linear diffusion equation with diffusion coefficient $k(K+4\mu _s/3)/\mu _l$.

5.1. Finding the osmotic modulus: the steady state

Once the system has reached equilibrium at any step $n$, when the height of the gel is $h_n$ and its radial extent is $a_n$, there are no gradients in pervadic pressure. Given that $p=0$ at the water–gel interface by continuity of pervadic pressure, $p \equiv 0$ throughout the gel. Adding (5.1a) and (5.1b) and using the fact that $\epsilon _{xx} + \epsilon _{zz} = 0$, we find that the steady-state force exerted per unit length on the plate, $\sigma _f$, is given by

(5.2)\begin{equation} \sigma_f = 2\varPi(\phi_n) = 2\varPi_n, \end{equation}

where $\phi _n$ is the steady-state polymer fraction. This polymer fraction can be deduced from polymer conservation, with $\phi _0 h_0 a_0 = \phi _n h_n a_n$. Note also that $\varPi$ is, by definition, zero when $\phi = \phi _0$.

5.2. Finding the shear modulus: the initial force

If the height of the gel is then stepped down from $h_n$ to $h_{n+1} = h_n(1-\varepsilon )$, for $\varepsilon \ll 1$, there is no immediate redistribution of water and the polymer fraction therefore remains equal to $\phi _n$. In order to conserve polymer, therefore, the horizontal extent increases to $-a_n' \leqslant x \leqslant a_n'$, where $a_n' = a_n h_n/h_{n+1}$. Working to first order in the small parameter $\varepsilon$, and letting $K_n = K(\phi _n)$, $\mu _{sn} = \mu _s(\phi _n)$ and $k_n = k(\phi _n)$,

(5.3)\begin{equation} \epsilon_{zz} = {\mathsf{e}}_{zz} + \left(\frac{\phi}{\phi_0}\right)^{1/2} - 1, \end{equation}

from (2.11). Given that the total vertical strain relative to the swollen state $\phi \equiv \phi _0$ is $e_{zz} = h_{n+1}/h_0 - 1$, we find that

(5.4)\begin{equation} \epsilon_{zz} = \frac{h_{n+1}}{h_0} + \left(\frac{\phi_n}{\phi_0}\right)^{1/2} - 2. \end{equation}

Now subtracting (5.1a) from (5.1b) and again using $\epsilon _{xx} + \epsilon _{zz} = 0$, we find that the stress exerted on the top plate immediately after compression is given by

(5.5)\begin{equation} \sigma_i = 4\mu_{sn}\left[2 - \frac{h_n}{h_0}\left(1-\varepsilon\right) - \left(\frac{\phi_n}{\phi_0}\right)^{1/2}\right], \end{equation}

allowing us to use a stress measurement to deduce the value of $\mu _{sn}$ because all of the other parameters here are already known.

5.3. Finding the permeability: the transient evolution

It is the permeability of the gel that dictates the rate at which fluid can either flow into or out of the porous polymer scaffold, through the transport equation (4.16). Start by considering the volume-averaged flux $\boldsymbol {q}$, which only has a horizontal component q, given our conceptual framework. Therefore, $\boldsymbol {\nabla }\boldsymbol{\cdot}\boldsymbol {q} = 0$ reduces to $\partial q / \partial x = 0$, and symmetry dictates that $q=0$ at $x=0$. Thus, $q \equiv 0$ through the gel.

As discussed previously, the gel will relax back to a uniform state with $\phi \equiv \phi _{n+1}$ as water is expelled. This value of $\phi _{n+1}$ is set by considering the boundary condition $\sigma _{xx} = p = 0$ at $x=a(t)$, which, using (5.1a) and (5.1b), shows that

(5.6)\begin{equation} \varPi(\phi_{n+1}) = 2\mu_s(\phi_{n+1})\epsilon_{zz}. \end{equation}

Using (5.4) for the vertical deviatoric strain and remarking that, for a sufficiently small incremental step, $\varPi (\phi _{n+1}) \approx \varPi _n + (K_n/\phi _n)(\phi _{n+1} - \phi _n)$,

(5.7)\begin{equation} \frac{K_n}{\phi_n}(\phi_{n+1} - \phi_n) - 2\mu_{sn}\left[\frac{h_n}{h_0}\varepsilon - \frac{1}{2}\left(\frac{\phi_n}{\phi_0}\right)^{1/2}\frac{\phi_{n+1}-\phi_n}{\phi_n}\right] = 0, \end{equation}

providing an implicit expression for the steady-state polymer fraction $\phi _{n+1}$.

We linearise around the mean state $\hat {\phi } = (\phi _n + \phi _{n+1})/2$, with constant mean material parameters $\hat {K} = (K_n + K_{n+1})/2$, $\widehat {\mu _s} = (\mu _{sn} + \mu _{s(n+1)})/2$ and $\hat {k} = (k_n + k_{n+1})/2$. Thus (4.16) becomes a linear diffusion equation

(5.8)\begin{equation} \frac{\partial\phi}{\partial t} = \frac{\hat{k}}{\mu_l}\left[\hat{K} + \left(\frac{\hat{\phi}}{\phi_0}\right)^{1/2}\widehat{\mu_s}\right]\frac{\partial^2\phi}{\partial x^2}. \end{equation}

The horizontal extent of the gel at time $t$, $a(t)$, is set by polymer conservation

(5.9)\begin{equation} h_{n+1}\int^{a(t)}_0{\phi\,\mathrm{d} x} = h_n a_n \phi_n. \end{equation}

The boundary conditions are given by symmetry at $x=0$, imposing $\partial \phi / \partial x = 0$ at $x=0$ and by requiring $\sigma _{xx}=p=0$ at $x=a(t)$. This latter condition is the same as that applied throughout the entire gel in § 5.1, and so sets the interfacial polymer fraction $\phi$ to equal the steady-state polymer fraction $\phi _{n+1}$ by equating $\varPi (\phi )$ with $4\mu _s(\phi )\epsilon _{zz}$, a condition which arises from subtracting (5.1b) from (5.1a), as illustrated in (5.7).

The initial conditions are simply $\phi = \phi _n$ with $a(0) = a_n' = a_n(1+\varepsilon ) + O(\varepsilon ^2)$. This system of equations can be solved to find the evolution of the polymer fraction in time, and, therefore, also the frontal position $a(t)$ and the force per unit area on the top plate using (5.4), as shown in the following.

5.3.1. Non-dimensional system

In order to solve this problem, we first define non-dimensional variables and parameters

(5.10a–e)\begin{equation} X = \frac{x}{a_n};\quad \varPhi = \frac{\phi}{\phi_n};\quad T = \frac{\hat{k} \hat{K} \phi_n}{\phi_0 a_0^2 \mu_l}t;\quad\mathcal{M} = \left(\frac{\hat{\phi}}{\phi_0}\right)^{1/2}\frac{\widehat{\mu_s}}{\hat{K}};\quad A(T) = \frac{a(t)}{a_n}, \end{equation}

that reduce the diffusion equation (5.8) to

(5.11)\begin{equation} \frac{\partial\varPhi}{\partial T} = \left(1+\mathcal{M}\right)\frac{\partial^2\varPhi}{\partial X^2}. \end{equation}

This is to be solved with the boundary conditions and initial conditions

(5.12)\begin{equation} \left.\begin{gathered} \frac{\partial \varPhi}{\partial X} = 0 \ \text{at $X=0$} \quad{\rm and}\quad \varPhi = \frac{\phi_{n+1}}{\phi_n} \ \text{at $X=A(T)$,} \\ \varPhi(X,0) = 1 \quad{\rm and}\quad A(0) = 1+ \varepsilon, \end{gathered}\right\} \end{equation}

whereas the horizontal extent $A(T)$ is determined from

(5.13)\begin{equation} \int_0^{A(T)}{\varPhi\,\mathrm{d}X} = 1 + \varepsilon. \end{equation}

We can differentiate this expression with respect to time and substitute in the expression for $\partial \varPhi /\partial T$ in (5.11) to find the explicit evolution equation for $A(T)$,

(5.14)\begin{equation} \frac{\phi_{n+1}}{\phi_n}\frac{\partial A}{\partial T} =-\left(1+\mathcal{M}\right)\left.\frac{\partial \varPhi}{\partial X}\right|_{X=A(T)}. \end{equation}

We solved this problem numerically, as detailed in Appendix D and illustrated in figure 4, to determine characteristic behaviour that agrees qualitatively with experiments carried out by Li et al. (Reference Li, Hu, Vlassak and Suo2012). The stress on the top plate is seen to ‘spike’ initially under the incompressible deformation, and then to relax in time to a constant value as water is expelled. Fitting the theoretical results obtained to experimental observations allows us to deduce a value for $\hat {k}$ because all of the other parameters for the problem are known.

Repeating the experiment described previously, reducing the separation distance $h_n$ incrementally from an original value of $h_0$, it is possible to determine constitutive relations for each of the three required material parameters $\varPi (\phi )$, $\mu _s(\phi )$ and $k(\phi )$. The two-dimensional rheometer described here needs elaboration for a three-dimensional realisation, but the principles illustrated are that $\mu _s$ is determined from short-time mechanical responses, $\varPi$ is determined from long-time equilibria and $k$ is determined from transient responses. Importantly, constitutive relationships for all three material parameters can be determined from macroscopic measurements.

6.1. Time evolution of the system

The swelling of the gel is driven by the flow of interstitial water through the hydrogel matrix, governed by gradients in pervadic pressure. As with the rheometer, this swelling behaviour is described by (4.16) with $n=2$, $\boldsymbol {q} = \boldsymbol {0}$ and only derivatives in the $z$ direction, giving

(6.6)\begin{equation} \frac{\partial\phi}{\partial t} = \frac{\partial }{\partial z}\left\lbrace \frac{k(\phi)}{\mu_l}\left[K(\phi) + \mu_s(\phi)\left(\frac{\phi}{\phi_0}\right)^{1/2}\right] \frac{\partial \phi}{\partial z}\right\rbrace. \end{equation}

This equation describes the most general time evolution of such a swelling problem irrespective of the constitutive relations for $k(\phi )$, $K(\phi )$ and $\mu _s(\phi )$.

The position of the swelling front $z=a(t)$ is determined by polymer conservation

(6.7)\begin{equation} \int^{a(t)}_{-\infty}{\left(\phi-\phi^*\right)\mathrm{d}z} + \phi^* a(t) = 0. \end{equation}

6.2. Linearisation

In order for ${\boldsymbol {\epsilon }}$ to remain small throughout the swelling process, $\phi$ must be close to $\phi ^*$ at all times and everywhere in the gel. Therefore, we linearise around the state $\phi = \phi ^*$ and (6.6) becomes a simple linear diffusion equation with diffusion coefficient $D$ given by

(6.8)\begin{equation} D = \frac{k(\phi^*)}{\mu_l}\left[K(\phi^*) + \mu_s(\phi^*)\left(\frac{\phi^*}{\phi_0}\right)^{1/2}\right]. \end{equation}

The interfacial polymer fraction is set explicitly by linearising the condition (6.4), using $\varPi (\phi ) \approx \varPi (\phi ^*) + K(\phi ^*)(\phi -\phi ^*)/\phi ^*$,

(6.9)\begin{equation} \phi_1 = \left[1-\frac{\varPi(\phi^*)}{K(\phi^*) + \mu_s(\phi^*)\left(\phi^*/\phi_0\right)^{1/2}}\right]\phi^* = \left(1-\mathcal{P}\right)\phi^* \end{equation}

with the non-dimensional parameter $\mathcal {P}$ introduced for brevity, and defined to be $\varPi /[K+\mu _s(\phi ^*/\phi _0)^{1/2}]$, evaluated at $\phi = \phi ^*$. This definition implies that $0 < \mathcal {P} < 1$ since $\phi _1$ must take a positive value less than $\phi ^*$. Furthermore, requiring $\varPi / \mu _s \ll 1$ for the assumption of small deviatoric strains to hold enforces $0 < \mathcal {P} \ll 1$. Differentiating the polymer conservation condition (6.7) with respect to time to find an evolution equation for the frontal condition then gives the full set of equations

(6.10a)$$\begin{gather} \frac{\partial \phi}{\partial t} = D\frac{\partial^2 \phi}{\partial z^2} \quad \text{with}\ \phi \to \phi^* \ \text{as $z\to -\infty$} \end{gather}$$

and

(6.10b)$$\begin{gather}\phi = \phi_1 = \left(1-\mathcal{P}\right)\phi^* \quad \text{at $z=a(t)$} \end{gather}$$

which evolves as

(6.10c)$$\begin{gather}\phi_1 \frac{{\rm d}a}{{\rm d} t} =-D\left.\frac{\partial \phi}{\partial z}\right|_{z=a(t)}. \end{gather}$$

This set of equations is entirely analogous to the Stefan problem considering the solidification of a pure melt at a boundary (see, for example, Langer Reference Langer1980; Worster Reference Worster2000; Davis Reference Davis2001). The solution to this swelling problem is

(6.11)\begin{equation} \phi\left(z,t\right) = \phi^*\left[1-\mathcal{P}\frac{\operatorname{erfc}\left(-\eta\right)}{\operatorname{erfc} \left(-\lambda\right)}\right] \quad \text{where} \ \eta = \frac{z}{2\sqrt{Dt}}, \end{equation}

and $a(t) = 2\lambda \sqrt {Dt}$, with the parameter $\lambda$ set implicitly by

(6.12)\begin{equation} \sqrt{\rm \pi}\lambda e^{\lambda^2} \operatorname{erfc}\left(-\lambda\right) = \frac{\mathcal{P}}{1-\mathcal{P}}. \end{equation}

Figure 6 shows how the parameter $\lambda$ depends on the value of $\mathcal {P}$, showing zero growth when $\mathcal {P} = 0$, corresponding to $\phi ^* = \phi _0$, whilst growth rates increase to infinity as $\mathcal {P} \to 1$ and the gel swells to an ever greater degree (i.e. $\phi _1$ decreases). Of course, $\mathcal {P}$ must remain small to satisfy the assumption of small deviatoric strains, and in this case (6.12) can be approximated by $\lambda \approx (\mathcal {P}/\sqrt {{\rm \pi} })[1+(1-2/{\rm \pi} )\mathcal {P} + \dots ]$, as illustrated in figure 6.

Figure 6.

The growth rate $\lambda$ for the linearised planar growth problem, as defined by the implicit relation (6.12), showing how the position of the hydrogel front $a(t)$ grows more rapidly for larger values of $\mathcal {P}$. The small-$\mathcal {P}$ limit is illustrated with a dashed black line.

7.1. Post-hoc justification of small deviatoric strains

Our model assumes from the outset that deviatoric strains remain small at all times, and we can linearise around a state of isotropic swelling. It is straightforward to find conditions on the parameters of the model described by (7.8) under which these strains will indeed remain small, and to check whether our assumption that large isotropic volume changes can occur without ever introducing significant deviatoric strains is borne out in reality. As this model is one-dimensional, it is straightforward to solve for the radial displacement, which must satisfy

(7.9)\begin{equation} \frac{\partial \xi}{\partial r} + \frac{2\xi}{r} = 3\left[1-\left(\frac{\phi}{\phi_0}\right)^{1/3}\right] \quad {\rm and}\quad \xi = 0 \quad \text{at $r=0$}, \end{equation}

from which we can calculate the three components of ${\boldsymbol {\epsilon }}$,

(7.10a,b)\begin{equation} \epsilon_{rr} = 2\left[1-\frac{\xi}{r}-\left(\frac{\phi}{\phi_0}\right)^{1/3}\right];\quad \epsilon_{\theta\theta}=\epsilon_{\varphi\varphi} = \left[\frac{\xi}{r}+\left(\frac{\phi}{\phi_0}\right)^{1/3}\right]. \end{equation}

The plots in figure 9 show that the deviatoric strain is greatest at the interface and, therefore, the deviatoric strains are small in our non-dimensional model provided that

(7.11)\begin{equation} 2\left(A(T)^{-1}-\varPhi_1^{1/3}\right) \ll 1 \quad \text{or} \quad \frac{\left(\phi_1/\phi_0\right) - 1}{2\mu_s/K_0} \ll 1, \end{equation}

through the use of the expression (7.8b) for the interfacial polymer fraction. We see the traditional linear case arises when $\varPhi _1 - 1 \ll 1$, but note here that our approach is still valid with $\varPhi _1 - 1$ finite provided that $\mathcal {M} = \mu _s / K_0$ is large.

Figure 9.

The deviatoric strain $\epsilon _{rr}$ (the component of ${\boldsymbol {\epsilon }}$ with the largest magnitude) plotted against non-dimensional radius at various times $0 \leqslant T \leqslant 1$ when $\varPhi ^* = 8$ and $\mathcal {M}=40$; in this case, the initial deviatoric strains are shown in black. In spite of the gel swelling to eight times its initial volume, deviatoric strains stay below $10\,\%$ in magnitude throughout, the approximate threshold quoted by Landau & Lifshitz (Reference Landau and Lifshitz1986) for linear theory to hold.