Flow and Deformation due to Periodic Loading in a Soft Porous Material

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,

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$.

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.

We work in an Eulerian reference frame, in which the solid displacement field is $\mathbf{u_{s}}=\mathbf{x}-\mathbf{X}(\mathbf{x},t)$, with $\mathbf{X}(\mathbf{x},t)$ the reference position of the material point that at time $t$ occupies the position $\mathbf{x}$. We choose $\mathbf{X}(\mathbf{x},0)\equiv\mathbf{x}$ so that $\mathbf{u_{s}}(\mathbf{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

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 $\mathbf{\hat{e}_{x}}$ is the unit vector in the $x$-direction.

In 1D, 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

Continuity for this 1D system can be written

which together imply that the total flux $q=\phi_{f}v_{f}+(1-\phi_{f})v_{s}$ is uniform in space, $\partial{q}/\partial{x}=0$.

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

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):

where $k_{0}\equiv{}k(\phi_{f,0})$ is the reference permeability. Kozeny-Carman permeability is a common choice for gels and soft tissues Sacco et al. (2014); Malandrino et al. (2014); Rahbari et al. (2017); Gao and 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 behavior.

Combining equations (4) and (5), we arrive at the nonlinear flow equations:

where the total flux $q$ is again

The fluid velocity and the solid velocity are then given by

and the local fluid flux is

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^{\prime}}$,

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

In 1D, equation (12) implies that

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

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 1D, as considered here, any elasticity law can be written in the form $\sigma^{\prime}=\sigma^{\prime}(\phi_{f})$. Thus, equation (7) can be rewritten as a nonlinear advection-diffusion equation:

where the nonlinear composite constitutive function

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 and Chen (2002) and sometimes for soft biological tissues Marchesseau et al. (2010); Fraldi et al. (2018). 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 behavior, 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. Ehlers et al. (2009); Sengers et al. (2004)); 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 and MacMinn (2018)

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

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

For comparison with the fully nonlinear model, we linearise the relations above 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, equation (7) reduces to the linear-poroelastic diffusion equation,

where Hencky elasticity reduces to linear elasticity,

and $D_{f,0}=k_{0}\mathcal{M}/\mu$ is the constant linear-poroelastic diffusivity. With appropriate initial and boundary conditions, Equation (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 equations (17) and (19), respectively. The linearised boundary conditions for the left boundary are

where $a(t)$ is again given by equation (1). Note that these conditions are linearized relative to equations (18) by virtue of being applied at $x=0$ rather than at $x=a(t)$.

We make the above problem dimensionless via the scaling

where $T_{\mathrm{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 §II.7.3), the nonlinear flow equation can be rewritten in dimensionless form as

with nonlinear-poroelastic diffusivity

elasticity law

initial conditions

left boundary conditions

and right boundary conditions

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

The above 1D 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}$.

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

The average of $\langle\phi_{f}\rangle(t)$ in time over any integer number of loading cycles is then

for any $n\geq 0$ and integer $m\geq 1$. Note that both the spatial and overall averages are negative because the loading has a nonzero mean (i.e., $\bar{a}=A/2$), so the material is on average compressed.

The linear problem posed in section II.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,

The corresponding porosity field can be derived from equation (33) via equation (3), and is given by

This solution is valid until the deformation spans the domain. Introducing the penetration length $\delta_{\mathrm{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_{\mathrm{et}}(t)\sim{}2\sqrt{t}$. We confirm this reasoning in Appendix B. Thus, the above solution is valid for $\delta_{\mathrm{et}}\lesssim{}1\,\to\,t\lesssim{}1/4$.

Having originated from linear poroelasticity, the above solution is limited to small deformations ($A\ll 1$). Naively, 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_{\mathrm{et}}(t)$ at early times. As a result, we expect the maximum strain near the left boundary to be roughly of size $\mathrm{max}[a(t)/\delta_{\mathrm{et}}(t)]\sim{}A/\sqrt{T}$. More precisely, equation (34) can be used to show that the extreme values of ${\phi_{f,\mathrm{et}}}$ will always occur at the left boundary and that the evolution of ${\phi_{f,\mathrm{et}}}$ at the left boundary is given by

where

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

and

where $\mathcal{I}_{\mathrm{max}}=\mathcal{I}(698/1909)\approx{}0.9491$ and $\mathcal{I}_{\mathrm{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 $\mathrm{min}({\phi_{f,\mathrm{et}}})=0$ and $\mathrm{max}({\phi_{f,\mathrm{et}}})=1$, respectively. In practise, 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 eq. 23). 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}$.

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

The corresponding porosity field can be derived from equation (39) via equation (3) and is given by

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

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,\mathrm{lpe}}}$ can be divided into three parts:

1. 1.

   a uniform quasi-static part proportional to $a(t)$, which is the linearised form of the nonlinear quasi-static solution derived below;

2. 2.

   an early-time transient that decays exponentially in time at a rate that is independent of $A$ and $T$; and

3. 3.

   a periodic forced response with period $T$.

The early-time transient captures the early-time solution derived above, 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.

As noted above, 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

and

where the first two terms in $u_{s,\mathrm{vf}}$ and in $\phi_{f,\mathrm{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.

In the full linear-poroelastic solution above, the porosity and displacement fields (equations 39 and 40) converge to the quasi-static limits ${\phi_{f,\mathrm{lpe}}}(x,t)\to\phi_{f,0}-(1-\phi_{f,0})a(t)$ and $u_{s,\mathrm{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

and

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

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,

This mean change in porosity accounts for the non-zero mean compression.

The various results above 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,

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

and

An appropriate scale for the fluid flux is therefore

We next compare the numerical solution to the linear-poroelastic and nonlinear quasi-static analytical solutions described in §III, 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 (RMS) 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).

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 (§III.2–III.4). We explore this localisation in more detail in §V.

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 equation (40). Based on these results, we distinguish between “slow loading” (SL; $T\lesssim{}0.1\pi$), where spatial variations in porosity are relatively small, and “fast loading” (FL; $T\gtrsim{}\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 §III.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 §III.4).

To visualise the distinct poromechanical responses for SL and FL, 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 SL (left two columns) and two for FL (right two columns) (figure 5).

The first row of figure 5 shows that, for FL (i.e., $T=0.03\pi$ and $0.1\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\pi$ and $10\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 5, second row, left column). 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_{f}v_{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 equation (51) (figure 5, third row).

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 5, fourth row). 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 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 §III.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 minimia 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.

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$; fig. 7).

The left panel in figure 7 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}$ (eq. 48). 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 captures 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\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.

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 captures 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}$ at the left and at 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.

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\pi$) and SL ($T=10\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}}$ against captures 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.

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

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_{\mathrm{pe}}$ and the relative dimensionless loading period $\tilde{T}=T/T_{\mathrm{pe}}$ for each example to understand the ranges of relevance for $\tilde{A}$, $\tilde{T}$ and $\phi_{f,0}$.

Table 1 shows 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 §V.3 that this value has little impact on the normalised mechanical response of the material. The range of poroelastic timescales $T_{\mathrm{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.