A stabilized finite element method for finite-strain three-field poroelasticity

We construct a stabilized finite-element method to compute flow and finite-strain deformations in an incompressible poroelastic medium. We employ a three-field mixed formulation to calculate displacement, fluid flux and pressure directly and introduce a Lagrange multiplier to enforce flux boundary conditions. We use a low order approximation, namely, continuous piecewise-linear approximation for the displacements and fluid flux, and piecewise-constant approximation for the pressure. This results in a simple matrix structure with low bandwidth. The method is stable in both the limiting cases of small and large permeability. Moreover, the discontinuous pressure space enables efficient approximation of steep gradients such as those occurring due to rapidly changing material coefficients or boundary conditions, both of which are commonly seen in physical and biological applications.


Introduction
Poroelasticity theory assumes a superposition of solid and fluid components to capture complex interactions between a deformable porous medium and the fluid flow within it, and was originally developed to study geophysical applications such as reservoir geomechanics [26,28,41]. Fully saturated, incompressible poroelastic models have since been used to model a variety of biological tissues and processes. Biological examples include the coupling of flow in coronary vessels with the mechanical deformation of myocardial tissue to create a poroelastic model of coronary perfusion [13,15]. Other examples include modelling tissue deformation and the ventilation in the lungs [6], protein-based hydrogels embedded within cells [22], brain oedema and hydrocephalus [39,56], microcirculation of blood and interstitial fluid in the liver lobule [36], and interstitial fluid and tissue in articular cartilage and intervertebral discs [21,24,42].
When using the finite element method to solve the poroelastic equations the main challenge is to ensure convergence of the method and prevent numerical instabilities that often manifest themselves in the form of spurious oscillations in the pressure field. It has been suggested that this problem is caused by the saddle point structure in the coupled equations resulting in a violation of the famous Ladyzhenskaya-Babuska-Brezzi (LBB) condition, thus highlighting the need for a stable combination of mixed finite elements [23].
In addition, there is a need for methods that do not give rise to localised pressure oscillations when seeking to approximate steep pressure gradients in the solution. For example, when modelling the diseased lung, abrupt changes in tissue properties and heterogeneous airway narrowing are possible. This can result in a patchy ventilation and pressure distribution [51]. In this situation methods that solve the poroelastic equations using a continuous pressure approximation strug-gle to capture the steep gradients in pressure and produce localised oscillations in the pressure [46]. Steep pressure gradients can also result from imposed Dirichlet pressure boundary conditions such as those in Terzaghi's problem [43,54]. The method presented here is able to overcome these types of pressure instability.

Two variable versus three variable formulations
The linear (infinitesimal strain) poroelastic equations are often solved in a displacement and pressure formulation from which the fluid flux can be recovered [43,54]. The stability and convergence of this displacement and pressure (u/ p) formulation was analysed in [43] and error bounds for inf-sup stable combinations of finite element spaces (e.g. Taylor-Hood elements) were obtained. In the current work we maintain the fluid flux as a variable, resulting in a three-field, displacement, fluid flux and pressure formulation. Retaining the fluid flux as a primary variable has the following advantages.
1. It allows for greater accuracy in the fluid velocity field.
This can be of particular interest when a poroelastic model is coupled with an advection diffusion equation, e.g., to account for gas exchange, thermal effects, contaminant transport or the transport of nutrients or drugs within a porous tissue [30]. 2. Physically meaningful boundary conditions can be applied at the interface when modelling the interaction between a fluid and a poroelastic structure [5]. 3. It allows for an easy extension of the fluid model from a Darcy to a Brinkman flow model, for which there are numerous applications in modelling biological tissues [30]. 4. It avoids the calculation of the fluid flux in postprocessing.
A clear disadvantage of a three-field formulation is the increased number of degrees of freedom of the linear system arising from the FEM discretisation, although this difficulty is mitigated with the proposed element, see comment 7 in §1.4.

Previous results: infinitesimal strain
Error estimates for finite element solutions of the linear three-field problem, using continuous piecewise linear approximations for displacements and mixed low-order Raviart-Thomas elements for the fluid flux and pressure variables, are presented in [44,45]. However this method was found to be susceptible to spurious pressure oscillations [47].
In an effort to overcome these pressure oscillations, a discontinuous linear three-field method was analysed in [38] with moderate success, and a linear non-conforming three-field method was analysed in [58]. However no implementation of these methods in 3D has yet been presented. Due to the size of the discrete systems resulting from a three field approach, there has been considerable work on operating splitting (iterative) approaches in which the poroelastic equations are separated into a fluid problem and deformation problem [19,28,32,53]. These methods are often able to take advantage of existing finite element software for elasticity and fluid flow. Matrix assembly for discontinuous and non-conforming finite elements in 3D can be complicated and calculating stresses using these methods can be particularly challenging. Methods that use standard and simple to implement elements are very appealing [55]. In [28], a linear three-field mixed finite element method using lowest order Raviart-Thomas elements was shown to overcome Dirichlet boundary type pressure instabilities.
The finite volume method has been used by [31,32] to discretize the flow. This results in a discontinuous approximation of the fluid pressure which is able to overcome localised pressure oscillations due to steep pressure gradients in the solution.
Introducing a displacement stress field, e.g. see [49], reduces the regularity requirements on the displacement field, thereby allowing for the implementation of a four-field conforming Raviart-Thomas element, but consequently greatly increasing the overall size and complexity of the problem.

Previous results: finite strain
Monolithic approaches for solving the quasi-static two-field incompressible finite-strain deformation equations are outlined in [1]. Two different approaches are advocated, a mixed-penalty formulation, in which the continuity condition is imposed using a penalty approach, and a mixed solid velocity-pressure formulation, where the linear momentum for the fluid is used to eliminate the fluid velocity in the remaining equations. The solid velocity-pressure formulation is similar to the commonly used reduced (u/ p) formulation in [4]. Two-field formulations require a stable mixed element pair such as the popular Taylor-Hood element to satisfy the LBB inf-sup stability requirement. However, using a continuous pressure element means that jumps in material coefficients may introduce large solution gradients across the interface, requiring severe mesh refinement or failing to reliably capture jumps in the pressure solution [54]. An operator splitting (iterative) approach for a near incompressible model is described by [13].
A three-field (displacement, fluid flux, pressure) formulation has been outlined in [37], however this method uses a low-order mixed finite element approximation without any stabilisation and therefore is not inf-sup stable. A three-field finite element using a continuous pressure approximation has been implemented in [52].
For both two-field and three-fields formulations and for any choices of finite element, implementation, construction and linearization of the nonlinear equations and convergence of the nonlinear mechanics problem using Newton's method or other iterative procedure is also nontrivial [50].

Contributions of the current work
In [7], we developed a stabilized, low-order, three-field mixed finite element method for the fully saturated, incompressible, small deformation case for which a linearly elastic model is sufficient. Low-order finite element methods are relatively easy to implement and allow for efficient preconditioning [20,27,55]. Rigorous theoretical results for the stability and optimal convergence rate for linear poroelasticity were presented in [7]. The stabilization term requires only a small amount of additional computational work and can be assembled locally on each element using standard finite element information, leading to a symmetric addition to the original system matrix and preserving any existing symmetry. The effect of the stabilization on the conservation of mass is minimal in 3D, and decreases as the mesh is refined, see [7].
In the current work we present a monolithic mixed finite element method to solve the fully incompressible three-field finite-strain poroelasticity equations. We use a low order approximation, namely, continuous piecewise-linear approximation for the displacements and fluid flux, and piecewiseconstant approximation for the pressure. The finite-strain case is a non-trivial extension of the infinitestimal-strain implementation and requires several challenges to be overcome including: (1) problem formulation and construction of the weak equations, (2) accurate integration over the deformed domain, (3) linearization of the weak equations, (4) construction and convergence of a Quasi-Newton iteration, and (5) time step selection. The main contributions of this work are therefore as follows.
1. A method to solve finite-strain fully incompressible poroelasticity using a stabilized discontinuous pressure approximation. [Note that for the linear (infinitesimal strain) equations other methods that use a discontinuous pressure approximation have been previously presented [7,28,38,58].] 2. A method for finite-strain fully incompressible poroelasticity that is both inf-sup stable and is able to overcome localized pressure oscillations. 3. A finite element method that is robust within all modelling regimes. Large differences in permeability within the computational domain can result in regions in which Darcy flow dominates over elastic effects and regions in which elastic effects are dominant. This low-order ele-ment is reliable in both scenarios, providing an effective numerical approach for problems in which heterogeneity presents computational challenges. 4. A finite element method that results in a discrete system with blocks arising from simple linear finite elements allowing fast solver approaches and preconditioning techniques to be easily implemented. 5. A finite element method for a finite-strain poroelastic model that resolves steep pressure gradients without localized oscillations.
We present a quasi-static finite-strain incompressible poroelastic model in Sect. 2 and develop the stabilized nonlinear finite-element method in Sect. 3. Implementation details are provided in Sect. 4 and in the appendices. In Sect. 5, we present a range of numerical experiments to verify the accuracy of the method and to demonstrate its ability to reliably capture steep pressure gradients.

Poroelasticity theory
Two complementary approaches have been developed for modelling a deformable porous medium. Mixture theory, also known as the Theory of Porous Media (TPM) [8,10,11], has its roots in the classical theories of gas mixtures and makes use of a volume fraction concept in which the porous medium is represented by spatially superposed interacting media. An alternative, purely macroscopic approach is mainly associated with the work of Biot. A comprehensive development of the macroscopic theory appears in [16]. Relationships between the two theories are explored in [14,17]. As is most common in biological applications, we use the mixture theory for poroelasticity as outlined in [8] and recently summarized in [52].

Kinematics
Let the volume Ω(0) be the undeformed Lagrangian (material) reference configuration and let X indicate the position of a particle in Ω(0) at t = 0. The position of a particle in the deformed configuration Ω(t) at time t > 0 is given by x, with x = χ (X, t) as shown in Fig. 1. The deformation map, χ(X, t), is a continuously differentiable, invertible mapping from Ω(0) to Ω(t). Thus the inverse of the deformation map, χ −1 (x, t), is such that X = χ −1 (x, t). The displacement field is given by The deformation gradient tensor is p(X, t) p (x, t) and the symmetric right Cauchy-Green deformation tensor is The Jacobian is defined as and represents the change in an infinitesimal control volume from the reference to the current configuration, i.e., Note that J > 0.

Volume fractions
We will only consider saturated porous media in which the fluid accounts for volume fractions φ 0 (X, t = 0) and φ(x, t) of the total volume in the reference and current configurations respectively, where φ is known as the porosity, defined through the Jacobian of the deformation (7). The fractions for the solid (or skeleton) are therefore 1 − φ 0 and 1 − φ in the reference and current configuration, respectively. For the mixture, ρ is the density in the current configuration given by where ρ s and ρ f are the densities of the fluid and solid, respectively. We assume that both the solid and the fluid are incompressible so that ρ s = ρ s 0 and ρ f = ρ f 0 . Although both the solid and fluid are assumed to be incompressible, the control volume can expand or contract due to fluid entering or leaving the region, and

The model
We define the boundary for the fluid, with an outward pointing unit normal n. We seek deformation χ(X, t), fluid flux z(x, t) and pressure p(x, t) such that in Ω(0).
where v f and v s are the velocities of the fluid and solid components respectively, χ t denotes ∂χ(X,t) ∂t , u D , q D , p D are given boundary conditions, f is a general external body force, g is a general source or sink term and σ e is the stress tensor given by where W (C), with C = F T F, denotes a strain-energy law (hyperelastic Helmholtz energy functional) dependent on the deformation of the solid. The permeability tensor is given by where k 0 (C) is the permeability in the reference configuration, which may be chosen to be some (nonlinear) function dependent on the deformation. Examples of deformation dependent permeability tensors for biological tissues can be found in [24,34,35]. Details of the derivation of (8) appear in "Appendix 1". It is important to recognize that ∇(·) = ∂/∂ x(·) denotes the partial derivative with respect to the deformed configuration. We will use ∇ to denote the spatial gradient in Ω(t) rather than the more explicit ∇ x=χ(X,t) . The latter more clearly indicates the dependency of the gradient operator on the deformation χ(X, t) and highlights the inherent nonlinearity that arises due to the fact that the deformation χ (X, t) is one of the unknowns. Similarly the deformed domain Ω(t) in which equations (8) pertain, is a function of the deformation map χ , and therefore incorporates another important nonlinearity.

The stabilized finite element method
We extend the method of [7] from the linear, small deformation poroelastic case to finite-strain poroelasticity. For ease of presentation, we will assume all Dirichlet boundary conditions are homogeneous, ie., u D = 0, q D = 0, p D = 0.

Weak formulation
We respectively define the following spaces for the deformed location, fluid flux and pressure, , Here ∇ S v = 1 2 ∇v + (∇v) T for some vector v.

The fully discrete model
Let T h be a quasi-uniform partition of Ω(t) into nonoverlapping elements K , where h denotes the size of the largest element in T h . We then define the following finite element spaces, where P 0 (K ) and P 1 (K ) are the spaces of constant and linear polynomials on K respectively. We define the combined solu- Δt . The fully discrete weak problem is: , The stabilization term is given by where Υ is a stabilization parameter that is independent of h and Δt. Here h ∂ K denotes the size (diameter) of an element edge in 2D or face in 3D, and · is the jump across an edge or face (taken on the interior edges only). The stabilization term has been introduced here to add stability and ensure a well-posed fully-discrete model. It has been shown that the convergence is insensitive to Υ , e.g. see in [7,12,27] .
). The nonlinear system of Eq. (12) can be recast in the form: Find where Given an approximate solution u n h , we approximate (13) by and solve for the Newton step δu h , where DG is the directional derivative of G, at u n h , in the direction δu h .

Approximation of DG n .
In biphasic tissue problems, it is common to approximate directional derivative of G by assuming the nonlinear elasticity term is the dominant nonlinearity and ignoring the other nonlinearities [50,54]. Let For Newton's method we require the directional derivative of E n ((χ n h , p n h ), v h ) at a particular trial solution (χ n h , p n h ) in the direction δχ h , given by (see [57, where n h is a fourth-order tensor and σ n e,h is the effective (elastic) stress tensor, both evaluated at a trial solution χ n h . Further, any variable with a bar above it will correspond to it being evaluated at a trial solution. The fourth-order spatial tangent modulus tensor is described in "Appendix 2". For a detailed explanation and derivation see [9,57]. The approximate linearization of the nonlinear problem is thus given by Using (14), (18) and Eq. (15) the Newton solve becomes: = , , (Ω(t n )).

Matrix assembly for the Newton iteration
Let φ k denote a vector-valued linear basis function for the (P1) d space, and Similarly let ψ i denote a basis function for the space P0, hence Now let u n i := (χ n i , z n i , p n i ) ∈ R n u +n z +n p denote the fully discrete solution at the ith step within the Newton method at time t n . The Newton algorithm at a particular time step n, is given in Algorithm 1.

Algorithm 1 Newton algorithm at
Update the mesh, Ω(t n ) i+1 = χ n i i = i + 1 end while At each Newton iteration we are required to solve the linear system This system can be expanded as where the elements in the matrices in (21) are given by The saddle point system given in Eq. (21) can be iteratively solved using similar approaches to those seen in [29,48], where the action of the inverse of the preconditioner fundamentally requires only the inverse of the linear blocks K e and M. Details of the matrices D and E appear in "Appendix 2". Note that the matrix equations are integrated in the deformed configuration obtained from the previous Newton step. This update Lagrangian approach overcomes complex linearisation otherwise needed when using a total Lagrangian approach [50]. The Newton iteration was found to be robust with respect to the stabilization parameter. More precisely, in all calculations fewer than four Newton iterations were required independently of the size of the stabilization parameter. Tables 2 and 3 in Sect. 5.1 show the Newton convergence for two choices of the stabilization parameter two orders of magnitude apart for a 3D stress relaxation text problem. In practice, the stabilization parameter was chosen so as to be as small as possible without producing oscillations, unless otherwise stated.

Stabilization matrix assembly
Let K ∈ T h be an element and D(K ) be the pressure degree of freedom associated with element K . We define A(K ) to be the set of elements L ∈ T h neighboring K .

Fluid-flux boundary condition
When solving the equations for Darcy flow using the Raviart-Thomas element (RT-P0), the fluid-flux boundary condition is enforced naturally by this divergence free element. Unfortunately this is not possible using our proposed P1-P1-P0-stabilized element. However, solving the poroelastic equations (8) using a piecewise linear approximation for the deformation and Raviart-Thomas element for the fluid (P1-RT-P0) does not satisfy the discrete inf-sup condition and can yield spurious pressure oscillations, see [46,47] for details.
To enforce the flux boundary condition z · n = q D along the boundary Γ f (t) we introduce a Lagrange multiplier Λ h , where Λ h ∈ W f h (t), the discrete space of piecewise constant functions defined on all element surfaces with non-zero intersection with Γ f (t). The resulting modified continuous weak-form is The discretization and implementation of this additional constraint is straightforward and results in a discrete system with additional degrees of freedom for every node on Γ f . The terms (Λ h , w · n) Γ f and (z · n, l) Γ f are nonlinear since the normal is a function of the (nonlinear) displacement. Note, within all the simulations we have undertaken, we found that treating these terms as linear terms did not prevent the convergence of the Newton algorithm. Alternatively these terms could be linearized as has been described in detail for the traction boundary condition, see [57, section 4.2.5] and [4].

Numerical results
We present four numerical examples to test the performance of the proposed stabilized finite element method. The first two examples are biological and geomechanical applications respectively, and the third is a swelling example that undergoes significant, large deformations. For the implementation we used the C++ library libmesh [33], and the multi-frontal direct solver mumps [2] to solve the resulting linear systems. For the strain energy law we chose a Neo-Hookean law taken from [57, eqn. (3.119)], with the penalty term chosen such that 0 ≤ φ < 1, namely For further discussion of strain energy laws for porelasticity we refer to [14] and [52]. The material parameters μ and λ in (24) can be related to the Young's modulus E and the Poisson ratio ν by μ = E/(2(1 + ν)) and λ = (Eν)/((1 + ν)(1 − 2ν)). Details of the effective stress tensor and fourth-order spatial tangent modulus for this particular law can be found in "Appendix 2". For the permeability law we chose

3D unconfined compression problem
This first example tests the correctness of the implementation by comparing the numerical solution of the finite-strain system of equations to one of the very few available analytical solutions for poroelastic problems, all-be-it for small deformations. Similar unconfined compression problems have been used to test other large deformation poroelastic software such as FEBio [40]. A cylinder of poroelastic materials is subjected to a prescribed displacement in the axial direction. The material is allowed to relax in the radial direction by constraining the fluid pressure to be zero at the outer radial surface and assuming the outer radial boundary is permeable and free-draining. The upper and lower fluid boundaries are assumed to be impermeable and frictionless. The original experiment involved a specimen of articular cartilage being compressed via impervious smooth plates as shown in Fig. 2. Both the radius and height of the cylinder are 1 mm, whereas the axial compression is 0.01 mm, hence the axial compression is small compared to the size of the height of the cylinder. Large deformation effects are therefore expected to be negligible. The parameters used for the simulation can be found in Table 1  where u is the radial displacement, 0 is the amplitude of the applied axial strain and a is the radius of the cylinder. Here α n are the solutions to the characteristic equation, given by where J 0 and J 1 are Bessel functions. The characteristic time of diffusion t g is given by t g = a 2 /Mk, where M = λ + 2μ is the P-wave modulus of the elastic solid skeleton, k is the permeability. For small axial compression the computed radial displacement shown in Fig. 3 is in good agreement with the analytical solution, indicating that the nonlinear poroelastic model is accurate in the small strain limit. As the axial compression becomes large, the numerical finite strain solution departs from the analytical linear small deformation solution as expected.
The effect of the stabilization parameter on the numerical solution is investigated in Fig. 4, and shown to be robost for a broad range of values, since the stabilization parameter can be chosen to be very small in 3D. Without stabilization the Newton method solving the non-linear equations diverged rapidly due to spurious pressure modes present at each Newton step and a solution could not be obtained. Further tests to investigate this type of loss of stability are given in Sect.

5.3.
Tables 2 and 3 illustrate convergence of the Newton iteration for the unconfined compression problem for the first time step (the most demanding due to the initial displacement boundary condition), with Υ = 10 −1 and Υ = 10 −3 ,

Terzaghi's problem
Terzaghi's problem is a common test problem within the geomechanics community that has an analytical solution. It has been used to investigate the origins of non-physical pressure oscillations arising in some finite element solutions near the boundary [43,54]. The domain consists of a porous column of unit height, bounded at the sides and bottom by rigid and impermeable walls. The top is free to drain ( p D = 0) and has a downward traction force, p 0 , applied to it. The boundary and initial conditions for this 1D problem can be written as The analytical pressure solution, in non-dimensional form is given by For a detailed explanation and derivation of this solution see [16, section 5.2.2]. We discretized the column using 60 hexahedral elements and solved the problem using both the stabilized low-order finite element method and a higher-order inf-sup stable finite element method with piecewise linear pressure approximation. The material parameters used for the simulation can be found in Table 4. The simulation results of the pressure for the two methods at t = 0.01 s and t = 1 s are shown in Fig. 5. At t = 0.01 s the piecewise linear (continuous) approximation, which is inf-sup stabilized using a Brinkman term [18], fails to approximate the thin boundary layer in the pressure field and suffers from overshooting (Fig. 5a). The stabilized low-order method does not suffer from this problem and accurately captures the pressure field near the boundary (Fig. 5c). At t = 1 s the boundary layer has grown and both the piecewise linear pressure approximation (Fig. 5b) and the piecewise constant pressure approximations (Fig. 5d) yield satisfactory results. Note that even when using discontinuous pressure interpolation, pressure oscillations are present without stabilization at t = 0.01 s, see Fig. 5e. Again this pressure oscillation disappears as the pressure boundary layer grows with time and the lack of inf-sup stability is not obvious from the solution at t = 1 s, see Fig. 5f.

Swelling test
Given a unit cube of material, a fluid pressure gradient is imposed between the two opposite faces at X = 0 and X = 1. The pressure p D on the inlet face X = 0 is increased very rapidly from zero to a limiting value of 10kPa, i.e., p D = 10 4 (1 − exp(−t 2 /0.25)) Pa). On the outlet face X = 1, the pressure is fixed to be zero, p D = 0. There are no sources of sinks of fluid. A zero flux condition is applied for the fluid velocity on the four other faces (Y = 0, 1, Z = 0, 1). Normal displacements are required to be zero on the planes X = 0, Y = 0 and Z = 0. The permeability of the cube 0 < X < 0.5, 0.5 < Y < 1, 0 < Z < 0.5, i.e., 1/8 of the volume, is smaller than in the rest of the unit cube by a factor of 500. The computational domain is shown in Fig. 6a, highlighting the region of reduced permeability. The parameters chosen for this test problem are given in Table  5. This problem is similar to the one in [13] and highlights the method's ability to reliably capture steep gradients in the pressure solution due to rapid changes in material parameters.
Fluid enters the region from the inlet face and the material swells like a sponge, undergoing large deformation as shown  in Fig. 6b. The evolution of the pressure and the Jacobian at the points at (0, 0, 1), (0.5, 0, 1) and (1, 0, 1) in the reference configuration are shown in Fig. 7a and b respectively. These positions are indicated by the red, blue and green balls in Fig.  6a. The pressure and volume change at the point (0, 1, 0) (black ball in Fig. 6a) is also shown in Fig. 7a and b. Due to its reduced permeability, this region is much slower to swell and achieve its ultimate equilibrium state and the fluid mainly flows around the region of reduced permeability, see Fig. 6b. The steep pressure gradients at the boundary of the less permeable region seen in Fig. 6b are well approximated by the piecewise constant (discontinuous) pressure space even on this relatively coarse discretization, and the no-flux boundary condition is enforced correctly along the deformed boundary. Continuous pressure spaces would require a much finer discretization in this region. Figure 8 shows the pressure solution for this test problem with (Υ = 10 −4 in Fig. 8a) and without stabilization (Υ = 0 in Fig. 8b) at t = 0.02 s. The computation was performed using 512 hexahedral elements. Note that without any stabilization pressure an instability in the pressure is observed.
To further investigate and demonstrate that a lack of stabilization will result in a loss of inf-sup stability and thus result in a spurious chequer board pressure solution, we run the same swelling test problem, but with a homogeneous permeability permeability set at k 0 = 10 −5 . Furthermore we solve this problem on a 12,288 element tetrahedral mesh, which has a smaller ratio of fluid and displacement nodes to pressure nodes, compared to the previously used hexahedral mesh, thus worsening the inf-sup instability properties [18]. Figure 9 shows the pressure solution with (Υ = 10 −4 in Fig.  9a) and (Υ = Υ = 10 −12 in Fig. 9b) at t = 0.02 s. Note that without sufficient stabilization the lack of inf-sup stability can clearly be observed in the form of a spurious pressure checkerboard solution. These numerical examples demonstrate that the stabilization scheme is very robust to ensuring inf-sup stability by allowing a large range of stabilization parameters to be used, before spurious pressure oscillations caused by a loss of inf-sup stability are observed. However when wishing to capture pressure boundary layer type solutions, care does need to be taken to ensure that the pressure solution is not overly smoothed. As a practical guide we recommend first choosing a large stabilization parameter and then repeatedly lowering the stabilization parameter by an order of magnitude until the lowest parameter is found that leads to a pressure solution without any oscillations.

Conclusions
Stabilized low-order methods can offer significant computational advantages over higher order approaches. In particular, one can employ meshes with fewer degrees of freedom, fewer Gauss points, and simpler data structures. The additional stabilization terms can also improve the convergence properties of iterative solvers.
The main contribution of this paper is to extend the local pressure jump stabilization method [12] already applied to three-field linear poroelasticity in [7], to the finite strain case. Thus, the proposed scheme is built on an existing scheme for which rigorous theoretical results addressing the stability and optimal convergence have been proven, and for which numerical experiments have demonstrated its ability to overcome spurious pressure oscillations. Owing to the discontinuous pressure approximation, sharp pressure gradients due to changes in material coefficients or boundary layers can be captured reliably, circumventing the need for severe mesh refinement. The addition of the stabilization term introduces minimal additional computational work, can be assembled locally on each element using standard element information, and leads to a symmetric addition to the original system matrix, thus preserving any existing symmetry. As the numerical examples have demonstrated, the stabilization scheme is robust and leads to high-quality solutions. (30) where v f is the velocity of the fluid and v s is the velocity of the solid given by v s (x, t)| x=χ(X,t) = ∂χ(X, t) ∂t . (31) and g is a general source or sink term. In differential form, or, First noting that ρ s and ρ f are constants and can be factored out and then adding Eqs. (32) and (33), provides the mass balance or continuity equation of the mixture,

Conservation of momentum
For α = s, f the conservation of linear momentum for solid and fluid components is given by where σ α is the Cauchy stress tensor for the α = s, f , f is a body force,p α are interaction forces representing frictional interactions between the solid and fluid (see 1) and β α is the constituent source term. Here β s = 0 and β f = ρ f g. Applying Reynolds Transport Theorem, we rewrite the integral conservation law in differential form and obtain Expanding the LHS, Using (34) and (35) to replace the second term above, ρ α a α = ∇ · σ α +ρ α f +p α in Ω(t), where a s (x, t)| x=χ(X,t) = ∂ 2 χ (X, t)

Constitutive relationships
Constitutive relationships for the interaction forces, permeability tensor and solid and fluid stress tensors are provided below. We choose expressions for the constitutive laws that appear frequently in the literature.

Interaction forces
The interaction force is given bŷ where k is the (dynamic) permeability tensor and p is the fluid pressure [16]. The first term, p∇φ, accounts for the pressure effect resulting from the variation of the section offered to the fluid flow, and the second term, φ 2 k · (v f − v s ), describes the viscous resistance opposed by the shear stress to the fluid flow from the drag at the internal walls of the porous network. This particular choice for the interaction force means that the momentum balance for the fluid flow can later be reduced to the well known Darcy law.

Permeability tensor
The permeability tensor is given by where k 0 (C) is the permeability in the reference configuration, which may be chosen to be some (nonlinear) function dependent on the deformation. Examples of deformation dependent permeability tensors for biological tissues can be found in [24,34,35]. A common isotropic assumption is where κ 0 is the permeability in the reference configuration and Π (J ) is some function dependent on the volume change.
For example, in [34], the following isotropic constitutive law for the permeability of lung tissue is proposed where κ 0 is the permeability in the reference configuration.

Solid stress tensor
The solid stress tensor is given by [8], where σ s e is the effective stress tensor given by Here W (C) denotes a strain-energy law (hyperelastic Helmholtz energy functional) dependent on the deformation of the solid.

Fluid stress tensor
The fluid stress tensor can be written as [8], where σ f vis denotes the viscous stress tensor of the fluid, given by where μ f is the dynamic viscosity of the fluid.

The general poroelasticity model
Summing the conservation laws for solid and fluid and applying the constitutive relations, the conservation of linear momentum for the mixture iŝ ρ s a s +ρ f a f = ∇ · (σ e + σ vis − p I) + ρ f in Ω(t). (50) The momentum equation for the fluid flow alone iŝ We define the boundary ∂Ω(t) = Γ d (t) ∪ Γ n (t) for the mixture and ∂Ω(t) = Γ p (t)∪Γ f (t) for the fluid, with an outward pointing unit normal n. The problem for the mixture theory model is : Find χ(X, t), v f (x, t) and p(x, t) such that ρ s a s +ρ f a f = ∇ · (σ e + σ vis − p I) + ρ f in Ω(t), See [9, chapter 5] and [57, chapter 3] for further details.
To simplify the implementation of the spatial tangent modulus we make use of matrix Voigt notation. The matrix form of β is given by D, which can be written as (see [ We also make use of the following implementation friendly matrix notation for ∇ S φ k