Cyclic loading of a heterogeneous non-linear poroelastic material

Abstract

Cyclic loading is a common feature in poroelastic systems, the material response depending non-trivially on the exact form of boundary conditions, pore structure, and mechanical properties. The situation becomes more complex when heterogeneity is introduced in the properties of the poroelastic material, yet heterogeneity too is common in physical poroelastic structures. In this paper, we analyse the behaviour of a soft porous material in response to a uniaxial cyclic stress or displacement, with a focus on understanding how this response is affected by continuous heterogeneity in the stiffness or permeability. Our work is motivated by observed altered material properties of the diseased tendon, but the framework we develop and analyse is generically applicable. We construct a one-dimensional non-linear poroelastic model, assuming Darcy flow through the pores of the solid skeleton which we assume has neo-Hookean elasticity. The system is driven by an applied uniaxial cyclic stress or a uniaxial cyclic displacement at one boundary. Heterogeneity in the stiffness or permeability profile is imposed via a Gaussian bump function. By exploring a range of loading frequencies together with magnitudes and locations of heterogeneity, we characterise the effect of heterogeneity on the response of the material, and show that the response of the system to an applied stress is qualitatively distinct from the response to an applied displacement. For example, damage corresponding to a local decrease in stiffness locally increases strain under an applied stress and locally perturbs the flux under an applied displacement. In addition, the closer the damage is to the loaded boundary, the more sensitive the response is to the damage. Our analysis of this simple model provides a foundation for understanding how heterogeneity affects the poroelastic response to cyclic loading.

---

1 Introduction

Poroelasticity, though classically derived in the context of soil consolidation,^1–4 has found applications far beyond this, covering areas as varied as biological tissues, seismology, transport and hydrogeology, to name a few. Their common denominator is the coupling of fluid flow and deformation of a porous solid. Linear poroelasticity theory was first formalised by Biot⁴ and later extended to large deformations,⁵ laying the groundwork still used today for the study of soft porous materials. Detournay and Cheng⁶ provide a good summary of the origin of poroelasticity and developments following on from Biot's work.

Many poroelastic materials are subject to cyclic loads or deformation throughout their lifetime, ranging from biological to geological systems, with applications in areas such as medicine, bioengineering, construction engineering, climate and drug delivery. For example, soil may experience periodic loading due to passing traffic,^7,8 seismic activity (particularly relevant in construction engineering),^9–11 and waterwaves or tidal forcings on the seabed.^8,12–15 Biological soft tissues such as tendon and cartilage also experience periodic tensile and/or compressive loads during daily activities such as walking, swimming, running or playing an instrument.^16–18 Extensive research supports the hypothesis that tendon (and cartilage) structure is remodelled in response to cyclic mechanical loading,^19–22 similarly to hard tissue such as bone, where naturally occurring cyclic loading has been shown to play an important role in bone remodelling and healing.^23–27 In addition, there is a rich literature surrounding solute transport in response to cyclic loading of these biological soft tissues, where the effect of parameters such as loading amplitude, frequency, and size of solute have been investigated.^23,27–34 Brain tissue is commonly modelled as a poroelastic material which may undergo pulsating loads from the dilation of arterioles; in this area there is a particular focus on fluid movement and the clearing of metabolic waste.^35–37 The cytoplasm of living cells has also been shown to behave as a poroelastic material.³⁸ A model system for many of these materials are gels, which are typically easier to test experimentally and from which elastic properties, such as stiffness, Poisson's ratio and permeability, may be extracted in the context of poroelastic theory.^39–44

Although many poroelastic models assume the material to be uniform, poroelastic media commonly exhibit spatially heterogeneous material properties, such as non-uniform stiffness, permeability and Poisson's ratio. When modelling soil, heterogeneity is typically accounted for by introducing horizontal layers of different thickness, density or permeability, resulting in discrete jumps in the material properties between each layer.^9,10,13 Trefry et al.¹⁴ explore the effect of continuous spatial heterogeneity in groundwater systems undergoing cyclic loads (tidal forcings) by imposing a randomly assigned heterogeneous hydraulic conductivity, and find this leads to chaotic mixing. Concerning biological tissues, Zhang²⁸ investigates the effect of linear variation of stiffness and permeability, chosen to fit experimental data for cartilage depth-dependent material properties, in a cylindrical cartilage disc subject to cyclic compression. Similarly to the layered soil models, Kameo et al.⁴⁵ and Ruiz, et al.⁴⁶ assign unique material properties to vertical layers of trabecula (anatomical unit of bone) and intervertebral disc (IVD) respectively. Hydrogels also exhibit heterogeneity originating from their complex polymer structures, and various experiments explore the impact of this heterogeneity on network remodelling and bulk material properties such as fracture strength.^47–52 While these models all find that the various forms of heterogeneity significantly impact the response of the material, to date there is no systematic study investigating the effect of generic continuous material heterogeneity on the poroelastic response to cyclic loading. Additionally, the biological models consider material heterogeneity to be a characteristic of the healthy tissue, however in tissues such as tendon several studies have reported altered mechanical and material properties in the damaged tendon compared to its healthy counterpart.^53,54 Damaged tissue has primarily been explored as cracks or fractures in the material,^55,56 as in rock mechanics,⁵⁷ rather than modified material properties. In this study, we will refer to “damage” as heterogeneity in the material properties of the tissue.

Our study is particularly motivated by tendinopathy, a common and painful condition affecting tendons and characterised by a modified structure of the tissue.¹⁶ Modelling the mechanics of such a complex structure forms a formidable challenge, for which a number of distinct approaches exist. Tendons have been modelled as purely elastic (these models are usually fibre reinforced),^58–60 viscoelastic,⁶¹ poroelastic,^22,62–67 or even poroviscoelastic.⁶⁸ Some models have also addressed the evolution of material properties in tendon. For example Wren et al.⁶² and Scott⁶⁷ incorporated phenomenological rules to update global tissue permeability or stiffness, and Lavagnino et al.²² built a submodel for the tendon cell to explore the relation between physical stimuli and protein production (mechanotransduction). A number of these models are finite element models; while they enable additional levels of biological detail to be captured,^{22,58,63–66,68} this additional complexity can cloud identification of the underlying driving physical mechanisms governing the response of the system. Tendons are able to withstand large strains,⁶⁹ motivating a non-linear framework; a one-dimensional formulation, as we will present here, allows for non-linear kinematics whilst remaining simple enough to draw physical and intuitive understanding of the underlying mechanics.

In this paper, we present a non-linear poroelastic model for a heterogeneous material subject to a uniaxial tensile cyclic stress or displacement of varying frequency. Of particular relevance to our study is the work of Fiori et al.,³³ and the thesis work of Fiori.⁷⁰ In, ref. 33 a systematic study was conducted exploring the effect of compressive cyclic applied displacement amplitude and frequency on the fluid response of a homogeneous poroelastic material in a uniaxial setting. They showed that as the period decreases and the amplitude increases the root-mean-square relative error between the linear and non-linear model grows, identifying the validity of the linear poroelastic model as , where A and T are the dimensionless amplitude and period of applied displacement, such that even at small amplitude A linear poroelasticity does not provide a good approximation for high frequency loading (small T). In addition, the authors found that decreasing loading period increases localisation of the deformation while varying loading amplitude and initial porosity has little qualitative effect on the poroelastic response. The results from Fiori's study for the fluid response of a homogeneous poroelastic material to cyclic compression qualitatively agree with the results from the uniform version of our model, despite employing a different solid constitutive law, though we might expect this choice to play a greater role when heterogeneity is present. The aim of the present paper is to systematically investigate and intuitively understand the role of localised material heterogeneity on the response of a poroelastic material to cyclic loading. This builds on the thesis work of Fiori⁷⁰ (Ch. 5), which considered the impact of layered permeability in the uniaxial setting, including the impact on solute transport, for applications in brain and cartilage. Our objective is to explore the effect of different forms of heterogeneity, and localised to different positions in the material. We also analyse how the consequences of heterogeneity vary with loading frequency, and we uncover both quantitative and qualitative differences in the material response to an applied displacement versus an applied stress.

Following MacMinn et al.⁷¹ and Coussy,⁷² we introduce a non-linear poroelastic model which in its reference state has uniform porosity and position-dependent permeability, and non-linear elasticity with position-dependent stiffness. Following Fiori,⁷⁰ which presents a similar set of equations as we derive below, we employ a Lagrangian framework that better facilitates incorporation of material heterogeneity. Assuming uniaxial deformation, we apply a tensile cyclic stress or displacement and fixed pressure at one end with no flow and displacement through the other end, and solve the full problem numerically using a finite volume method and stiff ode solver. We first impose damage as a local decrease in stiffness, and investigate the effect of this heterogeneity on the strain and flux response of the material for a fixed loading frequency and amplitude. We then characterise how the response varies with the magnitude and location of the decrease, as well as with loading frequency. We also consider the effect of a local decrease in permeability and compare to our results for heterogeneous stiffness. Finally, we discuss the physical interpretation of the model, including relevance of our findings to tendon and how our mechanistic insights might contribute to our understanding of tendinopathy development or diagnosis.

2 Model

We consider a fully saturated poroelastic material composed of a solid and a fluid phase, such that ϕ_f + ϕ_s = 1 where ϕ_f,s are the true fluid and solid volume fractions respectively. Here the word “true” is used interchangeably with “Eulerian”, meaning a quantity measured with respect to the current configuration. This is contrasted with “nominal” or “Lagrangian” quantities, which are measured with respect to the reference (initial) configuration. For example, the true porosity ϕ_f measures the current fluid volume fraction of the total current volume dv, while the nominal porosity Φ measures the current fluid volume fraction of the total reference volume dV₀, where capital letters are used to denote Lagrangian variables. The true and nominal porosities are related by Φ_f,s = Jϕ_f,s where J is the Jacobian determinant which measures local volume change i.e. dv = JdV₀, and is defined later in eqn (3). For the fully saturated poroelastic material in a Lagrangian description, we have Φ_f + Φ_s = J. Going forward, we drop subscripts and use only Φ ≡ Φ_f (ϕ ≡ ϕ_f) to denote nominal (true) porosity. Supposing that the material begins at time t = 0 in a stress free state, the reference or initial configuration refers to the material's configuration at time t = 0, while the current configuration refers to the material's configuration at current time t.

Since the fluid flow and solid deformation are coupled, we are presented with a choice of framework: Lagrangian or Eulerian. Most poroelastic descriptions are presented in an Eulerian frame, as this is most natural for fluid flow, but it is more natural to consider material heterogeneity in a Lagrangian frame, as this gives immediate access to material points and their properties. While we restrict to fixed heterogeneity in this paper, a Lagrangian frame could also more readily accommodate complexities such as tissue remodelling, via an evolution of properties at material points.

We therefore consider the uniaxial deformation of a poroelastic material in the Lagrangian domain

0 ≤ Z ≤ L₀ (1)

where L₀ is the length of the material, as shown in Fig. 1. We assume the boundaries of the system coincide with the boundaries of the solid. The material is pulled at Z = 0 where it is also subject to ambient fluid pressure p_A. The upper boundary Z = L₀ is fixed with no deformation or flow.

Fig. 1 We consider a poroelastic material of length L₀ which is pulled at Z = 0 (red arrow) where it is also subject to ambient fluid pressure p_A, and fixed at Z = L₀ with no deformation or flow. The blue arrows represent fluid flow. The material is shown in the reference configuration.

2.1 Kinematics

The deformation gradient tensor F is given bywhere X and x denote Lagrangian and Eulerian coordinates respectively, and x_i and X_j are components with i, j = 1, 2, 3. The Jacobian determinant is defined by

J = det(F). (3)

For uniaxial deformation, the deformation gradient tensor is then given by

The solid displacement U is defined as the difference between the current configuration x(X,t) and the reference (initial) configuration X of the solid:

U = x(X,t) − X, (5)

where x(X,t) maps the reference state to the current state. The material begins at t = 0 in the reference state, i.e. U(X,0) = 0. The deformation gradient tensor can then be expressed as

F = ∇_X(x) = I + ∇_X(U) (6)

where ∇_X denotes the Lagrangian gradient operator in Lagrangian coordinates (∂/∂X_I) and I is the identity tensor. For future use, ∇_x denotes the Eulerian gradient operator (∂/∂x_i). For uniaxial deformation then U = (0,0,U) so that strain

J = 1 + U_Z (7)

where X = (X,Y,Z) and we have used subscript Z as a shorthand for ∂/∂Z.

2.2 Continuity

We next derive local conservation equations to describe the deformation. It is more natural to write these first in terms of Eulerian quantities, which we then convert into a Lagrangian frame. We consider both the solid grains which form the solid skeleton and the fluid to be incompressible, so that when the poroelastic material is compressed the pore size must decrease (fluid is squeezed out), i.e. the solid fraction increases as the total volume decreases and vice versa. The fully saturated poroelastic material is therefore compressible through pore rearrangement. Mathematically, the incompressibility condition for the solid phase is expressed as (1 − ϕ)dv = (1 − ϕ₀)dV₀ or (J − Φ)dV₀ = (1 − Φ₀)dV₀ where ϕ₀ = Φ₀ is the initial porosity. The change in volume is thenwhere we have defined the normalised quantity [capital Phi, Greek, circumflex] ≡ Φ − Φ₀. Eqn (7) then becomes:

U_Z = [capital Phi, Greek, circumflex]. (9)

Since the fluid and solid are incompressible, and we assume uniform density for both phases, the fluid and solid mass densities are constant. Local continuity for the solid and fluid phases is then written

where v(t) is an arbitrary current volume, and D^f/Dt and D^s/Dt are respectively the material derivatives with respect to the fluid and solid (i.e. the time derivative that an observer attached to a fluid or a solid particle would derive⁷²).

Let v^s,f ≡ V^s,f be the solid and fluid velocities respectively. Note that although we write velocities (and subsequently pressure) in lowercase for Eulerian and uppercase for Lagrangian, they are the same in both of these frames – only the coordinate at which they are evaluated changes. The continuity eqn (10) can be rewritten in terms of arbitrary reference volume V₀ and using Euler's identity D^f,sJ/Dt = J∇_x·v^f,s:

where we have also defined the Eulerian relative flow q = ϕ(v^f − v^s). From Nanson's formula (see Appendix A) the Lagrangian relative flow Q is equal to Eulerian relative flow:where we used Φ = Jϕ. Eqn (11) are converted into Lagrangian coordinates using the deformation gradient tensor (4) to arrive at the Lagrangian equations for fluid conservation of mass as follows:

∇_X·Q + J(F⁻¹∇_X)·V^s = 0 (13b)

where we have denoted D^s/Dt ≡ d/dt. For uniaxial deformation we have V^s,f = (0,0,V^s,f) and Q = (0,0,Q), and the continuity eqn (13) simplify to:Note that the solid continuity eqn (14b) can also be replaced by eqn (8), as the incompressibility condition automatically satisfies conservation of solid mass. The continuity formulation (14b), however, will be useful when implementing boundary conditions.

2.3 Constitutive law for fluid flow

Let p(x,t) ≡ P(X,t) be the fluid pore pressure. We assume Darcy's law to govern the relative flow of the fluid with respect to the solid, neglecting gravity. In an Eulerian frame it is given bywhere k is the permeability tensor and μ is dynamic viscosity. The Lagrangian flux Q is thenLetting k denote permeability in the Z-direction, then under uniaxial deformation the flux is given byWe assume permeability is porosity dependent and adopt the normalised Kozeny–Carman formula,⁷¹ widely used for a range of materials. This formula readswhich we express in an Eulerian frame since it is physically more meaningful, and k₀ is the reference permeability value. In particular, it captures two important physical limits: as the true porosity vanishes (ϕ → 0), permeability vanishes, and as the true porosity tends to 1 (ϕ → 1), permeability diverges. This prevents the flow from driving true porosity below zero and above 1. An important consequence of this law is that permeability increases with porosity, which is consistent with observations in tendons.^66,73 In Lagrangian variables, and with the expression (8) for J, the formula reads

2.4 Mechanical equilibrium

Following MacMinn,⁷¹ the true total stress σ, which describes the current total force supported by the fluid and solid per current unit total area, can be divided into the true fluid stress σ^s and true solid stress σ^f (fluid and solid force per unit fluid and solid area respectively), assuming the phase area and phase volume fractions are equivalent: σ = (1 − ϕ)σ^s + ϕσ^f. We can write σ^f = −(p − p_A)I^71,72 where p_A is the ambient pressure as introduced at the start of Section 2. Since the fluid is coupled to the solid skeleton, the fluid exerts at all times a pressure on the solid, so that the total solid stress includes this component even when the fluid is at rest. This component of the stress cannot contribute to the deformation of the solid skeleton of incompressible grains, motivating the use of Terzaghi's effective stress σ′ ≡ (1 − ϕ)(σ^s + (p − p_A)I). The true total stress is then

σ = σ′ − (p − p_A)I. (20)

To satisfy mechanical equilibrium, the solid skeleton and the fluid must jointly support the mechanical stress. In the absence of body forces and inertia, this requires

∇_x·σ = ∇_x·σ′ − ∇_xp = 0 (21)

in an Eulerian frame. To describe mechanical equilibrium in the Lagrangian solid frame, let s be the nominal total stress, also known as the first Piola–Kirchoff stress tensor. It is related to the Cauchy stress tensor σ by

s = JσF^−T, (22)

and the effective nominal stress s′ is then s = s′ − JF^−TP. Combining (21) and (22), we arrive at the following equation for mechanical equilibrium in the Lagrangian frame:

∇_X·(s′) = ∇_X(JF^−TP). (23)

Under uniaxial deformation this may further be simplified towhere s′ denotes the ZZ-component of the stress tensor.

2.5 Constitutive law for the solid skeleton

Incompressibility of the solid grains means that the volume of solid contained within the Lagrangian element is fixed and cannot change, however the current volume of the material element can change freely through fluid flow and pore rearrangement. The constitutive law for the effective solid stress describes the “drained” behaviour of the solid skeleton – this corresponds to what would be measured for the “dry” skeleton, or if the “wet” skeleton were tested sufficiently slowly so that the fluid is stationary. Since the network forming the solid skeleton is compressible, a compressible constitutive law must be considered. Some care must be taken in choosing an appropriate effective stress law for a poroelastic solid, as highlighted by Zheng and Zhang,⁷⁴ who find that only the Green-strain tensor is compatible with the restrictions of thermodynamics. Accordingly, we consider a neo-Hookean model for a compressible solid:where W is the strain energy density, M is the oedometric stiffness and Λ is Lamé's second parameter, related to M by Poisson's ratio νAs a form of incorporating material heterogeneity, we consider the general case where M and Λ may be functions of position. The ZZ-component of effective stress is thus given byNote that for ν = 0.5 we have s′ = M(J − 1) = MU_Z, which is linear.

Eqn (8), (14a), (17), (19), (24) and (27) with specified boundary and initial conditions, may be combined to form a closed model for the evolution of porosity Φ in a Lagrangian framework, as will be highlighted in the summary 2.10.

2.6 Boundary and initial conditions

2.6.1 Lower boundary Z = 0. We subject the poroelastic material to cyclic tension via applying either a periodic stress or periodic displacement boundary condition at Z = 0. Adopting the convention that tension is positive and compression is negative, the tensile applied stress (AS) takes the formwhile the tensile applied displacement (AD) is given byHere A_s,d are the magnitudes of the applied stress and applied displacement, respectively, and ω is the frequency. Both scenarios generate tensile stress, since s*(t) ≥ 0 and a(t) ≤ 0 for all time. Note that these conditions are not equivalent: applying a stress with amplitude A_s does not result in a displacement of amplitude A_s at the boundary, and vice versa. We also define the loading period T = 2π/ω. The associated boundary conditions are

s(Z = 0,t) = s*(t) (29a)

or U(Z = 0,t) = a(t). (29b)

For an applied displacement the boundary condition can also be expressed as V^s(Z = 0,t) = da/dt ≡ ȧ (overdot denotes derivative with respect to time). It will be useful to differentiate between loading time periods when the material is being “pulled”, characterised by ṡ* > 0 or ȧ < 0, versus the unloading or relaxing time periods when the material is “let go” or “pushed” back, for which ṡ* < 0 or ȧ > 0.

In terms of the fluid, we suppose that the Z = 0 boundary is permeable and held at a fixed pressure, that is

P(Z = 0,t) = p_A. (30)

2.6.2 Upper boundary Z = L₀. The Z = L₀ boundary is assumed fixed and impermeable to fluid, so that:

U(Z = L₀,t) = 0, V^f(Z = L₀,t) = 0. (31)

Note that this also implies that V^s(Z = L₀,t) = 0 and so Q(L₀,t) = 0. Combined with solid continuity (14b) we find that

V^s = −Q (32)

everywhere, and from (12)so that the fluid and the solid always locally move in opposite directions, since J − Φ = 1 − Φ₀ > 0 from (8).

2.6.3 Initial condition. The material starts in the reference configuration, corresponding to a relaxed or undeformed state, with initial uniform porosity Φ₀

U(Z,0) = 0 and Φ(Z,t) = Φ₀. (34)

2.7 Non-dimensionalisation

We non-dimensionalise as follows:where k₀ and M₀ are the scalar reference permeability and oedometric stiffness and T_pe = μL₀²/k₀M₀ is the classical poroelastic timescale which corresponds to the characteristic diffusion timescale of the material, i.e. the time it takes for the response to an applied stress/displacement at Z = 0 to reach Z = L₀. The stiffer, shorter and more permeable the material, the quicker the response diffuses through the material. For a constant applied stress or displacement the poroelastic timescale corresponds to the time it takes for the system to reach steady state. For a cyclic applied stress/displacement, the loading frequency [small omega, Greek, tilde] falls into two different regimes: [small omega, Greek, tilde] ≤ 2π, for which the diffusion time is less than or equal to the loading period, and [small omega, Greek, tilde] > 2π, for which the diffusion time is greater than the loading period. The smaller [small omega, Greek, tilde] is, the more time the stress imposed at Z = 0 has to diffuse through the material, homogenising the response. On the other hand, the greater [small omega, Greek, tilde] is, the more localised the response becomes to Z = 0. These regimes may also be characterised by the Péclet number, as highlighted in Appendix C. Note that varying [small omega, Greek, tilde] may correspond to varying ω, L₀, k₀/μ or M₀. Going forward we consider the non-dimensional problem and drop the tilde notation, unless specified otherwise.

2.8 Parameter values

Our choice of parameter values is motivated by tendon, however these can easily be modified for other materials or tissues. Reported values for tendon are variable, with particular disagreement on the Poisson's ratio.^75–77 We choose the following:^22,78,79where M₀ is calculated from Young's modulus E₀ using

Tendons, in particular energy-storing tendons, undergo large deformation, with strain ranging from 5% to as large as 30%.⁶⁹ As we will discuss later, the system is diffusive so that for homogeneous stiffness strain U_Z decreases with Z: to bound strain throughout the material we simply need to bound strain at Z = 0. Setting 0.05 ≤ U_Z(0,t) ≤ 0.3 and recalling that J = 1 + U_Z, we calculate the corresponding neo-Hookean stress using eqn (27) with the above parameters (and uniform stiffness M = 1) to find an approximate range for A_s: 0.05 ≤ A_s ≤ 0.3. Unsurprisingly, A_s ≈ U_Z(0,t), since for ν = 0.3 we have s′ ≈ MU_Z for U_Z < 0.5. For an applied displacement, we look for a range for U(0,t) = a(t). Since strain decreases with Z we have . Accordingly, the range for U(0,t) is smaller than the range for U_Z(0,t) so we set 0.05 ≤ A_d ≤ 0.2. Fiori et al.³³ examined the impact of loading amplitude for a material subject to cyclic compression and demonstrated it had little to no qualitative effect on the poroelastic response. Accordingly, we do not expect loading amplitude to qualitatively affect the poroelastic response of the material subject to cyclic tension. In this study, we will therefore vary loading frequency and fix loading amplitude:

A_s = 0.2 and A_d = 0.1. (38)

While the quantitative results will be tied to these specific values of amplitude, we have confirmed that the results presented are qualitatively similar at the lower end of our amplitude (results not shown).

Since we are motivated by tendons, we are interested in loading frequencies relevant to physical activities that engage energy-storing tendons such as the Achilles tendon, like running, walking or resistance training. Given T_pe ≈ 14 s, resistance training corresponds to [small omega, Greek, tilde] ≤ 2⁸⁰ and running corresponds to ω ≥ 20,⁸¹ leaving 2 < [small omega, Greek, tilde] < 20 for activities in between such as walking or slow jogging. We will first look at the intermediate frequency [small omega, Greek, tilde] = 10 as a baseline frequency, before varying frequency in the range [0.5, 50].

2.9 Material heterogeneity

Our primary interest in this paper is the effect of material heterogeneity on the response to cyclic stress or displacement driven loading. Motivated by observed material changes in damaged biological tissues, where for instance diseased tendon exhibits a localised decrease in stiffness, we consider localised decreases in either stiffness or permeability. We impose heterogeneity in stiffness or permeability via a function f(Z) such that

M(Z) = f(Z) and k(Φ,Z) = k_KC(Φ)f(Z). (39)

We consider the following form for f:where l and d are the location and magnitude of damage such that f(l) = 1 − d, as shown in Fig. 2. We set the variance to c = 1/10. Conversely, an increase may also be considered by taking 2 − f(Z). Although our analysis will be focused on this form for f, the framework can easily accommodate for other functional forms of heterogeneity.

Fig. 2 Material damage is modelled with an inverted Gaussian f of magnitude d at Z = l.

2.10 Summary

Expressed in terms of [capital Phi, Greek, circumflex], the system evolves according to a non-linear advection-diffusion equation,wherewhich respectively represent the bulk advection and bulk diffusivity of the material. The advection term is only present when stiffness M is heterogeneous. The effective stress is given bywhere stiffness for the homogeneous and heterogeneous cases is respectively given byand similarly permeability is given by

k = k_KC or k = k_KCf(Z) (45)

with permeability law

The boundary condition on [capital Phi, Greek, circumflex] at Z = 0 is

where Φ* is calculated by inverting the relationship s′([capital Phi, Greek, circumflex]) for given s*. The boundary condition on [capital Phi, Greek, circumflex] at Z = 1 isFinally, we impose the initial condition

[capital Phi, Greek, circumflex](Z,0) = 0. (49)

As discussed, working in a Lagrangian framework (with respect to the solid) is the natural choice when considering material heterogeneity. This also presents a significant advantage when solving the system, as the boundaries remain stationary so that the equations can be solved on a fixed domain, unlike in an Eulerian framework which would involve a moving boundary. We solve the system (41)–(49) numerically for [capital Phi, Greek, circumflex] using a finite volume method in space and an implicit variable-step size method in time via the function ode15s in MATLAB (see Appendix B for details).

3 Results

We first consider heterogeneous stiffness for a fixed loading frequency ω and amplitude A_s,d. We compare the response of the system for a fixed damage magnitude and three damage locations to the response of the undamaged system over one steady cycle. By steady cycle we mean the code was run for sufficiently long so that the transient behaviour decays and the response has settled down to a periodic state. As shown by Fiori et al.,³³ under AD the transient response decays exponentially with time, and the time of decay is independent of ω. From our numerical simulations we visually determined the decay of the transient phase, as shown in Appendix C, and found the same is true under AS. Since the time it takes for the transient response to decay is independent of ω, the number of cycles in the transient phase increases with ω. We next define a metric to capture the effect of damage, allowing us to explore a range of damage magnitudes, damage locations and loading frequencies. We also briefly explore the effect of heterogeneous permeability in comparison to heterogeneous stiffness. Note that we will refer to normalised porosity as porosity for brevity.

3.1 Heterogeneous stiffness

Consider heterogeneous stiffness M = f(Z) and Λ = νf(Z)/(1 − ν). We only consider a local decrease in stiffness here, however we also display results for a local increase in stiffness in Appendix D.

3.1.1 Full steady response. In Fig. 3 we plot the strain U_Z(Z,t_i) and flux Q(Z,t_i) profiles, at 9 time points, equally spaced over one period, under an applied stress (AS) with (A_s, ω) = (0.2, 10) in the first two rows and under an applied displacement (AD) with (A_d, ω) = (0.1, 10) in the last two rows. We differentiate between the loading phase, ṡ* > 0 or ȧ < 0 (dark blue to light blue) and unloading phase, ṡ* < 0 or ȧ > 0 (light red to dark red), with the end of the cycle in dotted red. The first column is the homogeneous case, and columns 2–4 display profiles with imposed local decrease in stiffness where the location of the dip is increased from left to right (moving away from the moving boundary). Furthermore, A_s,d and are of the same order, and hence we are in a regime where we do not anticipate the linear model to provide an accurate approximation.³³ For ω = 10 the transient phase decays after about 1 cycle.

Fig. 3 Strain (first and third rows) and flux (second and fourth rows) plotted for 9 equally spaced values of the 20th cycle (38π/ω ≤ t ≤ 40π/ω), for uniform stiffness (first column) and locally decreased stiffness with d = 0.35 (columns 2–4). The first two rows are in response to an applied stress (A_s = 0.2, ω = 10) and the second two rows are in response to an applied displacement (A_d = 0.1, ω = 10). The applied stress s*(t) and displacement a(t), and their respective time derivatives ṡ*(t), ȧ(t) are displayed as insets in the first column. We differentiate between the loading phase (ṡ* > 0, ȧ < 0, dark blue to light blue) and unloading phase (ṡ* < 0, ȧ > 0, light red to dark red). The end of the cycle is in dotted red.

The governing eqn (41) is a non-linear advection-diffusion equation (the advection term is zero for homogeneous stiffness) for porosity, or strain since U_Z = [capital Phi, Greek, circumflex]; due to the constraint of uniaxial deformation and solid incompressibility, any increase in strain results in an increase in pore space. We refer to strain and porosity interchangeably in the discussion that follows.

A key feature of the non-linear model is porosity-dependent permeability: with the Kozeny–Carman expression for permeability (46) and neo-Hookean elasticity (43), the diffusion coefficient (42) takes the form

which is an increasing function of porosity, so long as Λ < M (and [capital Phi, Greek, circumflex] > −Φ₀ i.e. Φ > 0 which is always true). The Kozeny–Carman term k_KC increases with porosity faster than the stress term J⁻¹∂s′/∂[capital Phi, Greek, circumflex] decreases with porosity, so that diffusion as a whole is increased as porosity increases. Note that if we had constant k = 1, the diffusion coefficient would be a decreasing function of porosity.

The implications of the porosity-dependent permeability can be seen in Fig. 3. As the material is pulled (loading, blue lines, ṡ* > 0, ȧ < 0), fluid starts to enter (positive flow), although this is not an instantaneous response, to satisfy the incompressibility condition (rows 2 and 4 of Fig. 3), increasing porosity (rows 1 and 3) which increases permeability and thus diffusivity. The increased diffusivity means the stress imposed by AS or generated by AD at Z = 0, and the associated strain, diffuses farther into the material over the cycle than it would if permeability did not depend on porosity. This also causes porosity to increase further, creating a positive feedback loop where increased porosity increases permeability which increases diffusivity which increases porosity and so on. As tensile stress decreases, or the material is pushed back (unloading, red lines, ṡ* < 0, ȧ > 0), squeezing fluid out, diffusivity decreases (but not below 1), slowing the response and diffusion of strain/porosity through the material. Naturally, strain still decays to some extent as Z → 1, the extent of which depends on ω, however it would decay much faster if permeability k were constant. In addition, flux decays with Z and is zero at the upper boundary Z = 1 as required by eqn (48) (rows 2 and 4).

Under AS, strain/porosity at the lower boundary Z = 0 returns to zero at the end of the cycle, as depicted by the red dotted line in the first row of Fig. 3. This is because s′ ≈ M[capital Phi, Greek, circumflex]; in particular, this relation is exact for s′ = 0, i.e. when stress is zero, strain/porosity is zero. For Z > 0 however, [capital Phi, Greek, circumflex] > 0, indicating stress is non-zero, even at the end of the cycle. This is due to the diffusive nature of the system: the stress acquired as Z → 1 is not instantaneously lost. Additionally, whether the imposed stress at Z = 0 eventually reaches the material points as Z → 1 will depend on ω. Since the solid phase is incompressible, [capital Phi, Greek, circumflex] > 0 indicates an increase in volume of the material; the fluid that entered during the transient period (initial cycles) is trapped for the duration of loading, and will only exit again once the loading stops and the material can fully relax. For our 1D system, an increase in volume is equivalent to a negative displacement at Z = 0, so we can also write this as U(0,t) < 0 for Z > 0.

Under AD on the other hand (third row), we observe a negative strain near Z = 0 at the end of each cycle (red dotted line), highlighting the compressive stress required to push the material back to its initial position (displacement U(0,t) = 0 for t = 2nπ/ω, integer n, from eqn (29b)). The strain/porosity at the end of the cycle is also non-zero throughout the rest of the material, indicating that the material does not return to its initial configuration here either, despite returning to its initial volume. This is possible via a re-distribution of porosity throughout the material: we see that porosity is decreased near Z = 0 and increased near Z = 1 (such that total porosity as required). Again the diffusive nature of the system is at play: the resulting stress/strain from the imposed displacement is not instantaneously transferred to the rest of the material, and how far it will penetrate the material over the cycle depends on ω. The smaller ω, the less compressive strain will be generated near Z = 0 at the end of the cycle and the more homogeneous the strain will be in Z. For ω = 10 here, an intermediate frequency, the response is not strongly localised to Z = 0 but is not uniform in Z either.

As expected under AS (first row of Fig. 3), an overall decrease in stiffness results in greater strains, as the material is overall weaker, implying greater displacements. Strain is particularly increased at the point of maximum damage. For example, with damage imposed at l = 0.25 (second column), we find a local increase in strain corresponding to the local decrease in stiffness, and exhibiting a maximum at Z = 0.25. The flux profiles (second row) are much less affected by the heterogeneity, showing only a small increase around Z = l.

To understand how porosity can be increased whilst flux is negligibly affected, consider the expression (12) for flux, and let δΦ be a small increase in porosity. We can assume δΦ/J is small since δΦ ⪅ 0.1 here and J > 1. Let J_δ = J + δΦ and V^f,s_δ be the new fluid and solid velocities. We find Q_δ = (V^f_δ − V^s_δ)(Φ + δΦ)/J_δ = (V^f_δ − V^s_δ)Φ/J + O(δΦ/J). Assuming V^f,s_δ ≈ V^f,s then Q_δ ≈ Q. Since the velocity V^f,s is unchanged, this equates to the ratio of porosity to volume change staying approximately constant for small δΦ/J. This approximation breaks down as δΦ and V^f,s_δ grow, but for the small values dealt with here it is reasonable. Moreover, we note that this is Lagrangian flux, and the effect on true, Eulerian, flux q will be greater since q = JQ, i.e. q_δ ≈ q + QδΦ.

Under AD on the other hand, we find that (for these parameters) strain, and by association displacement, is unaffected by the local decrease in stiffness (third row), whilst flux, and by association stress, is heavily impacted (fourth row). As evident in the last row, the flow is locally perturbed around the damage location, characterised by a local increase or decrease in inflow or outflow either side of Z = l compared to the homogeneous case. For example for l = 0.25 (2nd column), inflow is increased to the left of l and decreased (or turned to outflow) to the right of l for most of the cycle, however towards the end of the cycle outflow is decreased (or turned to inflow) to the left of l and increased to the right of l. For l = 0.75 on the other hand (4th column), inflow is always increased to the left of l and outflow is always increased to the right of l. For all cases, the maximum magnitude of this increase in inflow is located at l − c and the maximum magnitude of increase in outflow is at l + c, where c is the variance of the Gaussian defined in (40). Note that flux is unchanged at Z = 0, as required by (32) and (47) under AD.

To understand the effect on flux under AD, consider now the expression (17) for Q. Flux is governed by stress gradients (recall ∂s′/∂Z = ∂P/∂Z) and permeability and porosity of the material (k/(1 + [capital Phi, Greek, circumflex])). Under an applied displacement, the position of the material at Z = 0 is imposed, irrespective of its stiffness (and thus the flux at Z = 0 too). If the material is locally weakened, a smaller stress is required to achieve the associated displacement throughout the material, and we see a local dip in the stress profile in the vicinity of the damage. This modifies the stress and thus pressure gradient, which modifies the flux profile accordingly, as illustrated in Fig. 4 which plots the time integrated stress and time integrated flux for the homogeneous case in yellow and the damaged cases in orange against Z, where t_f = 40π/ω. Since [capital Phi, Greek, circumflex] ≈ s′/M and the ratio s′/M remains approximately constant, strain/porosity is negligibly affected.

Fig. 4 Time integrated stress and flux response for uniform stiffness (light orange) and locally decreased stiffness (dark orange) against Z.

Additionally, we find that under both AS and AD the maximum strain/flux in the damaged area decreases with increasing l, however the magnitude of strain/flux in the damaged area is also more uniform in time for greater l (rows 1 and 4). For example under AS, for l = 0.25, the magnitude of strain at l varies between 0.08 and 0.27, whereas for l = 0.75, it only varies between 0.13 and 0.21. Similarly under AD, for l = 0.25, the magnitude of flux at l − c varies between −0.35 and 0.73, whereas for l = 0.75, it only varies between −0.25 and 0.06. This is once again due to the diffusive nature of the system which localises the stress–strain response to Z = 0 where the boundary condition is applied, so that stress and strain vary more near Z = 0 and less as Z → 1.

3.1.2 Defining a metric. The above results highlight that heterogeneous stiffness impacts the strain and flux of the material under AS/AD both temporally and spatially. To characterise the system response in a more concise form, in this section we define and analyse several metrics that capture the net response. To this end, we first define the following time integrated quantities:Note this definition includes the transient phase, and integrates up to a fixed time t_f. We choose to set t_f = 40π/ω which corresponds to 20 cycles. By integrating the modulus, these metrics do not differentiate between compressive and tensile strain, or between inflow and outflow. Under AS, the absolute value has no effect on strain since U_Z ≥ 0. Under AD on the other hand the compressive strain now has an additive effect with tensile strain. In both cases, outflow and inflow are also additive. These metrics are therefore cumulative, and the more cycles are applied, the more the effect of damage on strain and flux will cumulate, increasing the difference with the undamaged case. We refer to as “cumulative strain” and as “cumulative flux”.

The cumulative strain and flux are plotted in Fig. 5 for d = 0.35 and three locations l (increasing l from dark to light orange) as well as the uniform case in dotted grey. The top left panel highlights the increase in strain about Z = l under AS, and the bottom right panel highlights the disturbed flux about Z = l under AD. The figure additionally highlights the negligible effect heterogeneous stiffness has on strain under AD (bottom left) and the small effect on flux under AS (top right).

Fig. 5 Cumulative strain (left) and flux (right) against Z under AS (A_s = 0.2, ω = 10, top row) and AD (A_d = 0.1, ω = 10, bottom row) over 20 cycles as defined in (51), for homogeneous stiffness (dotted grey line) and for damaged stiffness with d = 0.35 and varying location l = 0.25, 0.5, 0.75.

Clearly, the spatial features present in the strain and flux fields in Fig. 3 are not necessarily preserved in the cumulative flux and strain profiles. For example, for l = 0.25, is decreased to the right of Z = l. For l = 0.5 and l = 0.75, the cumulative flux to the right of l is first decreased and then increased. Additionally, the cumulative strain under AS exhibits the same profile and magnitude for each value of l.

A natural extension to this spatial metric is to integrate over Z as well:

where the time integral captures the cumulative effect of damage on strain or flux magnitude as before, and the spatial integral averages the cumulative strain/flux over space. Again, we choose t_f = 40π/ω i.e. we integrate over 20 cycles. We refer to these quantities as “net strain” and “net flux” respectively. In Fig. 5 this simply corresponds to the area under the and curves. For this particular case, we find that although the maximum strain/flux about Z = l decreases with l, the net strain is approximately equal for all l under AS and the net flux slightly increases with l under AD. As these metrics output a single number for any given form of heterogeneity and loading conditions, they enable to more readily explore the system response to variation in parameters, which we do in the following subsection.

3.1.3 Parameter variation. In order to more thoroughly investigate the effect of damage magnitude d, damage location l, and loading frequency ω, we continuously vary d and ω for three values of l, computing the net strain and flux for each parameter set. We then compute the difference between the damaged and undamaged net strain/flux viawhere ·_d and ·_u denote “damaged” and “undamaged” respectively. If ΔU, ΔQ > 0 then damage increases net strain/flux, and conversely if ΔU, ΔQ < 0 then damage decreases net strain/flux. An increase in |ΔU| or |ΔQ| means an increase in the effect of damage.

Unsurprisingly, increasing damage magnitude d increases the effect of damage. For example the difference in net flux ΔQ under AD increases with d, as illustrated in Fig. 6(a). Here, the dependence of net flux on damage location l is not the same for all d: for d ≤ 0.5 net flux increases with l, in line with our observations from the previous section, whilst for d ≥ 0.5 the relationship reverses such that net flux decreases with l. This can be understood by noting that as d increases, the decrease with l in maximum flux at l ± c steepens, such that beyond d = 0.5 the net flux also decreases with l. In most cases however ΔU and ΔQ decrease with l, as illustrated in Fig. 6(c) for ΔU under AS or in Appendix E.

Fig. 6 Difference between damaged and undamaged net flux over 20 cycles as defined in (52) and (53) under AD (A_d = 0.1) against: (a) stiffness damage magnitude d (fixed ω = 10); (b) loading frequency ω (fixed d = 0.35), for three values of l. Difference between damaged and undamaged net strain over 20 cycles as defined in (52) and (53) under AS (A_s = 0.2) against: (c) stiffness damage magnitude d (fixed ω = 10); (d) loading frequency ω (fixed d = 0.35), for three values of l.

The effect of stiffness damage on net flux under AD increases exponentially with decreasing frequency, as illustrated in Fig. 6(b), whilst net strain remains small for all ω (see Appendix E.1). The simulations were run for 0.5 ≤ ω ≤ 50 however the curves stop at ω = 32.5, beyond which porosity locally vanishes and the code breaks. This is a known characteristic of a poroelastic material subject to cyclic applied displacement; we refer the reader to ref. 33, 34 and 82 for further details. Also noteworthy is the fact that our simulation time is based on a fixed number of loading cycles, therefore a smaller ω increases the total loading time. The same trend of increased effect at smaller frequency is seen for ΔU under AS in Fig. 6(d), as well as for ΔQ (see Appendix E.2), although smaller in magnitude (note vertical scale).

This increase in net strain/flux for decreasing ω is due to the increased loading period and the homogenisation of the response in Z. At large frequency, the strain and flux remain localised to the moving boundary, but for smaller frequency, as the loading increases, the strain tends to be more uniform in Z, homogenising the response, and the flux tends to linear in Z. The response of the material penetrates further into the material away from the boundary, leading to greater cumulative strain and flux, which is further amplified by the increased loading period. However, ΔU under AD shows a non-monotonic dependence on ω, though the magnitude is much smaller (ΔU ∼ 10⁻², see Appendix E.1). A summary of the response of ΔU and ΔQ under different parameter changes and loading regimes is presented in Appendix E.3.

Another parameter that is natural to vary is loading amplitude, but as previously discussed loading amplitude only has a quantitative and intuitive impact, so that increasing loading amplitude A_s or A_d will only amplify the effects of heterogeneity, and vice versa.

3.2 Heterogeneous permeability

So far we have only considered heterogeneous stiffness, however heterogeneity may be introduced in other material properties, such as permeability k = k_KCf(Z). Fig. 7 displays the effect of permeability damage magnitude d on the net strain and net flux and the effect of loading frequency on net flux under AS (see Appendix F.1 for the full solution – the spatial profiles plotted at time points over a loading cycle – with d = 0.8, ω = 10). Decreasing permeability naturally reduces the flow (cf. eqn (17)) as shown by ΔQ < 0. It also decreases diffusion, localising the response to the loaded boundary as was also observed by Fiori,⁷⁰ so that net strain is decreased under AS (ΔU < 0) and net compressive strain is increased under AD (ΔU > 0, not shown). Again, |ΔU| and |ΔQ| increase with d and decrease with l. The dependence on ω however is different to stiffness damage, displaying a non-monotonic profile as illustrated in the right panel of Fig. 7. This is a result of Q both increasing in magnitude and becoming more localised to Z = 0 as ω increases, maximising the effect of permeability damage at ω = ω*. In this case, ω* ≈ 1.5 and the value of ω* varies for AS, AD, ΔU and ΔQ. A summary of the response of ΔU and ΔQ under different parameter changes and loading regimes is presented in Appendix F.2.

Fig. 7 (a) Difference between damaged and undamaged net strain over 20 cycles as defined in (52) and (53) under AS (A_s = 0.2) against permeability damage magnitude d (fixed ω = 10). (b) and (c) Difference between damaged and undamaged net flux over 20 cycles as defined in (52) and (53) under AS (A_s = 0.2) against (b) permeability damage magnitude d (fixed ω = 10) and (c) loading frequency ω (fixed d = 0.35), for three values of l.

Whilst the order of magnitude of ΔQ is comparable between permeability and stiffness damage when varying d, net flux is considerably more sensitive to loading frequency with stiffness damage (subject to their different loading strategies). Additionally, the effect of permeability damage on net strain is one to two orders of magnitude smaller than with stiffness damage. When looking at the full solution over a steady cycle however (Appendix F), the strain's qualitative profile is considerably more affected than the flux's, highlighting the care that must be taken when interpreting results from a metric such as net strain/flux, particularly if trying to infer information about the full behaviour of the material. Note that to accentuate the qualitative effect of permeability damage, Appendix F plots the full solution for the larger value d = 0.8 than was used for stiffness damage.

4 Discussion and biological relevance

In this paper, we have investigated the effect of continuous material heterogeneity on the response of a poroelastic material to a uniaxial applied stress or applied displacement. While incorporating heterogeneity and cyclic loading creates a high-dimensional parameter space, we generally found that decreased stiffness increases strain under AS and increases flux under AD, whereas decreased permeability decreases flux under AS and increases absolute strain under AD, though the effect of permeability damage is one to two orders of magnitude smaller than stiffness damage. We also examined how the net magnitude of the strain and flux response varies with damage magnitude, location and frequency of AS/AD.

As our analysis was motivated by damage in tendon, it is worth interpreting our results in this context. In soft hydrated tissues such as tendon or cartilage, cells respond to their mechanical environment. Assuming for simplicity a symmetry between tension and compression (or flux direction), the cumulative (time integrated) strain or flux provide reasonable measures for the local mechanical environment experienced by a cell. Consider a tendon subject to cyclic applied loads, with damage creating a local decrease in stiffness, and suppose for instance that the tendon cells' response to their mechanical environment correlates with cumulative strain.^83,84 In this scenario, based on our analysis, the greater cumulative strain (compared to an undamaged tendon) would indicate a modified protein production (mechanotransduction) which in turn would modify the structure and thus material properties of the tendon, feeding back into how the tendon responds to the mechanical stress.

While our model of course presents a highly idealized system that does not include much of the biological or structural complexity of a real tendon, the response features we have uncovered could in principle constitute a basic diagnostic tool indicating how severe damage is and its location, given baselines values of strain and flux for healthy tissue. Our results suggest that the flux gives a better indication on location of damage, and that slower loading will result in a greater difference in net strain/flux between the healthy and unhealthy tendon. Additionally, since the maximum strain magnitude decreases with damage location l, a measurement of maximum strain could give an indication of the damage location. In the context of tendon mechanics, the regimes of fast (ω ≥ 20), moderate (2 < ω < 20) and slow loading (ω ≤ 2) correspond to different activities: for the Achilles tendon for example, fast loading would correspond to running, moderate loading to walking or slow jogging, and slow loading to exercises such as resistance training. Note that amplitude of loading also varies with these different exercises; although we did not vary amplitude in this study, its qualitative effect in our regime of interest is negligible for the homogeneous case³³ and we can expect a similar result for heterogeneous stiffness/permeability.

In our analysis we used the net strain and flux to characterise how the response varies with damage magnitude, location and frequency. This was primarily motivated by the form of the solution which, except for strain under AS, oscillates between positive and negative values over a cycle, and thus cancels out (to some extent) when integrating over time. The idea was to capture the cumulative effect of damage on the magnitude of the response. In practice, the choice of metric should be informed by the context of the problem. For example, the mechano-sensitive response of tendon cells may in fact correlate more closely with maximum strain or flux, or stress gradients, or net pressure. The sign of these measures might also be important, responding differently in relation to tensile and compressive strain, or to inflow versus outflow. In problems concerned with solute transport, the difference between positive and negative flux is clearly key. In practice, the choice of measure is also limited by the tools available. For instance in a medical exam measuring the strain at each Z might not be possible, however measuring macroscopic strain may be more accessible. This separates research aimed at improving our understanding of the development of diseased or damaged tissue, to inform prevention and healing, and research aimed at improving diagnosis.

The features of the model were primarily explored in the steady cyclic response, that is once the response has settled down to a periodic state. Many slow activities however may be performed once instead of cyclically, so that the initial response is important. Nevertheless, the transient response during slow loading is similar (or equal) to the steady cyclic response, so we expect the analysis to remain valid. Moreover, in defining the net and cumulative metrics we integrated over time including the initial cycles, so that the transient response is captured within the metric, though not explicitly teased out.

It is also worth noting that some care should be taken when comparing results to experimental data, as we have shown that applying a stress or applying a displacement generate very different responses, particularly when heterogeneity is imposed. For example, in vivo it might be more natural to apply a stress, or even a combination of stress and displacement due to the anatomical constraints, whereas ex vivo or in vitro experiments are often performed by means of an applied displacement as this is easier to manipulate. It should therefore not come as a surprise if in vivo and ex vivo/in vitro experiments yield different results, and comparison to theory should be conducted accordingly.

5 Conclusion

We have presented a fluid-saturated poroelastic model in a Lagrangian framework for a material undergoing uniaxial cyclic tensile stress or displacement, where the solid obeys neo-Hookean elasticity and permeability is porosity-dependent according to Kozeny–Carman's law. The novel contribution from our study is the elucidation of the effect of continuous heterogeneous material properties, referred to as “damage” and imposed as a local decrease in the stiffness or permeability by a magnitude d at a location l, on the strain and flux response to an applied stress or displacement. Although damage was found to affect the strain and flux response both spatially and temporally, we were able to reduce the response to a pair of scalar metrics to explore the effect of loading frequency and damage parameters d and l. While the metrics omitted the spatial and temporal detail of the response they gave a general sense of the system's sensitivity to these parameters and can serve as a useful guide for more specific studies. A further strength of our modelling approach is, given agreement between our predictions for net quantities and experimental data, the model also determines the underlying spatial and temporal features of the fields of interest.

The results from our analysis provide an insight into how material heterogeneity modifies the response of soft porous materials to cyclic loads or displacements. While the quantitative details will depend on modelling choices – in particular the form of heterogeneity and the metrics analysed – the qualitative principles we have uncovered are generally applicable to a wide range of physical systems, from rock mechanics to hydrated tissues. Possible future work is to test this theoretical model experimentally using gels, as their synthesis allows for good control of their poroelastic properties.

We were particularly motivated by the development of disease such as tendinitis and tendinopathy, in which the material properties of the tendon such as stiffness and permeability are altered. We chose to directly damage permeability, however we could also have indirectly manipulated it by damaging the initial porosity Φ₀ instead – to fully understand the relationship between modifying porosity and the associated permeability of a poroelastic material would require a multi-scale model. Although tendons primarily bear axial load, supported by their highly aligned collagen structure,⁶⁹ the assumption of uniaxial deformation is a simplification: the literature also suggests tendons are highly anisotropic and that fluid exchange is possible in the transverse directions, so that future work could be directed towards constructing a reduced three-dimensional model which takes into account these effects. Moreover, the neo-Hookean model and the Kozeny–Carman expression for permeability were chosen for their simplicity and wide use, and substitution with more relevant constitutive laws for the material in question in this model should be straightforward, provided thermodynamic constraints are respected. For instance, a tendon specific constitutive equation for the solid skeleton could be employed instead. The question of which mechanical cues are relevant to mechanotransduction remains open and is an active area of research, and will help further inform future work on this macroscopic model of the tendon. We highlight that the simplicity of this model provides us with a clear and intuitive understanding of the behaviour of heterogeneous poroelastic materials, which can be used to inform and build more complex, higher-dimensional models and incorporate additional biology.

Author contributions

Designed research: all authors; conducted research: Z. C. G.; wrote paper: all authors.

Conflicts of interest

There are no conflicts to declare.

Data availability

The code required to reproduce model output is available in the public repository: https://github.com/zoegdr/hetpore1D.

Appendices

A Surface transport equations

We derive two useful results which relate the flow of Eulerian quantities through a material surface da to the flow of Lagrangian quantities through a material surface dA. Nanson's formula states that a material surface dA with unit normal N transforms to the material surface da with unit normal n as

nda = JF^−TNdA or n_ida = JF_ji⁻¹N_jdA. (54)

Consider an arbitrary vector M in the reference configuration which transforms to the vector m attached to the current configuration, such that the flow of M through dA is equivalent to the flow of m through da:

m·nda = M·NdA. (55)

Using Nanson's formula this can be rearranged asMoreover, integrating (55) and applying the divergence theorem for arbitrary volume dv and dV₀:

∇_X·MdV₀ = ∇_x·mdv. (57)

Eqn (56) is useful when converting relevant Eulerian quantities to Lagrangian ones, and eqn (57) will be useful when converting from Eulerian gradient to Lagrangian gradient, provided eqn (55) holds.

B Numerical method

This numerical method is implemented to solve the advection-diffusion equation for porosity. Subsequent fields, such as pressure, displacement, stress and strain, are calculated after. Consider the finite volume mesh where the spatial domain [0, 1] is divided into N grid cells, where the cell walls are located at the black mesh points 0 ≤ i ≤ N, and the cell centres are located at the red mesh points , as depicted in Fig. 8. Let Z_i = iΔZ where ΔZ = 1/N so that Z₀ = 0 and Z_N = 1, and denote the unknowns [capital Phi, Greek, circumflex](Z_i,t) = [capital Phi, Greek, circumflex]_i(t). We solve for [capital Phi, Greek, circumflex]_i+1/2(t) at the cell centres, and boundary conditions are applied at i = 0 and i = N.

Fig. 8 Finite volume mesh for a one-dimensional poroelastic material, with stress or displacement applied at Z = 0 and no flux at Z = 1. Unknown values [capital Phi, Greek, macron]_i are solved at each of the red nodes.

Rewriting (41) in the following form,

where G([capital Phi, Greek, circumflex]) = k([capital Phi, Greek, circumflex])/(1 + [capital Phi, Greek, circumflex]) we integrate the above equation over cell i and divide by ΔZ:

Let the local cell average be given by

and the corresponding stress s̄_i+1/2 = s′([capital Phi, Greek, macron]_i+1/2), calculated using eqn (27) with J_i = 1 + [capital Phi, Greek, macron]_i. Let the numerical finite volume flux at the cell walls F_i(t) be given bywhere we have used a central average to approximate D(Φ_i) and a centered finite difference to approximate the derivatives for inner mesh points 1 ≤ i ≤ N − 1. The PDE (58) is now a system of ODEs for each timestep,which we can solve using ode15s in MATLAB with the appropriate boundary and initial conditions. The above formulation can also be used for the form (14) of the non-linear flow equation, giving

F_i(t) = −Q(Z_i,t) = V^s(Z_i,t) (63)

which is useful when implementing boundary conditions in the following section.

B.1 Implementation of boundary and initial conditions.
B.1.1 Zero flux and solid displacement. The zero flux boundary condition at i = N (31) can be rewritten as Q(Z_N,t) = 0, so that from (63)

F_N(t) = 0. (64)

B.1.2 Applied stress. For an applied stress s*(t) at Z = 0, we approximate the derivative ∂s′/∂Z using a forward finite difference such thatwhere the corresponding value Φ* := [capital Phi, Greek, circumflex](0,t) associated with s* is calculated by inverting (27) and choosing the positive root (this is the only physical choice for ν < 1/2). The finite volume flux at Z = 0 is then:
B.1.3 Applied displacement. For an applied displacement at Z = 0, we make use of the form (63) so that

F₀(t) = ȧ. (67)

B.1.4 Initial condition. The ODE time solver is initialised with [capital Phi, Greek, macron]_i(0) = 0.

B.2 Other quantities. Displacement is calculated from porosity by integrating eqn (9):Pressure is calculated by integrating eqn (24):

P(Z,t) = s′(Z,t) − s*(t) (69)

or in the case of applied displacement:which is approximated since s′(Z₀,t) cannot be computed, and stress is calculated from porosity using (43). Strain is simply given by

U_Z = [capital Phi, Greek, circumflex] (71)

and the relative flowwhich quantifies how much fluid is being transported through the material at a particular point.

B.3 Numerical convergence. To determine an appropriate step size, we solve the system under an applied stress and varying the grid step-size dZ = 1/N from 25 to 1000. We calculate the absolute relative error asevaluating displacement at the applied stress boundary (Z = 1/2N) as this involves integrating Φ over the whole domain. For applied displacement we perform the same error calculation but evaluating stress at Z = 1/2N.

The relative error, displayed in (Fig. 9), decreases with N and we compute our results with N = 400(log(N) ≈ 6) so as not to compromise on computation speed. We note that the method to enforce the stress boundary condition converges much faster than the displacement boundary condition.

Fig. 9 Relative error calculated for displacement decreasing the step-size (increasing N) for applied stress (A_s = 0.2, ω = 10) and applied displacement (A_s = 0.2, ω = 10).

C Transient phase decay and Péclet number

In Fig. 10 we plot the strain response under AS against time at Z = 0, 0.25, 0.5, 0.75, 1 for 5 different values of ω. The first four are plotted for 20 cycles, over which the transient phase response clearly decays and the response settles to steady periodic behaviour. For ω = 100 we plot the response over 100 cycles to emphasise the fact that the response indeed settles to periodic behaviour which persists over a large number of cycles, as indicated by the clear envelope in this plot. We note that the time it takes for the transient response to decay is the same for all ω, as was also shown by Fiori et al.,³³ meaning the number of cycles within the transient phase increases with ω as visible in these plots. Though this does not constitute a formal proof that the transient behaviour fully decays over large times, which is beyond the scope of this paper, it is suggestive that the system settles to a steady periodic state after a finite time.

Fig. 10 Strain under AS (A_s = 0.2) at Z = 0, 0.25, 0.5, 0.75, 1 (increasing light to dark blue) against t for increasing ω, over 20 cycles for the first four and 100 cycles in the final plot.

Fig. 10 also highlights the increased localisation of the deformation to Z = 0 with increasing ω, which may also be characterised by the Péclet number Pe, defined as the ratio of diffusion time to advection time. We identify the advection timescale with the loading period, and the diffusion time is simply 1 since we have non-dimensionalised by the poroelastic diffusion timescale T_pe, giving Pe = 1/T = ω/2π. For Pe ≪ 1 (ω ≪ 2π), the advection timescale is much larger than the diffusion timescale so that the fluid response is dominated by diffusion and thus primarily driven by pressure gradients; the deformation tends to homogeneous in Z. For Pe ≫ 1 (ω ≫ 2π), the opposite is true so that the fluid response is dominated by advection and thus primarily driven by solid deformation; the material behaves more like an elastic solid, and the deformation is localised to the loaded boundary. The transition from diffusion dominated to advection dominated flow has been studied theoretically and experimentally on gels within the framework of poroelasticity to investigate solute transport, with applications to drug delivery.^85–87

D Increased stiffness

Consider a local increase in stiffness instead of a decrease (to do so we simply change the sign of the Gaussian in eqn (40)). This means that the material is locally stiffer, thus requiring more stress to be pulled and compressed to the same displacement. Under AS we can expect smaller strain values at the point of damage due to the stiffer region. Under AD, we expect a flux disturbance which mirrors the one seen for the local decrease in stiffness.

In Fig. 11 we plot the time integrated absolute strain and flux under AS (first row) and under AD (second row) for a locally increased stiffness of magnitude d = 0.35 and locations l = 0.25, 0.5, 0.75 (to do so we simply change the sign of the Gaussian in eqn (40)). This mirrors the effect seen for a local decrease in stiffness, whereby strain is locally decreased under AS, and flux inversely perturbed under AD about Z = l. We note the magnitude of the response is smaller compared to decreased stiffness. Computing the areas under these curves, i.e. the net absolute strain and flux, we find the net strain to be approximately equal for the three locations under AS, though this time it is slightly increased for increasing l (whilst for decreased stiffness it was slightly increased for decreasing l). Under AD, the effect of location is not linear as before: is greatest for l = 0.75 but smallest for l = 0.5.

Fig. 11 Time integrated absolute strain and flux profiles against Z under AS (A_s = 0.2, ω = 10, left plots) and AD (A_d = 0.1, ω = 10, right plots), for homogeneous stiffness (dotted grey line) and for local increase in stiffness with d = 0.35 and varying location l = 0.25, 0.5, 0.75.

E Additional parameter variation for stiffness damage

E.1 Applied displacement. Fig. 12 displays the effect of damage magnitude d and loading frequency ω on the net strain under AD. Clearly, the effect of damage on net strain under AD is very small, with ΔU reaching a maximum magnitude of O(10⁻²). The non-monotonicity in ω observed in (b) is a result of the competition between localisation of the response and increase in compressive strain due to increasing ω.

Fig. 12 Difference between damaged and undamaged net strain under AD (A_d = 0.1) against (a) d and (b) ω for three locations l.

E.2 Applied stress. Fig. 13 displays the effect of damage magnitude d and loading frequency ω on net flux under AS. Net flux increases with d for l = 0.25 and l = 0.5 similarly to under AD. For l = 0.75 however net flux decreases with d. Net flux also decreases with ω, though comparatively less than under AD. In addition for ω < 3.5, ΔQ increases with l and for ω > 3.5 the relationship is reversed so that ΔQ decreases with l. Decreasing frequency homogenises the response so that maximum strain at Z = l is the same for all l; the same increase in strain but farther away from the boundary then requires a greater flux, so that ΔQ increases with l for small ω.

Fig. 13 Difference between between damaged and undamaged net flux against (a) d and (b) ω under AS (A_s = 0.2) for three locations l.

E.3 Summary table. We present a table summarising results from our parameter variation, including the sign of ΔU and ΔQ, their relationship with d, l and ω and their order of magnitude. Upward (downward) arrows are shorthand for increase (decrease), and refer to the absolute value of ΔU, ΔQ.

	Applied stress	Applied displacement
ΔU	>0	≈0
	↗ with d
	≈for all l
	↘ with ω
	O(10^−1 − 1)
ΔQ	>0	>0
	↗ with d	↗ with d
	↘ with ω	↘ with ω
	O(10^−1)	O(10^−1 − 10)

F Permeability damage

F.1 Full solution. Fig. 14 displays the full strain and flux solutions over one cycle for 9 equally spaced values of 38π/ω ≤ t_i < 40π/ω, under an applied stress (AS) with (A_s, ω) = (0.2, 10) in the first two rows and under an applied displacement (AD) with (A_d, ω) = (0.1, 10) in the last two rows. We differentiate between the loading phase, ṡ* > 0 or ȧ < 0 (dark blue to light blue) and unloading phase, ṡ* < 0 or ȧ > 0 (light red to dark red). The first column is the homogeneous case for reference and columns 2–4 display profiles with imposed local decrease in permeability with d = 0.8 where the location of the dip is increased from left to right (moving away from the moving boundary).

Fig. 14 Strain (first and third rows) and flux (second and fourth rows) plotted for 9 equally spaced values of cycle 38π/ω ≤ t < 40π/ω, for uniform stiffness (first column) and locally decreased permeability with d = 0.8 (columns 2–4). The first two rows are in response to an applied stress (A_s = 0.2, ω = 10) and the second two rows are in response to an applied displacement (A_d = 0.1, ω = 10). The applied stress s*(t) and displacement A(t), and their respective time derivatives ṡ*(t), ȧ(t) are displayed as insets in the first column. We differentiate between the loading phase (ṡ* > 0, ȧ < 0, dark blue to light blue) and unloading phase (ṡ* < 0, ȧ > 0, light red to dark red).

On the one hand under AS there is a clear decrease in strain for Z > l (first row), since diffusion is decreased. This is also accompanied with a decrease in flux (second row). Under AD on the other hand, strain is increased for Z < l (third row) and flux decays more rapidly with Z for Z > l (fourth row). We note that these results are for d = 0.8 whilst the profiles for d = 0.35 show very little change compared to the homogeneous case.

F.2 Summary table of parameter variation. We present a table summarising results from our parameter variation, including the sign of ΔU and ΔQ, their relationship with d, l and ω and their order of magnitude. Upward (downward) arrows are shorthand for increase (decrease), and refer to the absolute value of ΔU, ΔQ.

	Applied stress	Applied displacement
ΔU		>0
	↗ with d	↗ with d
	↘ with l	↘ with l
	O(10^−2 − 10^−1)	O(10^−2 − 10^−1)
ΔQ	<0	<0
	↗ with d	↗ with d
	↘ with l	↘ with l
	O(10^−1)	O(10^−1)

Acknowledgements

Z. C. G. gratefully acknowledges funding from the UKRI and the Paul Shreder scholarship. For the purpose of open access, the author has applied a CC BY public copyright licence to any author accepted manuscript arising from this submission. The authors thank Chris MacMinn and Matilde Fiori for helpful discussions.

References

1. K. Von Terzaghi, Sitzungsber. Akad. Wiss. Wien, Math.-Naturwiss. Kl., Abt. 2A, 1923, 132, 105–124.
2. L. Rendulic, Der Bauingenieur, 1936, 17, 559–564.
3. M. A. Biot, Ann. Soc. Sci. Bruxelles, Ser. B, 1935, 55, 110–113.
4. M. A. Biot, J. Appl. Phys., 1941, 12, 155–164.
5. M. A. Biot and G. Temple, Indiana Univ. Math. J., 1972, 21, 597–620.
6. E. Detournay and A. H.-D. Cheng, Analysis and Design Methods, Pergamon, Oxford, 1993, pp. 113–171.
7. B. Jin, Arch. Appl. Mech., 2004, 74, 277–287.
8. Y.-J. Hu, Y.-Y. Zhu and C.-J. Cheng, Soil Dyn. Earthquake Eng., 2011, 31, 491–501.
9. F. Genna and A. Cividini, Soil Dyn. Earthquake Eng., 1989, 8, 189–201.
10. R. Popescu, J. H. Prevost, G. Deodatis and P. Chakrabortty, Soil Dyn. Earthquake Eng., 2006, 26, 648–665.
11. A. Bonazzi, B. Jha and F. P. de Barros, Int. J. Numer. Anal. Methods Geomech., 2021, 45, 307–324.
12. T. Yamamoto, H. Koning, H. Sellmeijer and E. Van Hijum, J. Fluid Mech., 1978, 87, 193–206.
13. O. Madsen, Geotechnique, 1978, 28, 377–393.
14. M. G. Trefry, D. R. Lester, G. Metcalfe and J. Wu, Water Resour. Res., 2019, 55, 3347–3374.
15. M. Karim, T. Nogami and J. Wang, Int. J. Solids Struct., 2002, 39, 6011–6033.
16. N. L. Millar, K. G. Silbernagel, K. Thorborg, P. D. Kirwan, L. M. Galatz, G. D. Abrams, G. A. Murrell, I. B. McInnes and S. A. Rodeo, Nat. Rev. Dis. Primers, 2021, 7, 1.
17. M. Lavagnino, M. E. Wall, D. Little, A. J. Banes, F. Guilak and S. P. Arnoczky, J. Orthop. Res., 2015, 33, 813–822.
18. J. Bojsen-Møller and S. P. Magnusson, J. Appl. Physiol., 2019, 126, 1800–1807.
19. A. Arampatzis, K. Karamanidis and K. Albracht, J. Exp. Biol., 2007, 210, 2743–2753.
20. S. Docking, T. Samiric, E. Scase, C. Purdam and J. Cook, Muscles Ligaments Tendons, 2013, 3, 7.
21. F. Eckstein, M. Tieschky, S. Faber, K.-H. Englmeier and M. Reiser, Anat. Embryol., 1999, 200, 419–424.
22. M. Lavagnino, S. P. Arnoczky, E. Kepich, O. Caballero and R. C. Haut, Biomech. Model. Mechanobiol., 2008, 7, 405–416.
23. K. Piekarski and M. Munro, Nature, 1977, 269, 80–82.
24. D. Zhang and S. C. Cowin, J. Mech. Phys. Solids, 1994, 42, 1575–1599.
25. P. Manfredini, G. Cocchetti, G. Maier, A. Redaelli and F. M. Montevecchi, J. Biomech., 1999, 32, 135–144.
26. V.-H. Nguyen, T. Lemaire and S. Naili, Med. Eng. Phys., 2010, 32, 384–390.
27. F. Witt, G. N. Duda, C. Bergmann and A. Petersen, Tissue Eng., Part A, 2014, 20, 486–493.
28. L. Zhang, Int. J. Appl. Mech. Eng., 2011, 3, 507–524.
29. P. Riches, N. Dhillon, J. Lotz, A. Woods and D. McNally, J. Biomech., 2002, 35, 1263–1271.
30. R. L. Mauck, C. T. Hung and G. A. Ateshian, J. Biomech. Eng., 2003, 125, 602–614.
31. B. G. Sengers, C. W. Oomens and F. P. Baaijens, J. Biomech. Eng., 2004, 126, 82–91.
32. J. Cacheux, J. Ordonez-Miranda, A. Bancaud, L. Jalabert, D. Alcaide, M. Nomura and Y. T. Matsunaga, Sci. Adv., 2023, 9, eadf9775.
33. M. Fiori, S. Pramanik and C. W. MacMinn, J. Fluid Mech., 2023, 974, A2.
34. M. Fiori, S. Pramanik and C. W. MacMinn, J. Fluid Mech., 2025, 1009, A15.
35. G. Franceschini, D. Bigoni, P. Regitnig and G. Holzapfel, J. Mech. Phys. Solids, 2006, 54, 2592–2620.
36. R. T. Kedarasetti, P. J. Drew and F. Costanzo, Fluids Barriers CNS, 2022, 19, 34.
37. L. Bojarskaite, A. Vallet, D. M. Bjørnstad, K. M. Gullestad Binder, C. Cunen, K. Heuser, M. Kuchta, K.-A. Mardal and R. Enger, Nat. Commun., 2023, 14, 953.
38. E. Moeendarbary, L. Valon, M. Fritzsche, A. R. Harris, D. A. Moulding, A. J. Thrasher, E. Stride, L. Mahadevan and G. T. Charras, Nat. Mater., 2013, 12, 253–261.
39. Y. Hu, X. Chen, G. M. Whitesides, J. J. Vlassak and Z. Suo, J. Mater. Res., 2011, 26, 785–795.
40. J. Delavoipiere, Y. Tran, E. Verneuil and A. Chateauminois, Soft Matter, 2016, 12, 8049–8058.
41. C.-Y. Hui, Y. Y. Lin, F.-C. Chuang, K. R. Shull and W.-C. Lin, J. Polym. Sci., Part B: Polym. Phys., 2006, 44, 359–370.
42. Y. Hu, X. Zhao, J. J. Vlassak and Z. Suo, Appl. Phys. Lett., 2010, 96, 121904.
43. G. D. Degen, Y.-T. Chen, A. L. Chau, L. K. Månsson and A. A. Pitenis, Soft Matter, 2020, 16, 8096–8100.
44. C. Kopecz-Muller, V. Bertin, E. Raphaël, J. D. McGraw and T. Salez, Proc. R. Soc. A, 2023, 479, 20220832.
45. Y. Kameo, Y. Ootao and M. Ishihara, Biomech. Model. Mechanobiol., 2016, 15, 361–370.
46. C. Ruiz, J. Noailly and D. Lacroix, J. Mech. Behav. Biomed. Mater., 2013, 26, 1–10.
47. Q. Wang, J. Herrmann, K. Worthington and E. Sander, Sci. Rep., 2025, 15, 5997.
48. B. Le Roi and J. M. Grolman, ACS Macro Lett., 2025, 14, 176–181.
49. Z. Xu, H. Chen, H.-B. Yang, X. Yao, H. Qin, H.-P. Cong and S.-H. Yu, Nat. Commun., 2025, 16, 400.
50. C. S. Wyss, P. Karami, P.-E. Bourban and D. P. Pioletti, Extreme Mech. Lett., 2018, 24, 66–74.
51. D. Zhong, Z. Wang, J. Xu, J. Liu, R. Xiao, S. Qu and W. Yang, Nat. Commun., 2024, 15, 5896.
52. Y. Han, Y. Wu, F. Wang, G. Li, J. Wang, X. Wu, A. Deng, X. Ren, X. Wang and J. Gao, et al., Bioact. Mater., 2024, 35, 1–16.
53. S. Arya and K. Kulig, J. Appl. Physiol., 2010, 108, 670–675.
54. S. J. Obst, L. J. Heales, B. L. Schrader, S. A. Davis, K. A. Dodd, C. J. Holzberger, L. B. Beavis and R. S. Barrett, Sports Med., 2018, 48, 2179–2198.
55. V.-H. Nguyen, T. Lemaire and S. Naili, Biomech. Model. Mechanobiol., 2011, 10, 963–972.
56. G. A. Orozco, P. Tanska, A. Gustafsson, R. K. Korhonen and H. Isaksson, J. Mech. Behav. Biomed. Mater., 2022, 131, 105227.
57. I. Ambartsumyan, V. J. Ervin, T. Nguyen and I. Yotov, ESAIM:Math. Modell. Numer. Anal., 2019, 53, 1915–1955.
58. M. Bajuri, H. Isaksson, P. Eliasson and M. S. Thompson, Biomech. Model. Mechanobiol., 2016, 15, 1457–1466.
59. T. Shearer, J. Biomech., 2015, 48, 290–297.
60. T. Shearer, J. Biomech., 2015, 48, 3017–3025.
61. H. Gupta, J. Seto, S. Krauss, P. Boesecke and H. Screen, J. Struct. Biol., 2010, 169, 183–191.
62. T. A. Wren, G. S. Beaupre and D. R. Carter, J. Rehabil. Res. Dev., 2000, 37, 135–144.
63. B. N. Safa, E. T. Bloom, A. H. Lee, M. H. Santare and D. M. Elliott, J. Biomech., 2020, 109, 109892.
64. T. Atkinson, R. Haut and N. Altiero, J. Biomech. Eng., 1997, 400–405.
65. S. Butler, S. Kohles, R. Thielke, C. Chen and R. Vanderby, Med. Biol. Eng. Comput., 1997, 35, 742–746.
66. C.-T. Chen, D. S. Malkus and R. Vanderby Jr, Biorheology, 1998, 35, 103–118.
67. I. Scott, PhD thesis, University of Oxford, 2022.
68. H. Khayyeri, A. Gustafsson, A. Heuijerjans, M. K. Matikainen, P. Julkunen, P. Eliasson, P. Aspenberg and H. Isaksson, PLoS One, 2015, 10, e0126869.
69. C. T. Thorpe, H. L. Birch, P. D. Clegg and H. R. Screen, Tendon Regeneration, Elsevier, 2015, pp. 3–39.
70. M. Fiori, PhD thesis, University of Oxford, 2023.
71. C. W. MacMinn, E. R. Dufresne and J. S. Wettlaufer, Phys. Rev. Appl., 2016, 5, 2331–7019.
72. O. Coussy, Poromechanics, John Wiley & Sons, 2004.
73. A. Kent, S. L. Waters, J. Oliver and S. Chapman, SIAM J. Appl. Math., 2023, 83, 2118–2143.
74. P. Zheng and K. Zhang, Int. J. Mech. Sci., 2019, 161, 105074.
75. R. Gatt, M. V. Wood, A. Gatt, F. Zarb, C. Formosa, K. M. Azzopardi, A. Casha, T. P. Agius, P. Schembri-Wismayer and L. Attard, et al., Acta Biomater., 2015, 24, 201–208.
76. S. P. Reese, S. A. Maas and J. A. Weiss, J. Biomech., 2010, 43, 1394–1400.
77. C. T. Thorpe, G. P. Riley, H. L. Birch, P. D. Clegg and H. R. Screen, Acta Biomater., 2014, 10, 3217–3224.
78. J. H. Yoo, S. R. Yi and J. H. Kim, Surg. Radiol. Anat., 2007, 29, 623–628.
79. T. Kurihara, R. Sasaki and T. Isaka, ISBS-Conference Proceedings Archive, 2012.
80. R. Beyer, M. Kongsgaard, B. Hougs Kjær, T. Øhlenschlæger, M. Kjær and S. P. Magnusson, Am. J. Sports Med., 2015, 43, 1704–1711.
81. A. Rowlands, M. R. Stone and R. G. Eston, Med. Sci. Sports Exercise, 2007, 39, 716.
82. D. R. Hewitt, D. T. Paterson, N. J. Balmforth and D. M. Martinez, Phys. Fluids, 2016, 28, 24.
83. R. Nakamichi and H. Asahara, J. Bone Miner. Res., 2024, 39, 814–820.
84. J. Humphrey, J. Biomech. Eng., 2001, 123, 638–641.
85. B. Gardiner, D. Smith, P. Pivonka, A. Grodzinsky, E. Frank and L. Zhang, Comput. Methods Biomech. Biomed. Eng., 2007, 10, 265–278.
86. F. Urciuolo, G. Imparato and P. A. Netti, AIChE J., 2008, 54, 824–834.
87. B. L. Vaughan, P. A. Galie, J. P. Stegemann and J. B. Grotberg, Biophys. J., 2013, 105, 2188–2198.

---