New stabilized discretizations for poroelasticity and the Stokes’ equations

# New stabilized discretizations for poroelasticity and the Stokes’ equations

## Abstract

In this work, we consider the popular P1-RT0-P0 discretization of the three-field formulation of Biot’s consolidation problem. Since this finite-element formulation does not satisfy an inf-sup 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 finite-element 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 P1-RT0-P0 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 P1-P0 discretization). Numerical tests confirm the theory for both poroelastic and Stokes’ test problems.

###### keywords:
Stable finite elements, poroelasticity, Stokes’ equations

## 1 Introduction

The interaction between the deformation and fluid flow in a fluid-saturated porous medium is the object of study in poroelasticity theory. Such coupling has been modelled in the early one-dimensional work of Terzaghi (1). A more general three-dimensional 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 finite-difference schemes (6); (7) and finite-volume methods (see (8) for recent developments). Finite-element 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 finite-element schemes are constructed by either choosing discrete spaces satisfying appropriate inf-sup (or LBB) conditions, or applying suitable stabilization techniques to unstable finite-element pairs. For the two-field (displacements-pressure) formulation of Biot’s problem, the classical Taylor-Hood 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 three-field formulations (including the Darcy velocity), a stable finite-element method based on non-conforming Crouzeix-Raviart finite elements for the displacements, lowest order Raviart-Thomas-Nédélec elements for the Darcy velocity, and piecewise constants for the pressure, was proposed in (14). For a four-field 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 three-field 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 three-field formulation based on piecewise linear elements for displacements, Raviart-Thomas-Nédélec elements for the fluid flux, and piecewise constants for the pressure. This element, however, does not satisfy an inf-sup 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 well-known unstable pair based on piecewise linear and piecewise constant elements for velocities and pressure, respectively (P1-P0). This yields a stable finite-element 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 three-field 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 finite-element 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 well-posedness of the resulting scheme, as well as the corresponding error analysis are also provided here. In Section 5, we present the Stokes-stable finite-element method based on P1-P0 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 quasi-static 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: −divσ′+α∇p=ρg,inΩ, (1) constitutive equation: σ′=2με(u)+λdiv(u)I,inΩ, (2) compatibility condition: ε(u)=12(∇u+∇ut),inΩ, (3) Darcy’s law: w=−1μfK(∇p−ρfg),inΩ, (4) continuity equation: ∂∂t(1Mp+α divu)+divw=f,inΩ, (5)

where and are the Lamé coefficients, is the Biot modulus, and is the Biot-Willis 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 vector-valued 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 so-called three-field 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:

 −divσ′+α∇p=ρg,whereσ′=2με(u)+λdiv(u)I, (6) K−1μfw+∇p=ρfg, (7) ∂∂t(1Mp+α divu)+divw=f. (8)

This system is often subject to the following set of boundary conditions:

 p=0,forx∈¯¯¯¯Γt,σ′n=0,forx∈Γt, (9) u=0,forx∈¯¯¯¯Γc,w⋅n=0,forx∈Γc, (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,

 (1Mp+αdivu)(x,0)=0,x∈Ω, (11)

which yields the following mixed formulation of Biot’s three-field consolidation model:
For each , find such that

 a(u,v)−(αp,divv)=(ρg,v),∀ v∈V, (12) (K−1μfw,r)−(p,divr)=(ρfg,r),∀ r∈W, (13) (1M∂p∂t,q)+(αdiv∂u∂t,q)+(divw,q)=(f,q),∀ q∈Q, (14)

where,

 a(u,v)=2μ∫Ωε(u):ε(v)+λ∫Ωdivudivv, (15)

corresponds to linear elasticity. The function spaces used in the variational form are

 V={u∈H1(Ω) | u|¯¯¯Γc=0}, W={w∈H(div,Ω) | (w⋅n)|Γc=0}, Q=L2(Ω),

where is the space of square integrable vector-valued functions whose first derivatives are also square integrable, and contains the square integrable vector-valued functions with square integrable divergence.

We recall that the well-posedness of the continuous problem was established by Showalter (4), and, for the three-field 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, finite-dimensional spaces,

 Vh⊂V,Wh⊂W,Qh⊂Q. (16)

While we specify two choices of the space later, we fix and as follows,

 Wh={wh∈W | wh|T=a+ηx, a∈Rd, η∈R,  ∀T∈Th}, Qh={qh∈Q | qh|T∈P0(T),  ∀T∈Th},

where is the one-dimensional 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 Raviart-Thomas-Nédélec space (RT0) and is the piecewise constant space (P0).

Finally, using backward Euler as a time discretization on a time interval with constant time-step size , the discrete scheme corresponding to the three-field formulation (12)-(14) reads:
Find such that

 a(umh,vh)−(αpmh,divvh)=(ρg,vh),∀ vh∈Vh, (17) τ(K−1μfwmh,rh)−τ(pmh,divrh)=τ(ρfg,rh),∀ rh∈Wh, (18) (1Mpmh,qh)+(αdivumh,qh)+τ(divwmh,qh)=(˜f,qh),∀ qh∈Qh, (19)

where , and,

 (umh,wmh,pmh)≈(u(⋅,tm),w(⋅,tm),p(⋅,tm)),tm=mτ, m=1,2,…

### 2.2 Effects of permeability on the error of approximation

For , we start with a popular finite-element approximation for (12)–(14) by choosing

 Vh=Vh,1,withVh,1 := Extra open brace or missing close brace

where is the space of linear polynomials on . Then, is the space of piecewise linear (with respect to ), continuous vector-valued 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 P1-P0 discretization of the Stokes’ equation. As it is well known, the element pair, , does not satisfy the inf-sup 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 divergence-free function from is the zero function. More precisely,

 dim(Qh)>dimVh>dimRange(divh),divh=div∣∣Vh.

These inequalities imply that is not an onto operator, and, hence, the pair of spaces violates the inf-sup 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 finite-element 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

 u(x,y,t) = curlφ=(∂yφ−∂xφ),φ(x,y)=[xy(1−x)(1−y)]2, p(x,y,t) = 1.

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.

## 3 Stabilization and perturbation of the bilinear form

To resolve the above issue, we introduce a well-known stabilization technique based on enrichment of the piecewise linear continuous finite-element space, , with edge/face (2D/3D) bubble functions (see ((27), pp. 145-149)). The discretization described below is based on a Stokes-stable pair of spaces with . As we show later, in Section 4, this stabilization gives a proper finite-element 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,

 ne=ne,T+=−ne,T−ife=T+∩T−,andT±∈Th, (20)

Further, with every face , , we associate a vector-valued function ,

 Φe=φene,withφe∣∣∣T±=φe,T±,andφe,T±=d+1∏k=1,k≠j±λk,T±, (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 finite-element space is defined as

 Vh=Vh,1⊕Vb,Vb=span{Φe}e∈Eo,t. (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:

 Πv=Π1v+∑e∈Eo,tveΦe,ve=1|e|∫e(I−Π1)v.

With this choice of , the variational form, (17)–(19), remains the same and we have the following block form of the discrete problem:

 A⎛⎜ ⎜ ⎜⎝UbUlWP⎞⎟ ⎟ ⎟⎠=b,  with  A=⎛⎜ ⎜ ⎜ ⎜⎝AbbAbl0GbATblAll0Gl00τMwτGGTbGTlτGT−Mp⎞⎟ ⎟ ⎟ ⎟⎠, (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:

 a(ubh,vbh)→Abb,a(ulh,vbh)→Abl,a(ulh,vlh)→All, −(αph,divvbh)→Gb,−(αph,divvlh)→Gl,−(ph,divrh)→G, (K−1μfwh,rh)→Mw,(1Mph,qh)→Mp,

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 Stokes-Biot 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 :

 ∥(uh,wh,ph)∥τ=(∥uh∥21+τ∥wh∥2+τ2∥divwh∥2+∥ph∥2)1/2, (24)

Further, we associate a composite bilinear form on the space, ,

 B(uh,wh,ph;vh,rh,qh) := a(uh,vh)−(αph,divvh)+τ(K−1μfwh,rh)−τ(ph,divrh) −(1Mph,qh)−(αdivuh,qh)−τ(divwh,qh).

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 Stokes-Biot stable, then:

• is continuous with respect to ; and

• the following inf-sup condition holds.

 sup(vh,rh,qh)∈Vh×Wh×QhB(uh,wh,ph;vh,rh,qh)∥(vh,rh,qh)∥τ≥γ∥(uh,wh,ph)∥τ, (25)

with a constant independent of mesh size and time step size .

###### Proof.

Using the conditions in Definition 3.1, the two items in the statement of this theorem follow from arguments identical to the ones given in the proof of Theorem 2 in (28). ∎

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 non-zeroes 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 fill-in. This leads to a stable P1-RT0-P0 discretization for the Biot’s model and, consequently, to a stable P1-P0 discretization for the Stokes’ equation.

First, consider a natural decomposition of :

 u=ul+ub=Π1uul+∑e∈Eo,tueΦeub, (26)

and the local bilinear forms for , , and :

 aT(u,v)=2μ∫Tε(u):ε(v)+λ∫Tdivudivv. (27)

For the restriction of onto the space spanned by bubble functions , we have

 ab(ub,vb):=a(ub,vb)=∑T∈Thab,T(ub,vb)=∑T∈Th∑e,e′∈∂Tueve′aT(Φe′,Φe).

On each element, , then introduce

 db,T(u,v)=(d+1)∑e∈∂TueveaT(Φe,Φe),db(u,v)=∑T∈Thdb,T(u,v). (28)

Replacing with gives a perturbation, , of :

### 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 :

 R={v=a+bx∣∣a∈Rd,b∈so(d)},

where is the algebra of skew-symmetric matrices. The dimension of is and its elements are component-wise linear vector-valued functions.

Next, recall the classical Korn inequality (29); (30) for for a domain , star-shaped with respect to a ball. As shown by Kondratiev and Oleinik in (31); (32),

 infm∈so(d)∥∇u−m∥L2(Y)≲∥ε(u)∥L2(Y), (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 shape-regular simplicial mesh covering . Then, the following inequality holds for any and :

 infm∈so(n)∥∇u−m∥L2(T)≲∥ε(u)∥L2(T), (31)

where the constant hidden in “” depends on the shape regularity constant of .

Defining the unscaled bilinear form, ,

 ˜db,T(u,v):=∑e∈∂TueveaT(Φe,Φe), (32)

we have the following local, spectral equivalence result.

###### Lemma 4.2.

For all the following inequalities hold:

 ηT˜db,T(u,u)≤ab,T(u,u)≤(d+1)˜db,T(u,u),for allu∈Vb, (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 Cauchy-Schwarz inequality:

 ab,T(u,u) = ∑e,e′∈∂Taee′ueue′≤∑e,e′∈∂T√aeeae′e′|ueue′|=(∑e∈∂T√aee|ue|)2 ≤ (d+1)∑e∈∂Taeeu2e=(d+1)˜db,T(u,u).

We prove the lower bound by establishing the following inequalities for :

 h−2T∥u∥2L2(T)≲ab,T(u,u),and˜db,T(u,u)≲h−2T∥u∥2L2(T). (34)

By definition for all and all rigid body modes , we have that and . The classical interpolation estimates found in ((33), Chapter 3) give

 ∥u∥2L2(T) = ∥u−r−Π1(u−r)∥2L2(T)≲h2T∥∇(u−r)∥2L2(T).

Taking the infimum over all and applying Korn’s inequality (Lemma 4.1) then yields

 h−2T∥u∥2L2(T) ≲ infr∈R∥∇(u−r)∥2L2(T)=infm∈so(d)∥∇u−m∥2L2(T)≲∥ε(u)∥2L2(T).

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

 ˜db,T(u,u)≲∑e∈∂Tu2e[∥∇Φe∥2L2(T)+λ∥divΦe∥2L2(T)]≲h−2T∑e∈∂Tu2e∥Φe∥2L2(T).

Recalling the definition of in (21) and the formula for integrating powers of the barycentric coordinates, gives

 Φe=φene,∫Tλβ11…λβd+1d+1dx=|T|β1!…βd+1!d!(β1+…+βd+1+d)!. (35)

It follows that and , with . As the Gram matrix is positive semi-definite,

 ∑e∈∂Tu2e∥Φe∥2L2(T) = cd|T|∑e∈∂Tu2e≤cd|T|⎡⎣∑e∈∂Tu2e+∑e,e′∈∂Tueue′(ne⋅ne′)⎤⎦ = ∥∥ ∥∥∑e∈∂TueΦe∥∥ ∥∥2L2(T)=∥u∥2L2(T).

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: