Flow and deformation due to periodic loading in a soft porous material

Abstract

Soft porous materials, such as biological tissues and soils, are exposed to periodic deformations in a variety of natural and industrial contexts. The detailed flow and mechanics of these deformations have not yet been systematically investigated. Here, we fill this gap by identifying and exploring the complete parameter space associated with periodic deformations in the context of a one-dimensional model problem. We use large-deformation poroelasticity to consider a wide range of loading periods and amplitudes. We identify two distinct mechanical regimes, distinguished by whether the loading period is slow or fast relative to the characteristic poroelastic timescale. We develop analytical solutions for slow loading at any amplitude and for infinitesimal amplitude at any period. We use these analytical solutions and a full numerical solution to explore the localisation of the deformation near the permeable boundary as the period decreases and the emergence of nonlinear effects as the amplitude increases. We show that large deformations lead to asymmetry between the loading and unloading phases of each cycle in terms of the distributions of porosity and fluid flux.

1. Introduction

Soft porous materials are common in nature and industry; examples include biological cells and soft tissues, soils and sediments, and paper products and fabrics. In many scenarios, these materials are exposed to periodic loading. For example, soft tissues in the body can experience pulsating loads from the surrounding blood vessels or, on a larger scale, can be cyclically loaded during their basic mechanical function. The former scenario has attracted great interest recently as a potential driver of transport in brain tissue (Franceschini et al. 2006; Kedarasetti, Drew & Costanzo 2020; Bojarskaite et al. 2023) and the latter is important for load-bearing and transport in cartilage (Riches et al. 2002; Mauck, Hung & Ateshian 2003; Ferguson, Ito & Pyrak-Nolte 2004; Sengers, Oomens & Baaijens 2004; Schmidt et al. 2010; Zhang 2011; Di Domenico, Wang & Bonassar 2017; Cacheux et al. 2023) and bone (Piekarski & Munro 1977; Zhang & Cowin 1994; Manfredini et al. 1999; Nguyen, Lemaire & Naili 2010; Witt et al. 2014). Periodic loads are also commonly applied in regenerative medicine to improve cell differentiation in scaffolds via mechanotransduction (Kim et al. 1999; Butler, Goldstein & Guilak 2000; Mauck et al. 2000; Grenier et al. 2005; Haj, Hampson & Gogniat 2009; Gauvin et al. 2011; Amrollahi & Tayebi 2015; Peroglio et al. 2018). Soils and sediments experience periodic loading due to seismicity (Genna & Cividini 1989; Li, Borja & Regueiro 2004; Popescu et al. 2006; Bonazzi, Jha & de Barros 2021), vehicle traffic (Hu, Zhu & Cheng 2011; Ni et al. 2015; Ni & Geng 2022) and ocean waves and tides (Madsen 1978; Yamamoto et al. 1978; Karim, Nogami & Wang 2002; Cheng 2016; Trefry et al. 2019).

From a poromechanical point of view, periodic loading is fundamentally different from steady loading because the long-time response is inherently oscillatory and therefore time-dependent. Most previous work on periodic loading has focused on internal stress and/or pressure profiles, on macroscopic observables such as surface motion or net inflow or outflow, or on solute concentration profiles.

Periodic loading due to seismicity is a classical topic in poroelasticity (Biot 1956a,b). Seismicity involves frequencies that are high enough for inertia to play a dominant role. As a result, seismic response is typically dominated by the propagation of compressional and shear acoustic waves (e.g. Biot 1956a,b; Li et al. 2004; Gajo & Denzer 2011; Liu et al. 2019). In contrast, ocean waves and tides are typically associated with a low enough frequency that inertia can be ignored (e.g. Madsen 1978; Yamamoto et al. 1978; Cheng 2016). Instead, these studies typically focus on the pressure and stress profiles within seabed sediments in response to periodic fluctuations in hydrostatic pressure. The sediment is usually taken to be semi-infinite, the associated deformations are assumed to be small, and the response is often dominated by compressibility. Tidal forcing in coastal aquifers often has similar features (e.g. Trefry et al. 2019).

Our primary motivation here is tissue mechanics, where inertia and compressibility are usually negligible but moderate to large deformations are common. In this regime, poroelasticity is physically diffusive with a characteristic poroelastic relaxation time. For bone and cartilage, linear poroelasticity has been used to model the macroscopic mechanical response and/or the distribution of pore pressure during small periodic deformations (e.g. Zhang & Cowin 1994; Manfredini et al. 1999; Riches et al. 2002; Kameo, Adachi & Hojo 2008; Yaogeng et al. 2018). Zhang & Cowin (1994) showed that the magnitude and distribution of pore pressure depend strongly on the loading period, introducing the ratio of the loading period to the poroelastic relaxation time as a key dimensionless control parameter. For cartilage and hydrogel scaffolds, both linear and nonlinear poroelasticity have been used to model the impact of periodic deformations on the transport of solutes, typically by comparing the concentration profile at the end of loading across a small set of different loading conditions (e.g. Mauck et al. 2003; Ferguson et al. 2004; Sengers et al. 2004; Gardiner et al. 2007; Urciuolo, Imparato & Netti 2008; Zhang 2011; Vaughan et al. 2013). Several of these studies noted that faster loading (shorter loading period) is associated with larger fluid velocities that are localised near the surface, whereas slower loading (longer loading period) is associated with intermediate fluid velocities that penetrate more deeply (Gardiner et al. 2007; Kameo et al. 2008; Urciuolo et al. 2008; Di Domenico et al. 2017; Vaughan et al. 2013). Gardiner et al. (2007) further noted that, for small deformations, the magnitude of the fluid velocity is roughly proportional to the loading amplitude, whereas the penetration depth is relatively insensitive to amplitude. Despite this extensive previous work, many basic features of flow and deformation due to periodic loading have not yet been systematically studied, in part because most of the above studies have focused on relatively narrow regions of the associated parameter space and/or on relatively small sets of specific results. For example, kinematic and constitutive nonlinearities (characteristic of large deformations and nonlinear constitutive behaviour, respectively) become increasingly important as the amplitude grows, but the emergence of these nonlinearities has not been investigated. Moreover, the impacts of the deformation on the magnitude and distribution of the fluid flux, directly relevant to the transport of solutes, have not been examined.

Here, we study the periodic loading of a soft porous material over a wide range of loading periods (from very slow to very fast) and amplitudes (from infinitesimal to moderate/large) in the context of a simple one-dimensional model problem. Following MacMinn, Dufresne & Wettlaufer (2016), our large-deformation poroelastic model is kinematically rigorous and includes both deformation-dependent permeability and nonlinear elasticity. We characterise the motion of the fluid and the solid throughout the loading cycle, focusing on the evolution of porosity and fluid flux, which are particularly relevant for the transport of solutes. We develop a series of analytical solutions that describe loading with small amplitude but arbitrary period (§§ 3.2–3.4) and loading with large period but arbitrary amplitude (§ 3.5). We then solve the full problem numerically and compare with our analytical results. We use these solutions to examine the transition from very slow loading to very fast loading and from infinitesimal to large amplitudes, as well as the role of the initial porosity. Finally, we discuss the relevance of these results to some specific biological and biomedical scenarios.

2. Theoretical model

Our theoretical model is based on large-deformation poroelasticity (also referred to as biphasic theory in biomedical communities), here following the development in MacMinn et al. (2016).

2.1. Model problem

We consider a one-dimensional sample of soft porous material of relaxed length $L$ and relaxed porosity (fluid fraction) $\phi _{f,0}$. The left boundary of the material is permeable and located at $x=a(t)$ (moving); the right boundary is impermeable and located at $x=L$ (fixed in place) (figure 1). We impose a periodic, displacement-driven loading via the position of the left boundary,

(2.1)\begin{equation} a(t)= \frac{A}{2} \left[1-\cos\left(\frac{2{\rm \pi} t}{T}\right)\right], \end{equation}

where $A$ and $T$ are the amplitude and period of the loading, respectively. Macroscopically, the deformation is strictly compressive in the sense that $a(t)\geq {}0$.

Figure 1. We consider a one-dimensional sample of soft porous material of relaxed length $L$, subject to a periodic, displacement-driven loading at its left boundary (white arrows). The left boundary is permeable, thus allowing fluid flow in or out (blue squiggles) to accommodate the loading. The right boundary is impermeable and fixed in place.

We consider deformations ranging from small to large macroscopic nominal strain ($-0.4\,\%$ to $-20\,\%$ or $0.004\leq {}A/L\leq {}0.2$), with commensurate macroscopic changes in bulk (total) volume. We assume that the fluid and the solid phases are individually incompressible, so that the total volume of solid is constant and any change in bulk volume must correspond to a change in total pore (fluid) volume via a rearrangement of the pore structure and an influx or an efflux of fluid at the left boundary.

2.2. Kinematics

We work in an Eulerian reference frame, in which the solid displacement field is $\boldsymbol {u_s}=\boldsymbol {x}-\boldsymbol {X}(\boldsymbol {x},t)$, with $\boldsymbol {X}(\boldsymbol {x},t)$ the reference position of the material point that at time $t$ occupies the position $\boldsymbol {x}$. We choose $\boldsymbol {X}(\boldsymbol {x},0)\equiv \boldsymbol {x}$ so that $\boldsymbol {u_{s}}(\boldsymbol {x},0)=0$, in which case the reference configuration is the relaxed configuration. We denote the true volume fractions of fluid and solid by $\phi _f$ and $\phi _{s}$, respectively, where $\phi _{f}+\phi _{s}=1$. The flow and deformation are uniaxial, such that

(2.2a–c)\begin{equation} \boldsymbol{u_s} = u_s(x,t)\boldsymbol{\hat{e}_x}, \quad \boldsymbol{v_s} = v_s(x,t)\boldsymbol{\hat{e}_x}, \quad \boldsymbol{v_f} = v_f(x,t)\boldsymbol{\hat{e}_x}, \end{equation}

where $u_s$, $v_s$ and $v_f$ are the $x$-components of the solid displacement field and the solid and fluid velocity fields, respectively, and $\boldsymbol {\hat {e}_x}$ is the unit vector in the $x$-direction.

In one dimension, the local state of deformation is fully characterised by the Jacobian determinant $J= (1-\partial {u_s}/\partial {x})^{-1}$, which measures the local current volume per unit reference volume. For incompressible constituents and uniform initial porosity $\phi _{f,0}$, the local change in volume is linked to the change in porosity according to

(2.3)\begin{equation} J(x,t) = \frac{1-\phi_{f,0}}{1-\phi_f} \quad\to\quad \frac{\partial{u_s}}{\partial{x}} =\frac{\phi_f-\phi_{f,0}}{1-\phi_{f,0}}. \end{equation}

Continuity for this one-dimensional system can be written

(2.4a,b)\begin{equation} \frac{\partial{\phi_f}}{\partial{t}} +\frac{\partial}{\partial{x}}(\phi_f v_f) = 0 \quad\mathrm{and}\quad \frac{\partial{\phi_f}}{\partial{t}} -\frac{\partial}{\partial{x}}{[(1-\phi_f)v_s]} = 0, \end{equation}

which together imply that the total flux $q=\phi _fv_f+(1-\phi _f)v_s$ is uniform in space, $\partial {q}/\partial {x}=0$.

2.3. Darcy's law

We assume that the movement of the fluid relative to the solid skeleton is described by Darcy's law,

(2.5)\begin{equation} \phi_f(v_f-v_s) ={-}\frac{k(\phi_f)}{\mu}\frac{\partial{p}}{\partial{x}}, \end{equation}

where $k(\phi _f)$ is the permeability of the solid skeleton, $\mu$ is the dynamic viscosity of the fluid and $p$ is the fluid (pore) pressure, and where we have neglected gravity. We have taken the permeability $k(\phi _f)$ to be a function of porosity only. For simplicity, we use a normalised Kozeny–Carman relation for deformation-dependent permeability (MacMinn et al. 2016):

(2.6)\begin{equation} k(\phi_f) = k_0\frac{(1-\phi_{f,0})^2}{\phi_{f,0}^3} \frac{\phi_f^3}{(1-\phi_f)^2}, \end{equation}

where $k_0\equiv {}k(\phi _{f,0})$ is the reference permeability. Kozeny–Carman permeability is a common choice for gels and soft tissues (Malandrino et al. 2014; Sacco et al. 2014; Rahbari et al. 2017; Gao & Cho 2022). We compare it with a simpler power-law formulation in Appendix A, showing that the two are qualitatively and quantitatively similar and are expected to produce similar behaviour.

2.4. Fluid flow

Combining (2.4a,b) and (2.5), we arrive at the nonlinear flow equations:

(2.7a,b)\begin{equation} \frac{\partial{\phi_f}}{\partial{t}} +\frac{\partial}{\partial{x}}\left[{\phi_f q}-(1-\phi_f) \frac{k(\phi_f)}{\mu}\frac{\partial{p}}{\partial{x}}\right]=0 \quad\mathrm{and}\quad \frac{\partial{q}}{\partial{x}}=0, \end{equation}

where the total flux $q$ is again

(2.8)\begin{equation} q\equiv\phi_f v_f + (1-\phi_f)v_s. \end{equation}

The Darcy-scale fluid velocity and solid velocity are then given by

(2.9a,b)\begin{equation} v_f =q-\frac{1-\phi_f}{\phi_f}\,\frac{k(\phi_f)}{\mu}\frac{\partial{p}}{\partial{x}} \quad\mathrm{and}\quad v_s=q+\frac{k(\phi_f)}{\mu}\frac{\partial{p}}{\partial{x}}. \end{equation}

and the local fluid flux is

(2.10)\begin{equation} q_f=\phi_f v_f. \end{equation}

2.5. Mechanical equilibrium

The true Cauchy total stress $\boldsymbol {\sigma }$ is supported jointly by the fluid phase and the solid phase. The total stress can be decomposed into a contribution from the fluid pressure $p$ and a contribution from Terzaghi's effective stress $\boldsymbol {\sigma '}$,

(2.11)\begin{equation} \boldsymbol{\sigma}=\boldsymbol{\sigma'}-p\boldsymbol{I}, \end{equation}

where we adopt the sign convention of tension being positive. Neglecting inertia and body forces, mechanical equilibrium can be written

(2.12)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma} = \boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{\sigma'}-\boldsymbol{\nabla}p=0. \end{equation}

In one dimension, (2.12) implies that

(2.13)\begin{equation} \frac{\partial{\sigma'}}{\partial x}=\frac{\partial{p}}{\partial{x}}, \end{equation}

where $\sigma ^\prime$ is the $xx$ component of $\boldsymbol {\sigma }^\prime$.

2.6. Elasticity law

The effective stress is the portion of the total stress that contributes to deformation of the solid skeleton. We take the solid skeleton to be elastic, with no viscous or dissipative behaviours. For confined compression in one dimension, as considered here, any elasticity law can be written in the form $\sigma ^\prime =\sigma ^\prime (\phi _f)$. Thus, (2.7a,b) can be rewritten as a nonlinear advection–diffusion equation:

(2.14a,b)\begin{equation} \frac{\partial{\phi_f}}{\partial{t}} +\frac{\partial}{\partial{x}}\left[{\phi_f q}-D_f(\phi_f)\frac{\partial{\phi_f}}{\partial{x}}\right]=0 \quad\mathrm{and}\quad \frac{\partial{q}}{\partial{x}}=0, \end{equation}

where the nonlinear composite constitutive function

(2.15)\begin{equation} D_f(\phi_f)=(1-\phi_f)\frac{k(\phi_f)}{\mu}\frac{\mathrm{d}\sigma^\prime}{\mathrm{d}\phi_f} \end{equation}

is the poroelastic diffusivity.

Hencky elasticity is a simple nonlinear hyperelasticity model that considers logarithmic strains (‘true strains’ or ‘Hencky strains’) to capture kinematic nonlinearity (Hencky 1931), and which is commonly used to model soft rubbers and foams (Hencky 1933; Anand 1979; Xiao & Chen 2002) and sometimes for soft biological tissues (Marchesseau et al. 2010; Fraldi et al. 2019). Hencky elasticity is convenient for our present purposes because (i) it reduces to linear elasticity for small strains and (ii) it uses the same two elastic parameters as linear elasticity. It is straightforward to replace Hencky elasticity in the present formulation with a different elastic behaviour, such as a neo-Hookean model, as appropriate for the problem/material of interest. Neo-Hookean elasticity is commonly used as a simple model for soft tissues (e.g. Sengers et al. 2004; Ehlers, Karajan & Markert 2009); in the present context, we expect Hencky elasticity to provide a qualitatively similar mechanical response (see Appendix A).

For a uniaxial deformation, the relevant component of the effective stress tensor for Hencky elasticity is (MacMinn et al. 2016; Auton & MacMinn 2018)

(2.16)\begin{equation} \sigma^\prime= \mathcal{M}\frac{\ln(J)}{J}=\mathcal{M} \left(\frac{1-\phi_f}{1-\phi_{f,0}}\right) \ln\left(\frac{1-\phi_{f,0}}{1-\phi_f}\right), \end{equation}

where $\mathcal {M}$ is the $p$-wave or oedometric modulus. With appropriate initial and boundary conditions, (2.14a,b), (2.15) and (2.16) form a closed model for the evolution of the porosity.

2.7. Initial and boundary conditions

Finally, we specify appropriate initial and boundary conditions for the solid skeleton and for the fluid. As noted previously, we locate the left and right boundaries of the solid at $x=a(t)$ and $x=L$, respectively.

2.7.1. Initial conditions

Equation (2.1) suggests that $a(0)=0$. The solid is therefore initially relaxed and the initial porosity is uniform and equal to the relaxed porosity,

(2.17a,b)\begin{equation} u_s(x,0)=0 \quad\mathrm{and}\quad \phi_f(x,0)=\phi_{f,0}. \end{equation}

2.7.2. Left boundary

For $t>0$, we apply a displacement-controlled mechanical loading at the left boundary, which is therefore a moving boundary (see (2.1)). The associated boundary conditions are

(2.18a)\begin{equation} u_s(a,t)=a(t) \quad\mathrm{and}\quad v_s(a,t)=\frac{\mathrm{d}a}{\mathrm{d}t}. \end{equation}

The left boundary is also permeable, so we take

(2.18b)\begin{equation} p(a,t)=0. \end{equation}

2.7.3. Right boundary

The right boundary is impermeable and fixed in place, such that

(2.19)\begin{equation} u_s(L,t)=v_s(L,t)=v_f(L,t)=0. \end{equation}

This condition and the requirement that $q$ be uniform in space (see the end of § 2.2 and (2.7a,b)) together imply that $q\equiv 0$, meaning that there is no net flow through any cross-section. Equation (2.8) then requires that

(2.20)\begin{equation} v_f={-}\frac{(1-\phi_f)}{\phi_f} v_s, \end{equation}

meaning that the fluid and the solid always locally move in opposite directions.

2.8. Linear poroelasticity

For comparison with the fully nonlinear model, we linearise the relations given previously to arrive at linear poroelasticity, which is valid for infinitesimal deformations, $(\phi _f-\phi _{f,0})/(1-\phi _{f,0})=\partial {u_s}/\partial {x}\ll {}1$. In this limit, (2.7a,b) reduces to the linear-poroelastic diffusion equation,

(2.21)\begin{equation} \frac{\partial{\phi_f}}{\partial{t}} -\frac{\partial}{\partial{x}}\left(D_{f,0}\frac{\partial{\phi_f}}{\partial{x}}\right)\approx 0, \end{equation}

where Hencky elasticity reduces to linear elasticity,

(2.22)\begin{equation} \sigma^\prime\approx \mathcal{M}\,\frac{\partial{u_s}}{\partial{x}}=\mathcal{M}\left(\frac{\phi_f-\phi_{f,0}}{1-\phi_{f,0}}\right), \end{equation}

and $D_{f,0}=k_0\mathcal {M}/\mu$ is the constant linear-poroelastic diffusivity. With appropriate initial and boundary conditions, (2.21) is a closed linear model for the evolution of the porosity.

In the linear poroelastic model, the initial conditions and the boundary conditions for the right boundary are again (2.17a,b) and (2.19), respectively. The linearised boundary conditions for the left boundary are

(2.23a,b)\begin{equation} u_s(0,t)\approx a(t) \quad\mathrm{and}\quad v_s(0,t) \approx \frac{\mathrm{d}a}{\mathrm{d}t}, \end{equation}

where $a(t)$ is again given by (2.1). Note that these conditions are linearised relative to (2.18) by virtue of being applied at $x=0$ rather than at $x=a(t)$.

2.9. Scaling and summary

We make the above problem dimensionless via the scaling

(2.24a–h)\begin{align} \tilde{x}=\frac{x}{L},\quad \tilde{u}_s=\frac{u_s}{L},\quad \tilde{t}=\frac{t}{T_{{pe}}},\quad \tilde{\sigma}^\prime=\frac{\sigma'}{\mathcal{M}},\quad \tilde{p}=\frac{p}{\mathcal{M}},\quad \tilde{k}=\frac{k(\phi)}{k_0},\; \tilde{v}_f=\frac{v_f}{L/T_{{pe}}}, \quad \tilde{v}_s=\frac{v_s}{L/T_{{pe}}}, \end{align}

where $T_{{pe}}=L^2/D_{f,0}=\mu {}L^2/(k_0\mathcal {M})$ is the classical poroelastic timescale, which is the characteristic diffusion time for the relaxation of pressure over a distance $L$. Now taking $q\equiv {}0$, as required by the boundary conditions (see § 2.7.3), the nonlinear flow equation can be rewritten in dimensionless form as

(2.25)\begin{equation} \frac{\partial{\phi_f}}{\partial{\tilde{t}}} -\frac{\partial}{\partial{\tilde{x}}}\left[\tilde{D}_f(\phi_f)\frac{\partial{\phi_f}}{\partial{\tilde{x}}}\right]=0 \end{equation}

with nonlinear-poroelastic diffusivity

(2.26)\begin{equation} \tilde{D}_f=\frac{D_f}{D_{f,0}}=(1-\phi_f)\tilde{k}(\phi_f)\frac{\mathrm{d}\tilde{\sigma}^\prime}{\mathrm{d}\phi_f}, \end{equation}

elasticity law

(2.27)\begin{equation} \tilde{\sigma}^\prime=\left(\frac{1-\phi_f}{1-\phi_{f,0}}\right) \ln\left(\frac{1-\phi_{f,0}}{1-\phi_f}\right), \end{equation}

initial conditions

(2.28a–c)\begin{equation} \tilde{a}(0)=0, \quad \phi_f(\tilde{x},0)= \phi_{f,0}, \quad \tilde{v}_f(\tilde{x},0)=\tilde{v}_s(\tilde{x},0)=0, \end{equation}

left boundary conditions

(2.29a,b)\begin{equation} \tilde{u}_s(\tilde{a},\tilde{t})=\tilde{a}(\tilde{t})= \frac{\tilde{A}}{2} \left[1-\cos\left(\frac{2{\rm \pi}\tilde{t}}{\tilde{T}}\right)\right] ,\quad \tilde{v}_s(\tilde{a},\tilde{t})=\frac{\mathrm{d} \tilde{a}}{\mathrm{d}\tilde{t}}\quad \mathrm{and}\quad \tilde{p}(\tilde{a},\tilde{t})=0, \end{equation}

and right boundary conditions

(2.30)\begin{equation} \tilde{u}_s(1,\tilde{t}) =\tilde{v}_s(1,\tilde{t}) =\tilde{v}_f(1,\tilde{t})=0, \end{equation}

where $\tilde {A}=A/L$ and $\tilde {T}=T/T_{{pe}}$. We consider only dimensionless quantities in the following, dropping the tildes for convenience.

This one-dimensional model describes flow and mechanics in a poroelastic material subject to periodic loading. The kinematics are rigorous and nonlinear, the elasticity law is Hencky elasticity, and the permeability law is the Kozeny–Carman relation. The full and linearised problems share the same three dimensionless control parameters: the dimensionless amplitude and period of the loading, $A$ and $T$, and the relaxed porosity, $\phi _{f,0}$.

3. Analytical solutions

We next develop three different analytical solutions to the linear-poroelastic problem, which are valid for small deformations ($A\ll {}1$), and to the full problem for slow deformations ($T\gg {}1$) at any amplitude (i.e. the quasi-static limit). As formulated in § 2.8, the linear-poroelastic problem implies linearised kinematics, linear elasticity and constant permeability and, thus, a constant and uniform poroelastic diffusivity. We also solve the full problem numerically in general.

3.1. Average porosity

We begin by deriving some basic kinematic results for the average porosity. The macroscopic total volume at any instant is $1-a(t)$, whereas the total volume of solid is constant and equal to $1-\phi _{f,0}$. As a result, the total volume of fluid is $\phi _{f,0}-a(t)$ and the spatially averaged porosity is

(3.1)\begin{equation} \langle \phi_f \rangle (t) = \frac{\phi_{f,0}-a(t)}{1-a(t)}. \end{equation}

The average of $\langle \phi _f\rangle (t)$ in time over any integer number of loading cycles is then

(3.2)\begin{equation} \overline{\langle\phi_f\rangle} = \frac{1}{mT}\int_{nT}^{(n+m)T} \langle{\phi_f}\rangle \,\mathrm{d}t = 1 - \frac{1-\phi_{f,0}}{\sqrt{1-A}} \end{equation}

for any $n\geq 0$ and integer $m\geq 1$. Note that both the spatial and overall averages are smaller than $\phi _{f,0}$ for $a>0$ and $A>0$ because the loading has a non-zero mean (i.e. $\bar {a}=A/2>0$), so the material is on average compressed.

3.2. Linear poroelasticity: early time solution

The linear problem posed in § 2.8 can be rewritten as a bounded linear diffusion problem for the displacement. When the loading begins, information about the motion of the left boundary propagates into the domain via poroelastic diffusion. At early times, before this information has had time to reach the right boundary, the response is the same as if the material were semi-infinite in the $x$ direction. The corresponding semi-infinite diffusion problem involves applying the right boundary conditions at $x\to \infty$. This early time (‘$\mathrm {et}$’) solution can be derived via Laplace transform and written as a convolution integral,

(3.3)\begin{equation} u_{s,{et}}(x,t)=\frac{{\rm \pi}{}A}{T}\int_0^t\,\mathrm{erfc}\left(\frac{x}{2\sqrt{\tau}}\right)\sin\left[\frac{2{\rm \pi}(t-\tau)}{T}\right]\,\mathrm{d}\tau. \end{equation}

The corresponding porosity field can be derived from (3.3) via (2.3), and is given by

(3.4)\begin{equation} {\phi_{f,{et}}}(x,t)=\phi_{f,0}-(1-\phi_{f,0})\frac{{\rm \pi}{}A}{T}\int_0^t\,\frac{1}{\sqrt{{\rm \pi}\tau}}\exp\left(-\frac{x^2}{4\tau}\right)\sin\left[\frac{2{\rm \pi}(t-\tau)}{T}\right]\,\mathrm{d}\tau. \end{equation}

This solution is valid until the deformation spans the domain. Introducing the penetration length $\delta _{et}(t)$ as the distance from the left boundary over which the change in porosity has exceeded a threshold, the nature of the problem and the structure of the solution suggest that $\delta _{et}(t)\sim {}2\sqrt {t}$. We confirm this reasoning in Appendix B. Thus, the above solution is valid for $\delta _{et}\lesssim {}1\,\to \,t\lesssim {}1/4$.

Having originated from linear poroelasticity, the above solution is limited to small deformations ($A\ll 1$). Naïvely, this solution is valid for any loading period $T$; however, sufficiently small values of $T$ also violate the assumption of small deformations because deformation of size $\sim {}A$ are localised to a region of size $\sim \delta _{et}(t)$ at early times. As a result, we expect the maximum strain near the left boundary to be roughly of size $\max [a(t)/\delta _{et}(t)]\sim {}A/\sqrt {T}$. More precisely, (3.4) can be used to show that the extreme values of ${\phi _{f,{et}}}$ will always occur at the left boundary and that the evolution of ${\phi _{f,{et}}}$ at the left boundary is given by

(3.5)\begin{equation} {\phi_{f,{et}}}(0,t)=\phi_{f,0}-(1-\phi_{f,0})A\sqrt{\frac{\rm \pi}{T}}\, \mathcal{I}(t/T), \end{equation}

where

(3.6)\begin{equation} \mathcal{I}(y)=C(2\sqrt{y})\sin(2{\rm \pi}{}y)-S(2\sqrt{y})\cos(2{\rm \pi}{}y) \end{equation}

and $C$ and $S$ are the Fresnel cosine and sine integrals, respectively. The extreme values of ${\phi _{f,{et}}}$ are then given by

(3.7)\begin{equation} \min({\phi_{f,{et}}}) =\phi_{f,0}-(1-\phi_{f,0})A\sqrt{\frac{\rm \pi}{T}}\,\mathcal{I}_{max}, \end{equation}

and

(3.8)\begin{equation} \max({\phi_{f,{et}}}) =\phi_{f,0}-(1-\phi_{f,0})A\sqrt{\frac{\rm \pi}{T}}\,\mathcal{I}_{min}, \end{equation}

where $\mathcal {I}_{max}=\mathcal {I}(698/1909)\approx {}0.9491$ and $\mathcal {I}_{min}=\mathcal {I}(310/353)\approx {}-\!0.5406$ are the maximum and minimum values of $\mathcal {I}$, respectively. Strictly, the above solution is non-physical for parameter combinations for which the porosity decreases to 0 or increases to 1, corresponding to $\min ({\phi _{f,{et}}})=0$ and $\max ({\phi _{f,{et}}})=1$, respectively. In practice, these solutions become inaccurate for much less extreme parameter combinations as kinematic and constitutive nonlinearities become increasingly important. Hewitt et al. (2016) showed that very fast monotonic compression of a soft porous material can lead to extreme localisation near the piston in the form of a ‘bloated’ low-porosity boundary layer, the formation and evolution of which depends sensitively on the particular constitutive functions $k(\phi _f)$ and $\sigma ^\prime (\phi _f)$. We avoid these extreme loading conditions in the present study, focusing instead on scenarios that are likely to have more biological relevance. The above results confirm that the maximum local strain is indeed of size $A/\sqrt {T}$, suggesting that the above solution and the linear model in general are valid for $A\ll {}\sqrt {T}$. A similar condition can be derived by noting that, in the linear-poroelastic case, the boundary conditions at the left are applied at $x=0$ rather than at $x=a(t)$ (see (2.23a,b)). Thus, the validity of the linear-poroelastic model requires that the error in this linearisation, which is $\sim {}A$, must be negligible relative to the poroelastic diffusion length associated with the deformation, which is $\sim \sqrt {T}$.

3.3. Linear poroelasticity: full solution

The original bounded linear diffusion problem can be solved analytically via separation of variables. The resulting linear-poroelastic (‘$\mathrm {lpe}$’) displacement field is

(3.9)\begin{align} u_{s,{lpe}}(x,t) &= a(t)(1-x) \nonumber\\ &\quad - \sum_{n=1}^{\infty} 2A \sin{(n{\rm \pi} x)} \dfrac{\left[2 e^{{-}n^2{\rm \pi}^2t}-2\cos{\left(\dfrac{2{\rm \pi} t}{T}\right)}+n^2{\rm \pi} T \sin{\left(\dfrac{2{\rm \pi} t}{T}\right)}\right]}{n {\rm \pi}(4+ n^4 {\rm \pi}^2 T^2)}. \end{align}

The corresponding porosity field can be derived from (3.9) via (2.3) and is given by

(3.10)\begin{align} &\phi_{f,{lpe}}(x,t)=\phi_{f,0}- (1-\phi_{f,0})\nonumber\\ &\qquad \times \left\{ a(t) + \sum_{n=1}^{\infty} 2A \cos{(n{\rm \pi} x)} \dfrac{\left[2 e^{{-}n^2{\rm \pi}^2t}-2\cos{\left(\dfrac{2{\rm \pi} t}{T}\right)}+n^2{\rm \pi} T \sin{\left(\dfrac{2{\rm \pi} t}{T}\right)}\right]}{ (4+ n^4 {\rm \pi}^2 T^2)}\right\} . \end{align}

The corresponding fluid velocity field can be derived by taking the time derivative of (3.9) to obtain the solid velocity and then using (2.20) to arrive at

(3.11)\begin{align} &v_{f,{lpe}}(x,t)={-}\dfrac{(1-\phi_{f,0})}{\phi_{f,0}} \nonumber\\ &\qquad \times \left\{ \dot a(t)(1-x) - \sum_{n=1}^{\infty} 4 A \sin{(n{\rm \pi} x)} \dfrac{\left[ - n^2{\rm \pi} e^{{-}n^2{\rm \pi}^2t}+\dfrac{2}{T}\sin{\left(\dfrac{2{\rm \pi} t}{T}\right)}+n^2{\rm \pi} \cos{\left(\dfrac{2{\rm \pi} t}{T}\right)}\right]}{n \left(4+ n^4 {\rm \pi}^2 T^2\right)}\right\}. \end{align}

Like the early time solution, this solution provides a good approximation for $A\ll 1$ and for $A\ll {}\sqrt {T}$. Also like the early time solution, this solution provides general insight into the poromechanical response of the system to periodic loading. The expression for ${\phi _{f,{lpe}}}$ can be divided into three parts:

1. (i) a uniform quasi-static part proportional to $a(t)$, which is the linearised form of the nonlinear quasi-static solution derived in the following;

2. (ii) an early time transient that decays exponentially in time at a rate that is independent of $A$ and $T$; and

3. (iii) a periodic forced response with period $T$.

The early time transient captures the early time solution derived previously, which spans the domain after a time $t\approx {}1/4$ and then decays exponentially relative to the periodic forced response that dominates the solution thereafter. Going forward, we focus on this periodic regime. We show in Appendix C that the same reasoning also applies for large deformations.

3.4. Linear poroelasticity: response to very fast loading (Stokes’ second problem)

As noted previously, the response at early times will be confined to a region of size $\sim \sqrt {t}$, spreading diffusively until the entire domain is engaged, at which point the response will evolve exponentially toward its periodic regime. However, the oscillations in the periodic regime will be confined to a region of size $\sim \sqrt {T}$ (or $\sim {}A$ if larger, but recall that linear poroelasticity requires that $A\ll {}\sqrt {T}$). Thus, if the period is sufficiently small (i.e. $\sqrt {T}\ll 1$), the material near the piston will oscillate while the far field exists in a state of static compression. This response to very fast loading is well known from Stokes's classical ‘second problem’, in which oscillations diffuse into a semi-infinite domain with an amplitude that decays exponentially in space, much like an evanescent wave. The corresponding analytical solution to linear poroelasticity for the periodic response to very fast loading (‘$\mathrm {vf}$’) is

(3.12)\begin{equation} u_{s,{vf}}(x,t) = \frac{A}{2}\left[ 1 -x -\exp\left({-}x\sqrt\frac{\rm \pi}{T}\right)\cos\left(\frac{2{\rm \pi}{}t}{T}-x\sqrt\frac{\rm \pi}{T}\right)\right] \end{equation}

and

(3.13)\begin{align} \phi_{f,{vf}}(x,t) &= \phi_{f,0} - (1-\phi_{f,0})\frac{A}{2}\left\{ 1-\sqrt{\frac{\rm \pi}{T}}\exp\left({-}x\sqrt\frac{\rm \pi}{T}\right)\right.\nonumber\\ &\left.\quad \times\left[\cos\left(\frac{2{\rm \pi}{}t}{T}-x\sqrt\frac{\rm \pi}{T}\right)-\sin\left(\frac{2{\rm \pi}{}t}{T}-x\sqrt\frac{\rm \pi}{T}\right)\right]\right\}, \end{align}

where the first two terms in $u_{s,{vf}}$ and in $\phi _{f,{vf}}$ give the static far-field compression, which is also the (linearised) overall mean compression. This solution confirms that the oscillations will be increasingly localised near the piston as $T$ decreases, featuring near-piston oscillations with an amplitude proportional to $1/\sqrt {T}$ that decay exponentially in space over a characteristic distance $\sqrt {T}$. This solution is illustrated and discussed further in Appendix B.

3.5. Response to very slow loading (quasi-static solution)

In the full linear-poroelastic solution given previously, the porosity and displacement fields ((3.9) and (3.10)) converge to the quasi-static limits ${\phi _{f,{lpe}}}(x,t) \to \phi _{f,0}-(1-\phi _{f,0})a(t)$ and $u_{s,{lpe}}\to {}a(t)(1-x)$ as $T\to \infty$. This limit is a uniform state of strain in which poroelastic transients are fast relative to the loading period, and are therefore negligible. The fully nonlinear problem can be solved analytically in the same limit by taking $\partial {\phi _f}/\partial {t}\to {}0$. The resulting quasi-static (‘$\mathrm {qs}$’) solution is

(3.14)\begin{gather} u_{s,{qs}}(x,t)= \frac{a(t)}{1-a(t)} (1-x), \end{gather}

(3.15)\begin{gather} {\phi_{f,{qs}}}(t)= \frac{\phi_{f,0}-a(t)}{1-a(t)}, \end{gather}

and

(3.16)\begin{equation} v_{f,{qs}}(x,t)={-}\left[\frac{1-\phi_{f,0}}{\phi_{f,0}-a(t)}\right]\left\{\frac{1-x}{[1-a(t)]^2}\right\}\dot{a}(t). \end{equation}

This solution is kinematically exact for arbitrarily large values of $A$, but is only valid for $T\gg {}1$. Note that this expression for $\phi _{f,{qs}}$ is the same as that in (3.1) for $\langle {\phi _f}\rangle (t)$ because the quasi-static porosity is spatially uniform and must therefore be equal to the average porosity.

3.6. Scaling quantities

In the absence of a net flow, fluid motion is directly related to changes in porosity. Hence, we present our solutions and results below in terms of the change in porosity with respect to the overall average porosity,

(3.17)\begin{equation} \Delta{\phi_f} \equiv \phi_f - \overline{\langle{\phi_f}\rangle}. \end{equation}

This definition accounts for the non-zero mean compression.

The various results given previously suggest simple scaling relationships for the magnitudes of $\Delta {\phi _f}$ and $v_f$ in terms of the dimensionless control parameters $A$ and $T$. In particular, the magnitude of $\Delta {\phi _f}$ is captured by the spatially averaged change in porosity at mid-cycle,

(3.18)\begin{equation} \Delta\phi_f^M= \left|\langle{\Delta \phi_f}\rangle(T/2)\right| =\frac{\phi_{f,0}-A}{1-A}-\overline{\langle\phi_f\rangle}. \end{equation}

The quasi-static solution suggests that appropriate scales for the magnitude of the solid and fluid velocity are

(3.19)\begin{equation} v_s^{{\ast}}=\frac{2A}{T}, \end{equation}

and

(3.20)\begin{equation} v_f^{{\ast}}= \left( \frac{1-\phi_{f,0}}{\phi_{f,0}} \right) \frac{2A}{T}. \end{equation}

An appropriate scale for the fluid flux is therefore

(3.21)\begin{equation} q_f^{{\ast}}= \phi_{f,0} v_f^{{\ast}}. \end{equation}

4. Numerical solution

We solve the full nonlinear problem (§ 2.9) numerically in MATLAB using a Chebyshev spectral method in space and an implicit Runge–Kutta method in time, as described in more detail in Appendix D. We provide an example code in Fiori, Pramanik & MacMinn (2023). In figure 2, we illustrate the basic phenomenology of the response of a high-porosity material ($\phi _{f,0}=0.75$) to periodic loading at large amplitude ($A=0.2$) and moderate period ($T=0.3{\rm \pi}$) for one cycle in the periodic regime.

Figure 2. Response of a high-porosity material ($\phi _{f,0}=0.75$) to periodic loading at high amplitude ($A=0.2$) and moderate period ($T=0.3{\rm \pi}$). In (a–d), we show the evolution of $u_s$, $\Delta \phi _f$, $v_s$ and $v_f$ (all normalised) at times $t=nT$ to $(n+1)T$ in increments of $0.1T$ (dark to light) during one cycle in the periodic regime, where $n$ is an integer. In (a–d), we plot all fields against the Lagrangian spatial coordinate $X=x-u_s(x,t)$ for visual clarity. In (e), we plot $\sigma ^{\prime }$ against $t$ at fixed values of $X$ from $0$ to $1$ (light to dark). In ( f), we plot a phase portrait of $\sigma ^{\prime }(a,t)$ against $a(t)$ for $t=0$ to $t=20T$, emphasising the last cycle (dotted black). In all plots but ( f), we distinguish between the loading half of the cycle ($\dot {a}(t)>0$; solid curves) and the unloading half of the cycle ($\dot {a}(t)<0$; dashed curves).

All quantities considered, displacement, change in porosity, solid velocity, fluid velocity and effective stress, are largest in magnitude at the left boundary, where the material is forced, and smallest in magnitude at the right boundary, where the material is stationary, with a magnitude envelope that decreases monotonically from left to right. The displacement and both velocities vanish at the right boundary, as required. Note that these features depend greatly on the boundary conditions for both the solid and the fluid. The flow and deformation will focus toward boundaries where inflow and outflow are permitted, which here is the left-hand side. Reversing the permeability of the two boundaries (i.e. an impermeable moving boundary and a permeable fixed boundary) would instead focus the flow and deformation toward the right-hand side, roughly reversing the spatial profile of $\Delta {\phi _f}$ and producing a much more uniform profile of $v_f$. However, the latter scenario is identical to the present one when viewed from a moving frame that follows the left boundary, $x^\prime =1+a(t)-x$.

Figure 2(e, f) shows the effective stress normalised by the nominal strain at maximum compression, $\sigma ^{\prime }/A$. The left column shows that the response is out of phase with the loading, with a phase shift that varies with the Lagrangian spatial coordinate $X=x-u_s(x,t)$. For example, the stress at the left boundary leads the motion of the left boundary by about $0.15T$ (i.e. the moment of maximum $|\sigma ^\prime (a,t)|$ occurs about $0.15T$ before the moment of maximum $a(t)$), whereas the stress at the right boundary lags the motion of the left boundary by a similar amount. In addition, the material near the left boundary experiences strong compression during most of the loading phase ($\sigma ^{\prime }<0$ for $\dot {a}>0$) and mild tension during much of unloading ($\sigma ^{\prime }>0$ for $\dot {a}<0$). This hysteresis is highlighted by the large area enclosed by the loop in the phase portrait (right column), and it originates in the stronger relative role of viscous drag during moderate to fast loading.

Most of the features illustrated in figure 2 are qualitatively consistent with linear poroelasticity, although the quantitative accuracy of linear poroelasticity depends on the deformation parameters as discussed in the next section. A key qualitative feature introduced by nonlinearity is that the response during loading is not necessarily symmetric with the response during unloading. For example, the minimum values of $\Delta {\phi _f}(a,t)$ and $v_f(a,t)$ are much larger in magnitude than their maximum values and the stress loop is not symmetric about any axis. We explore this asymmetry in more detail in § 5.

4.1. Comparison with analytical solutions

We next compare the numerical solution to the linear-poroelastic and nonlinear quasi-static analytical solutions described in § 3, each of which is appropriate for a specific range of $A$ and $T$. The aim of this comparison is to quantify these ranges of validity and to examine the convergence of the numerical results to each of these special cases. To do so, we calculate all three solutions over a wide range of $A$ and $T$ and then calculate the root-mean-square (r.m.s.) relative difference between the numerical and linear-poroelastic solutions (figure 3), and between the numerical and nonlinear quasi-static solutions (figure 4). We calculate these differences using $\phi _f(a,t)$ during one cycle in the periodic regime. In both figures, we plot these differences against $T$ for several values of $A$ (left panels) and then against $A$ for several values of $T$ (right panels).

Figure 3. The r.m.s. relative difference between the full numerical solution and the linear-poroelastic analytical solution (from (3.10)) based on $\phi _f(a,t)$ during one cycle in the periodic regime. On (a), we plot the difference against $T$ for fixed values of $A$ ranging from $0.02$ to $0.2$ (dark to light). On (b), we plot the same difference against $A$ for fixed values of $T$ ranging from $0.001{\rm \pi}$ to $10{\rm \pi}$ (dark to light).

Figure 4. The r.m.s. relative difference between the full numerical solution and the nonlinear quasi-static solution (from (3.15)) based on $\phi _f(a,t)$ during one cycle in the periodic regime. On (a), we plot the difference against $T$ for fixed values of $A$ ranging from $0.02$ to $0.2$ (dark to light). On (b), we plot the same difference against $A$ for fixed values of $T$ ranging from $0.001{\rm \pi}$ to $10{\rm \pi}$ (dark to light).

As expected, figure 3 shows good agreement between the numerical and linear-poroelastic solutions for small deformations, worsening as $A$ increases. The difference scales with $A^2$ for a given value of $T$ (right panel), as expected from linear poroelasticity, which is first order in strain. The difference is insensitive to $T$ for $T\gtrsim {}1$, but scales as $T^{-1}$ for faster periods. Decreasing $T$ leads to increasingly strong localisation at the left boundary, and hence increasingly large deformations, even for small $A$, as expected from the early time and very fast analyses above (§§ 3.2–3.4). We explore this localisation in more detail in § 5.

As expected, figure (4) shows good agreement between the numerical and nonlinear quasi-static solutions for large periods, worsening as $T$ decreases. The difference scales with $T^{-1}$ for $T\gtrsim {}1$ and with $T^{-1/2}$ for shorter periods (left panel); the difference also scales with $A$ (right panel), consistent with the scaling of the non-quasi-static terms in (3.10). Based on these results, we distinguish between ‘slow loading’ (SL; $T\gtrsim {}{\rm \pi}$), where spatial variations in porosity are relatively small, and ‘fast loading’ (FL; $T\lesssim {}0.1{\rm \pi}$), where spatial variations in porosity are relatively large. For very slow loading ($T\gg {}1$), the porosity is uniform and the response is quasi-static (see § 3.5). For very fast loading ($T\ll {}1$), the oscillations are localised near the left boundary and the right portion of the material is in a state of static compression (see § 3.4).

5. Parameter study

We next examine and compare the poroelastic response for SL and FL as a function of $T$, $A$ and $\phi _{f,0}$. We focus on the evolution of $u_s$, $\phi _f$ and $q_f=\phi _fv_f$ in space and in time, and on the phase behaviour of $\sigma ^\prime (a,t)$. We conclude by considering the parameter ranges that would be relevant to various biological examples.

5.1. Impact of loading period

To visualise the distinct poromechanical responses for FL and SL, and the transition between them, we plot $u_s$, $\Delta {\phi _f}$ and $q_f$ (all normalised) over one cycle in the periodic regime for $\phi _{f,0}=0.75$, $A=0.1$, and four different values of $T$: two for FL (left two columns) and two for SL (right two columns) (figure 5).

Figure 5. Evolution of normalised $u_s$ (a–d) and $\Delta {\phi _f}$ (e–h) at times $t=nT$ to $(n+1)T$ in increments of $0.1T$ (dark to light), where $n$ is an integer, during one cycle in the periodic regime for $A=0.1$, $\phi _{f,0}=0.75$, and $T=0.03{\rm \pi}$ (a,e,i,m), $0.1{\rm \pi}$ (b, f, j,n), ${\rm \pi}$ (c,g,k,o) and $10{\rm \pi}$ (d,h,l,p). As in figure 2, we plot all fields in (a–h) against the Lagrangian spatial coordinate $X=x-u_s(x,t)$ for clarity. In (i–l), we plot normalised $q_f$ against $t$ for 10 different values of $X$ from $0$ to $1$ (dark to light). We distinguish between the loading half of the cycle ($\dot {a}>0$; solid curves) and the unloading half of the cycle ($\dot {a}<0$; dashed curves). In (m–p), we plot phase portraits of $\sigma ^\prime (a,t)$ against $a(t)$ from $t=0$ to $t=20T$, emphasising the last cycle (dotted black).

Figure 5(a–d) show that, for FL (i.e. $T=0.03{\rm \pi}$ and $0.1{\rm \pi}$), the displacement is highly nonlinear in $X$ (and in $x$), with substantial differences between the loading and unloading phases. As $T$ increases, transitioning into SL (i.e. $T=1{\rm \pi}$ and $10{\rm \pi}$), the displacement is increasingly linear in $X$ and converges toward the quasi-static solution, which is fully determined by the instantaneous value of $a$ and is thus symmetric between loading and unloading. For all values of $T$, the displacement is of characteristic size $A$ at the left and vanishes at the right.

For SL, $\Delta {\phi _f}$ is uniform in space and varies only in time, per the quasi-static solution. The material is in a uniform state of compression, fully determined by the instantaneous value of $a$ and thus symmetric between loading and unloading and fully relaxed at the beginning/end of each cycle. As $T$ decreases, transitioning into FL, $\Delta {\phi _f}$ becomes increasingly localised near the left boundary and also increasingly asymmetric between loading and unloading (figure 5e). For FL, the material experiences a substantial amount of tension in the left portion of the domain during unloading, despite the overall mean compression. Tension emerges for FL because this regime is, by definition, one where the rate of loading is much faster than the poroelastic relaxation time, so the left boundary must pull the material to the left during unloading. The right portion of the domain experiences an overall more limited range of porosities and remains compressed throughout the cycle, never reaching a state of tension or even full unloading. The latter feature is also visible in the corresponding displacement fields.

The fluid flux $q_f=\phi _fv_f$ is particularly relevant to the transport of solutes since it drives advection. Fluid leaves the domain during the loading phase of the cycle ($q_f(x=a,t)=q_f(X=0,t)<0$ when $\dot {a}>0$) and enters the domain during unloading. For SL, $q_f$ is entirely in phase with the loading, but opposite in sign. For FL, the peak value of $q_f$ in the interior exhibits a lag relative to the peak value of $\dot {a}$, and this lag increases with $x$. Note that the fluid flux is orders of magnitude larger for FL than for SL because the rate of loading is orders of magnitude faster, but this variation is largely scaled out by normalisation with $q_f^\ast (\phi _{f,0},A,T)\propto {}T^{-1}$, per (3.21) (figure 5i–l).

The extreme values of $q_f$ at the left boundary are larger in loading than in unloading, but progressively more symmetric as $T$ increases. This asymmetry is a result of the kinematic and constitutive nonlinearity of large deformations. We show in Appendix E that this asymmetry originates in the nonlinear kinematics of large deformations, meaning that it emerges from the nonlinear model during large deformations even with linear elasticity and constant permeability. This asymmetry is then strongly amplified by deformation-dependent permeability (relative to constant permeability) and slightly suppressed by Hencky elasticity (relative to linear elasticity). The latter occurs because Hencky elasticity stiffens in compression, resulting in a larger poroelastic diffusivity and therefore weaker localisation during loading (see also Appendix A).

The asymmetry between loading and unloading is also highlighted by the evolution of $\sigma ^{\prime }(a,t)$ (figure 5m–p). For FL, $\sigma ^\prime (a,t)$ exhibits hysteresis: for the same value of $a(t)$, the value of $\sigma ^\prime (a,t)$ is considerably higher in magnitude during loading than during unloading (and tensile during much of the unloading phase). This hysteresis decreases as $T$ increases, such that, for SL, $\sigma ^\prime (a,t)$ is modest in magnitude, fully compressive and symmetric between loading and unloading (i.e. non-hysteretic); as a result of the latter, the phase portrait for SL is a single curve (rather than a loop). These phase portraits also illustrate the convergence to the periodic regime: in all cases, the overall response grows less extreme as the transient component decays exponentially over the first $\sim {}T^{-1}$ cycles.

Since a key difference between SL and FL is that the deformation is uniformly distributed in the former and localised toward the left in the latter, we study the emergence of this localisation to further quantify the transition from SL to FL. In figure 6, we plot the normalised maximum and minimum values of $\Delta {\phi _f}$, $q_f$ and $\sigma ^\prime$ at the left boundary ($X=0$) and then either at the right boundary ($X=1$) for $\Delta {\phi _f}$ and $\sigma ^\prime$ or at the material mid-point ($X=1/2$) for $q_f$. The latter is necessary because $q_f$ vanishes at $X=1$. We consider a range of periods $T$ at two fixed values of $A$ for $\phi _{f,0}=0.75$.

Figure 6. We illustrate the transition from FL to SL as $T$ increases by plotting the normalised maximum and minimum values $\Delta {\phi _f}$ (a), $q_f$ (b) and $\sigma ^\prime$ (c) against $T$ for $A=0.1$ (solid lines) and $A=0.001$ (dotted lines) during the periodic regime. We show the maximum values (dark colours) and minimum values (light colours) of all three quantities at the left boundary ($X=0$, blue curves), of $\Delta {\phi _f}$ and $\sigma ^\prime$ at the right boundary ($X=1$, red curves), and of $q_f$ at the material midpoint ($X=1/2$, red curves).

Figure 6 shows that, for SL, the deformation is uniform, varying only in time. The maxima and minima of both $\Delta {\phi _f}$ and $\sigma ^{\prime }$ remain separate, but their left and right values converge. The flux $q_f$ remains linear in space, even for SL (see § 3.5). As $T$ decreases, the deformation is progressively localised at the left, such that the right is increasingly static. At the right, the maximum and minimum values of $\Delta {\phi _f}$ and $q_f$ converge to zero while the maximum and minimum values of $\sigma ^\prime$ converge to weak compression, as also illustrated in figure 5.

We find that a small-amplitude deformation ($A=0.001$) and a large-amplitude deformation ($A=0.1$) exhibit qualitatively similar features. However, large deformations lead to an increasingly strong asymmetry between the maxima and minima of all quantities at the left boundary as $T$ decreases. In particular, the curves are biased downward, such with higher magnitudes reached in loading (minima) than in unloading (maxima). This asymmetry is does not occur for small deformations, for which the maxima and minima of all quantities are symmetric in magnitude. Note also that the maxima of $\Delta {\phi _f}$ and $\sigma ^{\prime }$ increase as $T$ decreases for both values of $A$, whereas the maximum of $q_f$ is independent of $T$ for small deformations but decreases with $T$ for large deformations, consistent with figure 5.

5.2. Impact of loading amplitude

We next examine the impact of loading amplitude $A$ on the poromechanical response. We first assess the interaction between amplitude and period by plotting $\Delta {\phi _f}$ against $X$ at mid-cycle for two different values of $T$ (one for SL and one for FL) and for five different values of $A$, ranging from small to large deformations ($A=0.002$ to $0.2$; figure 7).

Figure 7. Profiles of $\Delta {\phi _f}$ at mid-cycle during the periodic regime for $\phi _{f,0}=0.75$ and for five values of $A$ ranging from small to large deformations (dark to light), each for two different periods: $T=0.1{\rm \pi}$ (FL; solid) and $T=10{\rm \pi}$ (SL; dashed). Panels (a) and (b) are non-normalised and normalised, respectively, showing that this normalisation scales out the leading-order impact of $A$ on $\Delta {\phi _f}$.

Figure 7(a) shows that, for SL, the non-normalised value of $\Delta {\phi _f}$ is uniform in $X$ and increases with $A$, as expected. For FL, the non-normalised value of $\Delta {\phi _f}$ is increasingly non-uniform with $X$ as $A$ increases. Relative to the SL case, the deformation is increasingly amplified at the left boundary and suppressed at the right boundary (recall that the mean value of $\Delta {\phi _f}$ is independent of $T$). The right panel shows the same results, but now normalised by $\Delta \phi _f^M$ (3.18). By definition, all of the SL curves collapse onto $\Delta {\phi _f}/\Delta \phi _f^M=-1$ ($\Delta \phi _f^M$ is the magnitude of the quasi-static change in porosity at mid-cycle). The FL curves also nearly collapse onto a master curve, suggesting that this normalisation scales out the leading-order impact of $A$, even for large deformations for FL. The largest deviations from this collapse are near the left boundary, where nonlinearities are particularly pronounced.

In figure 8, we focus on FL by fixing $T=0.1{\rm \pi}$ and plotting the evolution of $\Delta {\phi _f}$, $q_f$ and $\sigma ^\prime (a,t)$ for three different values of $A$, ranging from small to large deformations.

Figure 8. Evolution of normalised $\Delta {\phi _f}$, $q_f$ and $\sigma ^\prime (a,t)$ for $\phi _{f,0}=0.75$, $T=0.1{\rm \pi}$, and $A=0.01$ (a,d,g), $0.1$ (b,e,h) and $0.2$ (c, f,i). In (a–c), we plot $\Delta {\phi _f}$ against the Lagrangian spatial coordinate $X=x-u_{s}(x,t)$ at times $t=nT$ to $(n+1)T$ in increments of $0.1T$ (dark to light) during one cycle in the periodic regime. In (d–f), we plot $q_f$ against $t$ for ten different values of $X$ from $0$ to $1$ (dark to light); we distinguish between loading ($\dot {a}>0$; solid curves) and unloading ($\dot {a}<0$; dashed curves). In (g–i), we plot phase portraits of $\sigma ^\prime (a,t)$ against $a(t)$ from $t=0$ to $t=20T$, emphasising the last cycle (dotted black lines).

We find that increasing $A$ amplifies the asymmetry in the extreme values of $\Delta {\phi _f}$ and $q_f$ at the left boundary in loading and unloading, as noted in figure 6. The normalisation scales out the leading-order impact of $A$ on all of these quantities, as noted in figure 7; the non-normalised values of all three quantities would vary by two orders of magnitude across this range of $A$. As $A$ increases, $\sigma ^\prime (a,t)$ exhibits more hysteresis and larger normalised magnitudes, in accordance with the amplified asymmetry.

We further explore the impact of $A$ in figure 9 by plotting the normalised maxima and minima of $\Delta {\phi _f}$, $q_f$ and $\sigma ^\prime$ on the left and on the right, as in figure 6, but now against $A$ for two different values of $T$, showing the smooth transition from small deformations ($A \lesssim 0.01$) to large deformations ($A \gtrsim 0.01$). For small deformations, the normalised maxima and minima of all quantities at the left and at the right become independent of $A$, and the extreme values of $\Delta {\phi _f}$ and $q_f$ become symmetric. For SL, $\Delta {\phi _f}$ and $\sigma ^\prime$ become uniform in space, such that their normalised left and right values are equal and depend only weakly on $A$ for the largest deformations shown here (e.g. $A\gtrsim {}0.1$). For FL, the values of all three quantities at the left become increasingly asymmetric and biased downward as $A$ increases.

Figure 9. We illustrate the transition from small to large deformations as $A$ increases by plotting the normalised maximum and minimum values $\Delta {\phi _f}$ (a), $q_f$ (b) and $\sigma ^\prime$ (c) against $A$ for $T=0.1{\rm \pi}$ (solid lines) and $T=10{\rm \pi}$ (dotted lines) during the periodic regime. We show the maximum values (dark colours) and minimum values (light colours) of all three quantities at the left boundary ($X=0$, blue curves), of $\Delta {\phi _f}$ and $\sigma ^\prime$ at the right boundary ($X=1$, red curves) and of $q_f$ at the material midpoint ($X=1/2$, red curves).

5.3. Impact of initial porosity

Finally, we consider the initial porosity $\phi _{f,0}$. In figure 10, we plot the distribution of $\Delta {\phi _f}$ at mid-cycle for three different values of $\phi _{f,0}$ for fixed $A=0.01$ and for two different values of $T$, one for FL ($T=0.1{\rm \pi}$) and SL ($T=10{\rm \pi}$). For both SL and FL, the non-normalised values of $\Delta {\phi _f}$ (left panel) increase in magnitude as $\phi _{f,0}$ decreases because the same change in total volume is a larger portion of the initial fluid volume. The right panel shows that normalisation by $\Delta {\phi _f^M}$ scales out the leading-order impact of varying $\phi _{f,0}$ on $\Delta {\phi _f}$, again with small deviations at the left boundary for FL, when nonlinearities are most pronounced.

Figure 10. Impact of $\phi _{f,0}$ on $\Delta {\phi _f}$ at mid-cycle in the periodic regime for $A=0.1$ and for two different periods: $T=0.1{\rm \pi}$ (FL; solid) and $10{\rm \pi}$ (SL; dashed). We show results for three values of $\phi _{f,0}$ (increasing from dark to light). Panels (a) and (b) are non-normalised and normalised, respectively, showing that this normalisation scales out the leading-order impact of $\phi_{f,0}$ on $\delta{\phi_f}$.

We next plot the evolution of $\Delta {\phi _f}$ and $\sigma ^\prime$ in time for $A=0.01$ and $T=0.1{\rm \pi}$ for the same three values of $\phi _{f,0}$ (figure 11), confirming that normalisation scales out the primary impact of $\phi _{f,0}$ on $\Delta {\phi _f}$ and on $\sigma ^\prime (a,t)$.

Figure 11. Impact of $\phi _{f,0}$ on the evolution of $\Delta {\phi _f}$ and $\sigma ^\prime (a,t)$ for $A=0.01$, $T=0.1{\rm \pi}$, and $\phi _{f,0}=0.25$ (a,d), $0.5$ (b,e) and $0.75$ (c, f). In (a–c), we plot $\Delta {\phi _f}$ against the Lagrangian spatial coordinate $X=x-u_s$ at times $t=nT$ to $(n+1)T$ in increments of $0.1T$ (dark to light) during one cycle in the periodic regime, and we distinguish between loading ($\dot {a}>0$; solid curves) and unloading ($\dot {a}<0$; dashed curves). In (d–f), we plot phase portraits of $\sigma ^\prime (a,t)$ against $a(t)$ from $t=0$ to $t=20T$, emphasising the last cycle (dotted black lines).

5.4. FL in biological examples

We conclude by considering appropriate parameter ranges for several examples of periodic loading in soft biological tissues (table 1). Based on these values, we calculate the poroelastic timescale $T_{{pe}}$ and the relative dimensionless loading period $\tilde {T}=T/T_{{pe}}$ for each example to understand the ranges of relevance for $\tilde A$, $\tilde T$ and $\phi _{f,0}$.

Table 1. Material and loading parameters for some examples of biological materials.

Table 1 indicates that, for a variety of soft tissues, deformations are in the range of the dimensionless amplitudes $\tilde {A}$ considered here, with nearly all being near the upper end. The range of $\phi _{f,0}$ is slightly wider than the range considered here, but we showed in § 5.3 that this value has little impact on the normalised mechanical response of the material. The range of poroelastic timescales $T_{{pe}}$ and respective loading frequencies suggest that dimensionless loading periods $\tilde {T}$ span a wide range, with many corresponding to FL. These values justify our analysis and underscore the importance of characterising the nonlinearity of large poromechanical deformations during FL in particular.

6. Conclusions

We have provided an analysis of the poromechanical coupling between large deformations and fluid flow in a periodically loaded soft porous material. To do so, we used a kinematically rigorous one-dimensional continuum model with Hencky elasticity and Kozeny–Carman permeability. In particular, we examined the roles of the three dimensionless control parameters: the initial porosity $\phi _{f,0}$, the loading amplitude $A$ and the loading period $T$.

We began by deriving several analytical solutions from linear poroelasticity, an early time solution, a full solution and a solution for the response to very fast loading (Stokes's second problem), as well as a quasi-static solution to the fully nonlinear problem. The former are valid for small deformations, which corresponds to $A\ll {}1$ and also, less obviously, to $A\ll {}\sqrt {T}$. The quasi-static solution is valid for very slow loading ($T\gg {}1$) but arbitrarily large amplitudes. We then compared these solutions with our numerical results, highlighting the existence of two mechanical regimes: slow loading (SL), where the loading is much slower than the poroelastic relaxation time $T_{pe}$, and fast loading (FL), where the loading is much faster than $T_{pe}$.

We then showed that the material response to SL ($T\gtrsim {\rm \pi}$) is an increasingly uniform deformation throughout the domain, approaching the quasi-static solution for very slow loading. For FL ($T\lesssim 0.1{\rm \pi}$), the deformation is non-uniform and increasingly localised near the left boundary. In the limit of very fast loading, this localisation is such that the left portion of the material oscillates while the right portion is in a state of static compression. We showed that FL is also characterised by asymmetry between loading and unloading, with a larger change in porosity and higher fluid flux magnitudes during loading (when fluid is squeezed out) than during unloading (when fluid is sucked back in). This asymmetry originates in the kinematic nonlinearity of large deformations and is amplified by the localisation of the deformation near the left boundary as $T$ decreases (itself a linear effect) and by the nonlinearity of deformation-dependent permeability as $A$ increases. Thus, this asymmetry emerges from the fast and large deformation of an elastic and initially homogeneous material, and is therefore purely poromechanical; it does not occur during slow loading. Asymmetry between quasi-static loading and unloading can instead result from wall friction in confined geometries and/or constitutive hysteresis due to microstructural changes, as has been observed experimentally for some soft porous materials (e.g. Sobac, Colombani & Forterre 2011; Hewitt et al. 2016; Lutz, Wilen & Wettlaufer 2021).

Through the analysis on the fluid flow at different fixed material points in the domain, we also showed that faster loading leads to an increasing delay in the interior of the domain relative to the motion of the left boundary. Although the motion of the fluid itself is reversible, these features are likely to have an important impact on irreversible phenomena. For example, these results have interesting implications for the deformation-driven transport and mixing of solutes, which we consider in detail in a companion study.

The evolution of the stress at the left boundary as a function of the piston position revealed that, for faster loading and larger deformations, the force required to compress the material is much larger than the force required to pull it back. This asymmetry occurs because the viscous drag associated with outflow compresses the material against the outlet, reducing the permeability and resisting further fluid expulsion, whereas the viscous drag associated with inflow stretches the material away from the outlet, increasing the permeability and easing further fluid influx.

Finally, we showed that a larger initial porosity $\phi _{f,0}$ leads to much lower fluid fluxes that can be accounted for by normalisation.

Our results elucidate the local and global poromechanical behaviour of soft porous media during periodic loading over a wide range of $\phi _{f,0}$, $A$ and $T$. Having used relatively generic constitutive models, we expect our qualitative insights to be robust across a wide range of materials; however, it is straightforward to adapt our approach to other constitutive models (see Appendix A).