From arteries to boreholes: steady-state response of a poroelastic cylinder to fluid injection

The radially outward flow of fluid into a porous medium occurs in many practical problems, from transport across vascular walls to the pressurization of boreholes. As the driving pressure becomes non-negligible relative to the stiffness of the solid structure, the poromechanical coupling between the fluid and the solid has an increasingly strong impact on the flow. For very large pressures or very soft materials, as is the case for hydraulic fracturing and arterial flows, this coupling can lead to large deformations and, hence, to strong deviations from a classical, linear-poroelastic response. Here, we study this problem by analysing the steady-state response of a poroelastic cylinder to fluid injection. We consider the qualitative and quantitative impacts of kinematic and constitutive nonlinearity, highlighting the strong impact of deformation-dependent permeability. We show that the wall thickness (thick versus thin) and the outer boundary condition (free versus constrained) play a central role in controlling the mechanics.

For a simple uniaxial deformation, Hencky elasticity reduces to where σ is the normal effective stress, M is the oedometric modulus, and λ = 1 + ∆L/L is the stretch, with ∆L the change in overall length and L the original length. Linear elasticity instead predicts We compare these behaviours in Figure A1.

B. Injection via fixed pressure drop
To enforce a constant pressure drop ∆p, we must derive an associated expression for the evolving flow rate q(t). To do so, we rearrange and integrate the expression for vs from equation (2.12b) to 10 -2 10 -1 10 0 10 1 10 2 6 0 0.5 1 1.5 2 j< 0 =Mj Figure A1. The absolute value of the dimensionless effective stress vs. the stretch for uniaxial deformation according to linear elasticity (red) and Hencky elasticity (blue). Hencky elasticity provides a stiffer response than linear elasticity in compression (λ < 1) and a softer response in tension (λ > 1). In tension, the stress predicted by Hencky elasticity reaches a maximum value of σ /M = 1/e for a stretch of λ = e (blue diamond) before decreasing asymptotically to zero. The two models agree in the limit of small strain, For linearised kinematics, we replace the latter expression with vs ≈ ∂us ∂t .
The above expressions also apply at steady state, for which vs ≡ 0.
C. Linear elasticity with constant permeability (L-k 0 and Q-k 0 ) Assuming linear elasticity, we solve equation (3.2) for constant permeability (k[φ f (us)] ≡ 1) to arrive at a general expression for the displacement, where B 1 and B 2 are determined by the boundary conditions. This result is solely mechanical and constitutive, and is therefore valid for both the L-k 0 and Q-k 0 models. The general expressions for the effective stresses are then From these expressions, we arrive at four distinct solutions by combining the two different treatments of the kinematics (rigorous Q and linearised L) with the two different sets of outer boundary conditions (an applied stress at the outer boundary (equation 2.27) and a fixed outer boundary (equation 2.28)). The two L-k 0 solutions are classical solutions from linear poroelasticity [16]. An approximate version of the Q-k 0 solutions was derived by Barry and Aldis [10] and Barry and Mercer [11], who applied boundary conditions at the moving boundary but linearised the relationship between φ f and us.
(a) Solution for L-k 0 with an applied effective stress at the outer boundary For an applied effective stress at the outer boundary, we derive expressions for B 1 and B 2 by applying the appropriate inner and outer boundary conditions (equations (2.25) and (2.27), respectively). We linearise the kinematics by applying these at r = a 0 (rather than at a) and at r = 1 (rather than at b), respectively. We obtain

fixed outer boundary
Similarly, for a fixed outer boundary, we apply the appropriate inner and outer boundary conditions (equations (2.25) and (2.28), respectively) at r = a 0 and at r = 1, respectively, to obtain and All other quantities can be derived from the expressions for us. Thus, we have complete explicit solutions following classical linear poroelasticity for the two different sets of outer boundary conditions. Note that, for linearised kinematics, φ f should be calculated from us according to equation (2.16).
(c) Solution for Q-k 0 with an applied effective stress at the outer boundary For an applied effective stress at the outer boundary, we now apply equations (2.25) and (2.28) at r = a and r = b, respectively, to the general elastic solution (equation C1). This leads to This solution is not explicit because the inner radius a and outer radius b are now determined by the two kinematic conditions (see equations (2.25) and (2.27)), leading to two coupled, implicit expressions for a and b. We solve these expressions numerically using a root-finding technique.

(d) Solution for Q-k 0 with a fixed outer boundary
For a fixed outer boundary, we now obtain and The problem is closed by applying the kinematic condition at the inner boundary (see equation (2.25)), leading to an implicit expression for a and the fixed outer boundary condition, equation (2.28). We again solve this numerically using a root-finding technique. As above, all other quantities can then be derived from the expressions for us. Note that, for rigorous kinematics, φ f should be calculated from us according to equation (2.7).

D. Numerical solution via Chebyshev spectral collocation
When the ODE presented in §3 cannot be solved analytically, it must instead be integrated using standard numerical methods for BVPs, such as direct finite differences or a shooting method. For a shooting method, one must guess the locations of the free boundaries, solve the ODE as an initial value problem (IVP) subject to two of the constraints, and then iterate on the guesses until the remaining constraints are satisfied. For direct finite differences, two approaches are possible. One may follow the same approach as for a shooting method, but solving the BVP directly using finite differences and root finding (e.g., Newton's method) rather than solving it as an IVP. Alternatively, one may solve the BVP and all constraints simultaneously using finite differences and root finding. Although straightforward to implement, these approaches are unreliable in the present context because the iteration process can easily lead to a nonphysical state that prohibits further iteration. To mitigate these difficulties, we instead use a direct method based on Chebyshev spectral collocation (i.e., a Chebyshev pseudospectral method) [e.g., 35, 36]. That is, we solve the BVP and all constraints simultaneously as described above, but replacing the sparse finite-difference differentiation matrix with a dense Chebyshev-pseudospectral differentiation matrix. This approach still requires Newton iteration, but is more robust than finite differences because the density of the pseudospectral differentiation matrix directly couples the solution at each discrete point to the solution at every discrete point. This approach also allows for the straightforward incorporation of additional unknowns and constraints, such as solving the problem for an imposed pressure drop ∆p rather than for an imposed flow rate q. We illustrate Repeat with new u.
Evaluate the discrete nonlinear differential operator F(u) (the residual).
Calculate the Jacobian matrix ∂F/∂u.
Update the solution using Newton's method.
Solution found.
NO YES Figure D1. Procedure for direct solution via Chebyshev spectral collocation method.
the overall structure of the method in Figure D1. Note that, for purposes of Newton iteration, we calculate the Jacobian analytically for the L and Q models and numerically for the N models. Spectral collocation methods involve discretising the solution domain into a set of N points (collocation points), defining a global function that interpolates the solution at these collocation points (the interpolant), and then approximating the derivatives of the solution as the derivatives of the interpolant. In a Chebyshev spectral collocation method, the collocation points are the N Chebyshev points x k ∈ [1, −1], which can be defined as [35] x k = cos (k − 1)π N − 1 , k = 1, . . . , N.
The basis functions from which the interpolant is composed are then a set of N polynomials of degree N − 1 satisfying the criterion that each is nonzero at exactly one distinct collocation point. Note that other definitions of the Chebyshev points are also commonly used [e.g., 34]. For the definition given in equation (D1), Weideman and Reddy [37] provide a suite of MATLAB functions that generate the Chebyshev points and differentiation matrices, and that perform interpolation.

E. Impact of geometry
In Figures 4 and 5 of the main text, we plot the evolution of various key quantities as contours of fixed ∆p against a 0 . It is useful for interpretation to present the same results in several different ways. Here, we show the results of Figure 4 as contours of fixed q against a 0 (figure E1), contours of fixed a 0 against q (figure E2), and contours of fixed a 0 against ∆p (figure E3). We also do the same for Figure 5 (figure E4).

F. Rheological effects
We focused in the main text on results from the Q-k KC model because it provides a good compromise between rigour, robustness, and computational efficiency. We examine this choice in more detail in Figure F1 by plotting q against a 0 for contours of fixed ∆p, as in the last row of Figure 4, for five different models: Q-k KC (first row), L-k KC (second row), Q-k 0 (third row), N-k KC (last row), and L-k 0 (all rows, grey lines). Note that q appears to be much more sensitive to the permeability law than other aspects of the deformation (cf. Figures 2 and 3), making it a useful quantity for this comparison. For unconstrained cylinders (left column), the Q-k KC and N-k KC models predict qualitatively similar behaviour, with the contours in the latter bending to the left somewhat more strongly. The latter model is also much more computationally expensive. The Q-k 0 model exhibits similar behaviour, but with much more extreme bending of the contours (note the different vertical scale). Our results for the L-k KC model are inconclusive because this model is much less robust than either of the Q models; our method fails to find a solution for even moderate values of q. This is likely because the L-k KC model is asymptotically inconsistent and does not correctly capture the kinematic relationship between porosity and displacement. The Q-k KC is much more rigorous in these regards, and is only slightly more computationally expensive in our pseudospectral collocation framework.
For constrained cylinders (right column), all three of the k KC models exhibit very similar behaviour despite the different elasticity laws (L and Q vs. N) and the different treatments of the kinematics (L vs. Q and N). This presentation does not constitutive a careful quantitative comparison, but it suggests that deformation-dependent permeability plays a key role in the mechanics of the problem, particularly for high pressures (right), whereas large-deformation kinematics are less important. This is somewhat unsurprising since constrained cylinders generally deform much less than unconstrained cylinders.
These results suggest that rigorous large-deformation kinematics (including the relationship between porosity and displacement) are important for model robustness and are central to the double-valued behaviour of unconstrained cylinders. Deformation-dependent permeability appears to moderate (but not eliminate) the double-valued behaviour of unconstrained cylinders, and to be central to the behaviour of constrained cylinders.
As noted in the main text, we have chosen Kozeny-Carman permeability and Hencky elasticity as relatively generic constitutive laws that capture the qualitatively important features of deformation-dependent permeability and large-deformation elasticity, respectively. Given the strong role of deformation-dependent permeability in our results, a comparison with results for other permeability laws would be an interesting topic for future work.

G. Solution for Q-k 0 in the thin-walled limit
We now derive an approximate solution to the Q-k 0 model in the limit of vanishing wall thickness, starting from equation (3.2) with k[φ f (us)] ≡ 1. We do this for the case of an applied effective stress at the outer boundary since the case of no displacement at the outer boundary is ultimately limited to small displacements and thus is well-captured by linear poroelasticity.
We begin by defining a new radial coordinate ≡ r − a such that ∈ [0, δ] where δ ≡ b − a 1 is the wall thickness. We then rewrite equation (3.2) in terms of and seek a solution under the assumption that 1. From these assumptions, and writing us(r) = U ( ), we obtain at leading order Note that these assumptions require for asymptotic consistency that q/a = O(1). Equation (G1) is a linear, second-order ODE with solution We now apply the relevant boundary conditions (equations (2.25) and (2.27) with σ r ≡ 0), which results in four equations for four unknowns: The two integration constant, A 1 and A 2 , and the inner and outer radii, a and a + δ, respectively. We use the two conditions at the inner boundary to derive expressions for A 1 and A 2 in terms of a, The two conditions at the outer boundary then give and (G4b) This is now a root-finding problem for the values of a and δ, which we solve numerically using the standard MATLAB function fsolve.
The pressure field is given by This then leads to ∆p = (q/a)δ and, since q/a = O(1), we have that ∆p = O(δ). This implies that a small pressure drop will drive a large flow rate when the walls are sufficiently thin.