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

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

Physically meaningful boundary conditions can be applied at the interface when modelling the interaction between a fluid and a poroelastic structure [5].

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

It avoids the calculation of the fluid flux in post-processing.

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 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 piecewise-constant 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.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.

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

A method for finite-strain fully incompressible poroelasticity that is both inf-sup stable and is able to overcome localized pressure oscillations.

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 element is reliable in both scenarios, providing an effective numerical approach for problems in which heterogeneity presents computational challenges.

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.

A finite element method for a finite-strain poroelastic model that resolves steep pressure gradients without localized oscillations.

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 be the undeformed Lagrangian (material) reference configuration and let indicate the position of a particle in at . The position of a particle in the deformed configuration at time is given by , with as shown in Fig. 1. The deformation map, , is a continuously differentiable, invertible mapping from to . Thus the inverse of the deformation map, , is such that . The displacement field is given byThe deformation gradient tensor isand the symmetric right Cauchy–Green deformation tensor isThe Jacobian is defined asand represents the change in an infinitesimal control volume from the reference to the current configuration, i.e.,Note that .

Fig. 1

Illustration of the solid deformation

Volume fractions

We will only consider saturated porous media in which the fluid accounts for volume fractions and 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 and in the reference and current configuration, respectively. For the mixture, is the density in the current configuration given bywhere and are the densities of the fluid and solid, respectively. We assume that both the solid and the fluid are incompressible so that and . 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 mixture and for the fluid, with an outward pointing unit normal . We seek deformation , fluid flux and pressure such thatThe fluid flux where and are the velocities of the fluid and solid components respectively, denotes , are given boundary conditions, is a general external body force, g is a general source or sink term and is the stress tensor given bywhere , with , denotes a strain-energy law (hyperelastic Helmholtz energy functional) dependent on the deformation of the solid. The permeability tensor is given bywhere 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 denotes the partial derivative with respect to the deformed configuration. We will use to denote the spatial gradient in rather than the more explicit . The latter more clearly indicates the dependency of the gradient operator on the deformation and highlights the inherent nonlinearity that arises due to the fact that the deformation is one of the unknowns. Similarly the deformed domain 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., .

Weak formulation

We respectively define the following spaces for the deformed location, fluid flux and pressure,where .

The continuous weak problem is: Find , and for any time such that Here for some vector .

The fully discrete model

Let be a quasi-uniform partition of into non-overlapping elements K, where h denotes the size of the largest element in . We then define the following finite element spaces,where and are the spaces of constant and linear polynomials on K respectively. We define the combined solution space .

The discretization in time is given by partitioning [0, T] into N evenly spaced non-overlapping regions , , where . For any sufficiently smooth function v(t, x) we define and the discrete time derivative by .

The fully discrete weak problem is: For , find , and such thatThe stabilization term is given bywhere is a stabilization parameter that is independent of h and . Here 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] .

Solution via quasi-Newton iteration at .

Let denote the solution vector at a particular time step, denote the solution increment vector, and where . The nonlinear system of Eq. (12) can be recast in the form: Find such thatwhereGiven an approximate solution , we approximate (13) byand solvefor the Newton step , where DG is the directional derivative of G, at , in the direction .

Approximation of .

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]. LetFor Newton’s method we require the directional derivative of at a particular trial solution in the direction , given by (see [57, section 3.5.3])where is a fourth-order tensor and is the effective (elastic) stress tensor, both evaluated at a trial solution . 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 byUsing (14), (18) and Eq. (15) the Newton solve becomes: Find , and such that

Implementation details

Matrix assembly for the Newton iteration

Let denote a vector-valued linear basis function for the (P1) space, andSimilarly let denote a basis function for the space P0, henceNow let denote the fully discrete solution at the ith step within the Newton method at time . The Newton algorithm at a particular time step n, is given in Algorithm 1.

At each Newton iteration we are required to solve the linear systemThis system can be expanded aswhere the elements in the matrices in (21) are given byThe 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 and M. Details of the matrices and 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.

Table 2

Convergence of the Newton iteration for the 3D unconfined compression problem with at

Newton iteration	||uin-ui-1n||	||R(uin,u^n-1)||
1	0.81	0.023202
2	2.81699e−04	0.011276
3	6.93986e−08	1.34048e−06
4	4.10726e−10	7.64882e−09

Table 3

Convergence of the Newton iteration for the 3D unconfined compression problem with at

Newton iteration	||uin-ui-1n||	||R(uin,u^n-1)||
1	0.81	0.0235609
2	1.28528e−05	1.99541e−04
3	7.71658e−08	1.49304e−06
4	4.6844e−10	8.17681e−09

Stabilization matrix assembly

Let be an element and be the pressure degree of freedom associated with element K. We define to be the set of elements 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 along the boundary we introduce a Lagrange multiplier , where , the discrete space of piecewise constant functions defined on all element surfaces with non-zero intersection with . 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 . The terms and 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 , namelyFor 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 and . 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. All computations were performed using 3080 tetrahedral elements.

Fig. 2

The 3D unconfined compression problem

Table 1

Parameters used for the 3D unconfined compression test problem

Parameter	Description	Value
ϕ_0	Initial fluid volume fraction	0.9
k_0	Dynamic permeability	10^-3m^3skg^-1
ν	Poisson’s ratio	0.15
E	Young’s modulus	1000 kgm^-1s^-2
Δt	Time step used in the simulation	4s
T	Final time of the simulation	1000s
Υ	Stabilization parameter	10^-3

A closed-form axisymmetric solution for small strains [3] iswhere u is the radial displacement, is the amplitude of the applied axial strain and a is the radius of the cylinder. Here are the solutions to the characteristic equation, given by where and are Bessel functions. The characteristic time of diffusion is given by , where 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.

Fig. 3

The 3D unconfined compression problem. Normalized radial displacement versus normalized time for vertical normalized displacements compared to the analytical, infinitesimal strain solution. All computations performed with

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.

Fig. 4

The 3D unconfined compression problem. Normalized radial displacement versus normalized time calculated using stabilized finite element method for using various values of and compared to the analytical, infinitesimal strain solution

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 and , respectively. The Newton convergence is minimally affected by the decrease in stabilization. The Newton iteration fails to converge if the stabilization parameter is further reduced to , and as a consequence no solution can be obtained. The linear system contains 8162 degrees of freedom, takes to assemble and to solve, using one Intel Xenon CPU.

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 () and has a downward traction force, , applied to it. The boundary and initial conditions for this 1D problem can be written asThe analytical pressure solution, in non-dimensional form is given byFor 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.

Table 4

Parameters used for Terzaghi’s problem

Parameter	Description	Value
ϕ_0	Initial fluid volume fraction	0.9
k_0	Dynamic permeability	10^-5m^3skg^-1
ν	Poisson ratio	0.25
E	Young’s modulus	100 kgm^-1s^-2
Δt	Time step used in the simulation	0.01s
T	Final time of the simulation	1s
Υ	Stabilization parameter	2×10^-5

The simulation results of the pressure for the two methods at and are shown in Fig. 5. At 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 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 , 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 , see Fig. 5f.

Fig. 5

Terzaghi’s problem. a Pressure at using a continuous linear pressure approximation. b Pressure at using a continuous linear pressure approximation. c Pressure at using a discontinuous piecewise constant approximation with . d Pressure at using a discontinuous piecewise constant approximation with . e Pressure at using a discontinuous piecewise constant approximation without stabilization. f Pressure at using a discontinuous piecewise constant approximation without stabilization.

Fig. 7

Swelling test. a Pressure, p, at locations [red], [blue], [green] and [black]. b Volume change, J (b) at locations [red], [blue], [green] and [black]. (These locations are shown using the colored balls in Fig. 6a)

Swelling test

Given a unit cube of material, a fluid pressure gradient is imposed between the two opposite faces at and . The pressure on the inlet face is increased very rapidly from zero to a limiting value of , i.e., . On the outlet face , the pressure is fixed to be zero, . There are no sources of sinks of fluid. A zero flux condition is applied for the fluid velocity on the four other faces (). Normal displacements are required to be zero on the planes and . The permeability of the cube , 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.

Fig. 6

Swelling test. a Initial configuration. The grey cube represents the region of reduced permeability. The colored balls indicate the position of the points used for tracking the pressure and volume changes shown in Fig. 7a and b. b The deformed cube after 1 s showing the pressure solution and fluid flux

Table 5

Parameters used for the swelling test problem

Parameter	Description	Value
ϕ_0	Initial fluid volume fraction	0.9
k_0	Dynamic permeability	10^-5m^3skg^-1
ν	Poisson ratio	0.3
E	Young’s modulus	8000 kgm^-1s^-2
Δt	Time step used in the simulation	0.02s
T	Final time of the simulation	20s
Υ	Stabilization parameter	10^-4

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 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 ( in Fig. 8a) and without stabilization ( in Fig. 8b) at . The computation was performed using 512 hexahedral elements. Note that without any stabilization pressure an instability in the pressure is observed.

Fig. 8

Swelling test problem with non-uniform permeability. Pressure field at using the stabilized finite element method. a , b

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 . 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 ( in Fig. 9a) and ( in Fig. 9b) at . 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.

Fig. 9

Swelling problem with uniform permeability. Tetrahedral mesh. a Pressure field at using the stabilized finite element method. a , b

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.