New stabilized discretizations for poroelasticity and the Stokes’ equations
Abstract
In this work, we consider the popular P1RT0P0 discretization of the threefield formulation of Biot’s consolidation problem. Since this finiteelement formulation does not satisfy an infsup condition uniformly with respect to the physical parameters, several issues arise in numerical simulations. For example, when the permeability is small with respect to the mesh size, volumetric locking may occur. Thus, we propose a stabilization technique that enriches the piecewise linear finiteelement space of the displacement with the span of edge/face bubble functions. We show that for Biot’s model this does give rise to discretizations that are uniformly stable with respect to the physical parameters. We also propose a perturbation of the bilinear form, which allows for local elimination of the bubble functions and provides a uniformly stable scheme with the same number of degrees of freedom as the classical P1RT0P0 approach. We prove optimal stability and error estimates for this discretization. Finally, we show that this scheme can also be successfully applied to Stokes’ equations, yielding a discrete problem with optimal approximation properties and with minimum number of degrees of freedom (equivalent to a P1P0 discretization). Numerical tests confirm the theory for both poroelastic and Stokes’ test problems.
keywords:
Stable finite elements, poroelasticity, Stokes’ equations1 Introduction
The interaction between the deformation and fluid flow in a fluidsaturated porous medium is the object of study in poroelasticity theory. Such coupling has been modelled in the early onedimensional work of Terzaghi (1). A more general threedimensional mathematical formulation was then established by Maurice Biot in several pioneering publications (see (2) and (3)). Biot’s models are widely used nowadays in the modeling of many applications in different fields, ranging from geomechanics and petroleum engineering, to biomechanics. The existence and uniqueness of the solution for these problems have been investigated by Showalter in (4) and by Zenisek in (5). Regarding the numerical simulation of the poroelasticity equations, there have been numerous contributions using finitedifference schemes (6); (7) and finitevolume methods (see (8) for recent developments). Finiteelement methods, which are the subject of this work, have also been considered (see for example the monograph by Lewis and Schrefler (9) and the references therein).
Stable finiteelement schemes are constructed by either choosing discrete spaces satisfying appropriate infsup (or LBB) conditions, or applying suitable stabilization techniques to unstable finiteelement pairs. For the twofield (displacementspressure) formulation of Biot’s problem, the classical TaylorHood elements belong to the first class (10); (11); (12), as well as the MINI element (13). On the other hand, a stabilized discretization based on linear finite elements for both displacements and pressure has been recently analyzed in (13), and belongs to the second type. Regarding threefield formulations (including the Darcy velocity), a stable finiteelement method based on nonconforming CrouzeixRaviart finite elements for the displacements, lowest order RaviartThomasNédélec elements for the Darcy velocity, and piecewise constants for the pressure, was proposed in (14). For a fourfield formulation of the problem, which includes the stress tensor, the fluid flux, the solid displacement, and the pore pressure as unknowns, a stable approach is proposed in (15). In that work, two mixed finite elements, one for linear elasticity and one for mixed Poisson, are coupled for the spatial discretization.
This paper focuses on the threefield formulation, which has received a lot of attention from the point of view of novel discretizations (16); (17); (18); (19), as well as for the design of efficient solvers (20); (21); (22). Because of its application to existing reservoir engineering simulators, one of the most frequently considered schemes is a threefield formulation based on piecewise linear elements for displacements, RaviartThomasNédélec elements for the fluid flux, and piecewise constants for the pressure. This element, however, does not satisfy an infsup condition uniformly with respect to the physical parameters of the problem. Thus, we propose a stabilization of this popular element which gives rise to uniform error bounds, keeping the same number of degrees of freedom as in the original method.
A consequence of our analysis is that a new stable scheme for the Stokes’ equations is derived. The resulting method can be seen as a perturbation of the wellknown unstable pair based on piecewise linear and piecewise constant elements for velocities and pressure, respectively (P1P0). This yields a stable finiteelement pair for Stokes, which has the lowest possible number of degrees of freedom.
The rest of the paper is organized as follows. Section 2 is devoted to describing Biot’s problem and, in particular, the considered threefield formulation and its discretization. A numerical example is given, illustrating the difficulties that appear. In Section 3, we introduce the stabilized scheme in which we consider the enrichment of the piecewise linear continuous finiteelement space with edge/face (2D/3D) bubble functions. Section 4 is devoted to the local elimination of the bubbles to maintain the same number of degrees of freedom as in the original scheme. The wellposedness of the resulting scheme, as well as the corresponding error analysis are also provided here. In Section 5, we present the Stokesstable finiteelement method based on P1P0 finite elements obtained by following the same strategy as presented in the previous sections for poroelasticity. Finally, in Section 6, we confirm the uniform convergence properties of the stabilized schemes for both poroelasticity and Stokes’ equations through some numerical tests.
2 Preliminaries: model problem and notation
We consider the quasistatic Biot’s model for consolidation in a linearly elastic, homogeneous, and isotropic porous medium saturated by an incompressible Newtonian fluid. According to Biot’s theory (2), the mathematical model of the consolidation process is described by the following system of partial differential equations (PDEs) in a domain , with sufficiently smooth boundary :
equilibrium equation:  (1)  
constitutive equation:  (2)  
compatibility condition:  (3)  
Darcy’s law:  (4)  
continuity equation:  (5) 
where and are the Lamé coefficients, is the Biot modulus, and is the BiotWillis constant. Here, and denote the drained and the solid phase bulk moduli. As is customary, stands for the absolute permeability tensor, is the viscosity of the fluid, and is the identity tensor. The unknown functions are the displacement vector and the pore pressure . The effective stress tensor and the strain tensor are denoted by and , respectively. The percolation velocity of the fluid, or Darcy’s velocity, relative to the soil is denoted by and the vectorvalued function represents the gravitational force. The bulk density is , where and are the densities of solid and fluid phases and is the porosity. Finally, the source term represents a forced fluid extraction or injection process.
Our focus here is on the socalled threefield formulation in which Darcy’s velocity, , is also a primary unknown in addition to and . As a result, we have the following system of PDEs:
(6)  
(7)  
(8) 
This system is often subject to the following set of boundary conditions:
(9)  
(10) 
where is the outward unit normal to the boundary, , with and being open (with respect to ) subsets of with nonzero measure. In the following, we omit the symbol “” over and as it will be clear from the context that the essential boundary conditions are imposed on closed subsets of . Other sets of boundary conditions are also of interest, such as full Dirichlet conditions for both and on all of .
The initial condition at is given by,
(11) 
which yields the following mixed formulation
of Biot’s threefield consolidation model:
For each , find
such that
(12)  
(13)  
(14) 
where,
(15) 
corresponds to linear elasticity. The function spaces used in the variational form are
where is the space of square integrable vectorvalued functions whose first derivatives are also square integrable, and contains the square integrable vectorvalued functions with square integrable divergence.
We recall that the wellposedness of the continuous problem was established by Showalter (4), and, for the threefield formulation by Lipnikov (23). Next, we focus on the behavior of some classical discretizations of Biot’s model.
2.1 Discretizations
First, we partition the domain into dimensional simplices and denote the resulting partition with , i.e., . Further, with every simplex , we associate two quantities which characterize its shape: the diameter of , , and the radius, , of the dimensional ball inscribed in . The simplicial mesh is shape regular if and only if uniformly with respect to .
With the partitioning, , we associate a triple of piecewise polynomial, finitedimensional spaces,
(16) 
While we specify two choices of the space later, we fix and as follows,
where is the onedimensional space of constant functions on . We note that the inclusions listed in (16) imply that the elements of are continuous on , the functions in have continuous normal components across element boundaries, and that the functions in are in . This choice of is the standard lowest order RaviartThomasNédélec space (RT0) and is the piecewise constant space (P0).
2.2 Effects of permeability on the error of approximation
For , we start with a popular finiteelement approximation for (12)–(14) by choosing
where is the space of linear polynomials on . Then, is the space of piecewise linear (with respect to ), continuous vectorvalued functions. For uniformly positive definite permeability tensor, , such choice of spaces has been successfully employed for numerical simulations of Biot’s consolidation model (see (18); (23)). However, the heuristic considerations that expose some of the issues with this discretization are observed in cases when . In such cases, and the discrete problem approaches a P1P0 discretization of the Stokes’ equation. As it is well known, the element pair, , does not satisfy the infsup condition and is unstable for the Stokes’ problem. In fact, on a uniform grid in , it is easy to prove that volumetric locking occurs, namely, that the only divergencefree function from is the zero function. More precisely,
These inequalities imply that is not an onto operator, and, hence, the pair of spaces violates the infsup condition (24) associated with the discrete Stokes’ problem. More details on this undesirable phenomenon for Stokes are found in the classical monograph (24) and also in ((25), pp. 45–100) and (26).
Here, we demonstrate numerically that for Biot’s model, the error in the finiteelement approximation does not decrease when the permeability is small relative to the mesh size. We consider , and approximate (12)(14) subject to Dirichlet boundary conditions for both and on the whole of . We cover with a uniform triangular grid by dividing an uniform square mesh into right triangles. The material parameters are , , , , and . We consider a diagonal permeability tensor with constant , and the other data is set so that the exact solution is given by
Finally, we set and , so that we only perform one time step.
As seen in Table 2.1 the energy norm ( for ) for the displacement errors and the norm for pressure errors do not decrease until the mesh size is sufficiently small (compared with the permeability). Thus for small permeabilities, this could result in expensive discretizations which are less applicable to practical situations.
\diagbox  

3 Stabilization and perturbation of the bilinear form
To resolve the above issue, we introduce a wellknown stabilization technique based on enrichment of the piecewise linear continuous finiteelement space, , with edge/face (2D/3D) bubble functions (see ((27), pp. 145149)). The discretization described below is based on a Stokesstable pair of spaces with . As we show later, in Section 4, this stabilization gives a proper finiteelement approximation of the solution of Biot’s model independently of the size of the hydraulic conductivity, .
3.1 Stabilization by face bubbles
To define the enriched space, following (27), consider the set of dimensional faces from and denote this set by , where is the set of interior faces (shared by two elements) and is the set of faces on the boundary. In addition, is the set of faces on the boundary and . Note, if (pure traction boundary condition), then and . For any face , such that , and , let be the outward (with respect to ) unit normal vector to . With every face , we also associate a unit vector which is orthogonal to it. Clearly, if we have . For the boundary faces , we always set , where is the unique element for which we have . For the interior faces, the particular direction of is not important, although it is important that this direction is fixed. More precisely,
(20) 
Further, with every face , , we associate a vectorvalued function ,
(21) 
where , are barycentric coordinates on and is the vertex opposite to the face in . We note that is a continuous piecewise polynomial function of degree .
Finally, the stabilized finiteelement space is defined as
(22) 
The degrees of freedom associated with are the values at the vertices of and the total flux through of , where is the standard piecewise linear interpolant, . Then, the canonical interpolant, , is defined as:
With this choice of , the variational form, (17)–(19), remains the same and we have the following block form of the discrete problem:
(23) 
where , , and are the unknown vectors for the bubble components of the displacement, the piecewise linear components of the displacement, the Darcy velocity, and the pressure, respectively. The blocks in the definition of correspond to the following bilinear forms:
where , , , and an analogous decomposition for .
Next, we define the following notion of stability for discretizations of Biot’s model needed for the analysis.
Definition 3.1.
The triple of spaces is StokesBiot stable if and only if the following conditions are satisfied:

, for all , ;

, for all ;

The pair of spaces is Poisson stable, i.e., it satisfies stability and continuity conditions required by the mixed discretization of the Poisson equation;

The pair of spaces is Stokes stable.
Here, and denote the standard norm and norm, respectively.
To show how stability for the Biot’s system follows from the conditions above, we introduce a norm on :
(24) 
Further, we associate a composite bilinear form on the space, ,
We then have the following theorem which shows that on every time step the discrete problem is solvable.
Theorem 3.2.
If the triple is StokesBiot stable, then:

is continuous with respect to ; and

the following infsup condition holds.
(25) with a constant independent of mesh size and time step size .
Proof.
Note that if we replace with any spectrally equivalent bilinear form on , the same stability result holds true. In the next section, we introduce such a spectrally equivalent bilinear form which allows for: (1) Efficient elimination of the degrees of freedom corresponding to the bubble functions via static condensation; and (2) Derivation of optimal error estimates for the fully discrete problem, following the analysis in (28).
4 Local perturbation of the bilinear form and elimination of bubbles
In this section, we show how the unknowns corresponding to degrees of freedom in can be eliminated. A straightforward elimination of the edge/face bubbles is not local, and, in general, leads to prohibitively large number of nonzeroes in the resulting linear system. To resolve this, we introduce a consistent perturbation of , which has a diagonal matrix representation. It is then easy to eliminate the unknowns corresponding to the bubble functions in with no fillin. This leads to a stable P1RT0P0 discretization for the Biot’s model and, consequently, to a stable P1P0 discretization for the Stokes’ equation.
First, consider a natural decomposition of :
(26) 
and the local bilinear forms for , , and :
(27) 
For the restriction of onto the space spanned by bubble functions , we have
On each element, , then introduce
(28) 
Replacing with gives a perturbation, , of :
(29) 
4.1 A spectral equivalence result
To prove that the form and are spectrally equivalent, we need several auxiliary results. First, recall the definition of the rigid body motions (modes), on :
where is the algebra of skewsymmetric matrices. The dimension of is and its elements are componentwise linear vectorvalued functions.
Next, recall the classical Korn inequality (29); (30) for for a domain , starshaped with respect to a ball. As shown by Kondratiev and Oleinik in (31); (32),
(30) 
where the constant hidden in depends on the shape regularity of , that is, on the ratio . For convenience when referencing (30) later, we state the following lemma, which gives a simpler version of the inequality defined on simplices, where .
Lemma 4.1.
Let be a shaperegular simplicial mesh covering . Then, the following inequality holds for any and :
(31) 
where the constant hidden in “” depends on the shape regularity constant of .
Defining the unscaled bilinear form, ,
(32) 
we have the following local, spectral equivalence result.
Lemma 4.2.
For all the following inequalities hold:
(33) 
where the constant is independent of and .
Proof.
Set and note that for all . The upper bound follows immediately by several (two) applications of the CauchySchwarz inequality:
We prove the lower bound by establishing the following inequalities for :
(34) 
By definition for all and all rigid body modes , we have that and . The classical interpolation estimates found in ((33), Chapter 3) give
Taking the infimum over all and applying Korn’s inequality (Lemma 4.1) then yields
This shows the first inequality in (34), and to prove the second inequality, we note that from the definition of and the inverse inequality, we have that
Recalling the definition of in (21) and the formula for integrating powers of the barycentric coordinates, gives
(35) 
It follows that and , with . As the Gram matrix is positive semidefinite,
Multiplying by on both sides of this inequality furnishes the proof of (34), completing the proof of the Lemma. ∎
Next, we show the spectral equivalence for the bilinear forms and .
Lemma 4.3.
The following inequalities hold:
where depends on the shape regularity of the mesh.
Proof.
Let , . From the definition of in (28), , and the lower bound follows immediately:
To prove the upper bound, we use the following local estimate, which is established using an inverse inequality, a standard interpolation estimate, and for all rigid body modes ,
Taking the infimum over all and applying the Korn’s inequality (Lemma 4.1) then yields