A comparison of two nonconforming finite element methods for linear three-field poroelasticity

Jun Guo; Yanchao Shi; Weihua Luo

doi:10.1515/dema-2024-0073

Artikel Open Access

A comparison of two nonconforming finite element methods for linear three-field poroelasticity

Jun Guo , Yanchao Shi und Weihua Luo

Veröffentlicht/Copyright: 11. Dezember 2024

Veröffentlicht von

Veröffentlichen auch Sie bei De Gruyter Brill

Informationen für Autor*innen Erkunden Sie dieses Fachgebiet

Aus der Zeitschrift Demonstratio Mathematica Band 57 Heft 1

Abstract

We present and analyze two kinds of nonconforming finite element methods for three-field Biot’s consolidation model in poroelasticity. We employ the Crouzeix-Raviart element for one of the displacement component and conforming linear element for the remaining component, the lowest order Raviart-Thomas element (or the first-order Brezzi-Douglas-Marini element) for the fluid flux, and the piecewise constant for the pressure. We provide the corresponding analysis, including the well-posedness and a priori error estimates, for the fully discrete scheme coupled with the backward Euler finite difference for the time discretization. Such scheme ensures that the discrete Korn’s inequality is satisfied without adding any stabilization terms. In particular, it is free of poroelasticity locking. Numerical results are presented to compare the accuracy and locking-free performance of the two introduced schemes.

Keywords: Biot’s consolidation model; nonconforming finite element method; poroelasticity locking; locking-free

MSC 2010: 65N30; 74S05

1 Introduction

Biot’s consolidation model in poroelasticity is pioneered by Terzaghi [1] and further promoted by Biot [2], which allows one to describe the fluid flow in a fluid-saturated porous skeleton. As a classical model in the subsurface consolidation processes, it still has applications today in solid waste disposal, food processing, geological carbon sequestration, perfusion of bones and soft living tissues, and so on; see the study by Anandarajah [3] for the details. Therefore, the research for mechanical properties of such model is of great practical significance.

Until now, there have been tremendous amount of work for Biot’s consolidation models [4–13]. However, since the parameters involved in such models can vary over many orders of magnitude [14], the success in accurately obtaining poroelasticity solutions is often blocked by so-called poroelasticity locking, i.e., volumetric (Poisson) locking when porous skeleton is nearly incompressible (the Lamé coefficient λ tends to infinity) [15], as well as nonphysical pressure oscillations (NPOs) for the case, where c 0 = 0 and the hydraulic conductivity of porous skeleton is very low [16].

Also, there have been different explanations for the cause of poroelasticity locking. Initially, NPOs were attributed to the violation of the discrete inf-sup (compatibility) condition by Murad and Loula [17]. Subsequently, Rodrigo et al. [16] gave the relevant numerical examples to illustrate that the discrete inf-sup condition was not sufficient. On the other hand, Phillips and Wheeler [18] concluded that the poroelasticity locking occurred in the initial deformation of porous solid with a low permeability, when the solid medium exhibited incompressibility ( λ → ∞ ), the fluid became nearly incompressible ( c 0 = 0 ), and/or a small time step was employed. Besides, from the perspective of the numerical algebra, Yi [15] found that NPOs were due to the incompatibility of the discrete displacement and pressure space if the discrete fluid variables space satisfied the compatibility condition, and that volumetric locking was similarly to locking phenomena in linear elasticity problem [19]. We also refer to [20,21] and the references therein for details.

Correspondingly, various finite element methods (FEMs) with different formulations have recently been proposed to alleviate poroelasticity locking. For the two-field formulation: by adding a perturbation to the mass balance equation, NPOs were removed in previous studies [16,22]; based on a hybrid high-order discretization for displacement and a symmetric weighted interior penalty discretization for the pressure, a nonconforming high-order method without NPOs was put forward on general meshes in Boffi et al. [4]; using the super-convergent HDG discretization for the pressure Poisson operator and the divergence-conforming HDG method for the elasticity operator, poroelasticity locking has also been avoided in the study by Fu [23]. For the three-field formulation: combining the Bernardi-Raugel element for displacements and the mixed FEM for fluid variables, a family of locking-free FEMs were proposed in [15]; making the mass lumping for the Raviart-Thomas element, a nonconforming FEM was also addressed to overcome NPOs in [24]. For the four-field formulation: based on the Hellinger-Reissner formulation for the mechanical subproblem and the standard mixed FEM for the fluid subproblem, a remedy for poroelasticity locking was discussed in [25]; poroelasticity locking could also no longer occur in [26], because it was robust for all parameters. Besides, for the discussion of poroelasticity locking with random coefficients (uncertain inputs), we also refer to [27,28].

The purpose of this article is to develop a locking-free nonconforming FEM for three-field Biot’s consolidation model in poroelasticity. To this end, we need to solve the following problems:

under the background of the nonconforming FEM, the discrete Korn’s inequality must be readily satisfied to ensure the stability of the proposed scheme [4,29];
to avoid NPOs, the discrete inf-sup condition of the displacement and pressure space must be satisfied, except ensuring the compatibility of the fluid variables [15,30];
to overcome volumetric locking, the discrete divergence operator must be carefully designed to ensure a suitable commutativity [29,31], so that the error analysis of the proposed scheme is independent of λ .

For such reasons, we shall employ the Crouzeix-Raviart element for one of the displacement component and conforming linear element for the remaining component, the lowest-order Raviart-Thomas element (RT0) for the fluid flux, and the piecewise constant (P0) for the pressure. The corresponding analysis, including the well-posedness and a priori error estimates, is provided for the fully discrete scheme coupled with the backward Euler finite difference for the time discretization. Different from the scheme proposed by Hu et al. [24], our scheme not only ensures that the discrete Korn’s inequality is satisfied without adding any stabilization terms but also is free of poroelasticity locking. At the same time, it also reduces the degree of freedom. In particular, we also use the first-order Brezzi-Douglas-Marini element (BDM1) for the fluid flux, but only sub-optimal L 2 -error estimations can be obtained in the theoretical analysis, which is why RT-element is popularly selected in most of the literatures [12,15,24,30]. However, through numerical experiments, a interesting discovery is that optimal L 2 -error estimations of the fluid flux can still be obtained in this case, and it also keeps the locking-free performance.

The rest of this article is organized as follows: recalling a classical three-field variational formulation with displacements, pressure, and fluid flux as its unknowns, Biot’s consolidation model is sketched in Section 2. In Section 3, we make a weak constraint on the grid and then introduce the discrete space of the unknowns. The fully discrete nonconforming finite element scheme is presented, and then its well-posedness is proved in Section 4. The emphasis in Section 5 is set on giving a priori error estimates. Numerical results are provided in Section 6. Section 7 concludes this article.

2 Biot’s consolidation model

Let Ω ⊂ R d , d = 2 , 3 be a convex, polytopal, and bounded Lipschitz domain. For a given finite time J > 0 , external body force f , and volumetric source/sink g , we consider the following quasi-static Biot’s consolidation problem:

(1a) − ( λ + μ ) ∇ ( ∇ ⋅ u ) − μ ∇ 2 u + α ∇ p = f , in Ω × ( 0 , J ] ,

(1b) ∂ ∂ t ( c 0 p + α ∇ ⋅ u ) − ∇ ⋅ ( K ∇ p ) = g , in Ω × ( 0 , J ] ,

where the primary unknowns involve the displacement of porous skeleton u and the pressure of fluid p . The Lamé constants λ and μ are given by

(2) λ = E ν ( 1 + ν ) ( 1 − 2 ν ) and μ = E 2 ( 1 + ν ) ,

where ν is the Poisson’s ratio and E the Young’s modulus. The Biot-Willis constant α is usually equal to one in practice. The specific storage coefficient c 0 ( ≥ 0 ) vanishes for incompressible fluids. The permeability tensor K is symmetric and isotopic, satisfying

(3) k max − 1 ξ T ξ ≤ ξ T K − 1 ( x ) ξ ≤ k min − 1 ξ T ξ , ∀ ξ ∈ R d and ∀ x ∈ Ω ,

where 0 < k min ≤ k max < ∞ . Meanwhile, by introducing the fluid flux ϕ , the mass conservation of fluid obeys Darcy’s law:

(4) ϕ = − K ∇ p .

The effective stress tensor σ is defined by

(5) σ ( u ) ≔ 2 μ ε ( u ) + λ ∇ ⋅ u I , ε ( u ) ≔ 1 2 ( ∇ u + ∇ u T ) ,

where I is the identity matrix. We also denote the total stress tensor by σ ˜ that expresses the internal forces on surface element, and σ ˜ = σ − α p I .

In addition, we presuppose the boundary and initial conditions. We consider the homogeneous boundary conditions as following:

(6a) p = 0 , on Γ p D × ( 0 , J ] ,

(6b) ϕ ⋅ n = 0 , on Γ p N × ( 0 , J ] ,

(6c) u = 0 , on Γ u D × ( 0 , J ] ,

(6d) σ ˜ n = 0 , on Γ u N × ( 0 , J ] ,

where n is the unit outward normal of boundary ∂ Ω ; { Γ p D , Γ p N } and { Γ u D , Γ u N } are, respectively, a couple of partitions of ∂ Ω such that ∂ Ω ¯ = Γ ¯ p D ∪ Γ ¯ p N = Γ ¯ u D ∪ Γ ¯ u N , Γ p D ∩ Γ p N = ∅ , and Γ u D ∩ Γ u N = ∅ . We further assume Γ p D and Γ u D are closed, Γ p N and Γ u N are open, and

(7) Meas ( Γ u D ) > 0 , Meas ( Γ p D ) > 0 .

Also, we prescribe the following initial conditions:

(8) p ( 0 ) = p 0 and u ( 0 ) = u 0 , in Ω .

For the solvability of (1), one can refer to [32].

For sake of simplifying the notation, only homogeneous boundary conditions are considered here, but we stress that, e.g., nonhomogenous cases for the load vector σ ˜ n will be incorporated in numerical experiment, due to classical lifting arguments [31].

2.1 Variational formulation

The emphasis of this section is set on introducing a classical three-filed variational formulation for Biot’s consolidation model. Let’s start with some simplifications to make the discussion as simple as possible while preserving all the major difficulties, and thus set

(9) α = 1 .

By treating the fluid flux ϕ as an unknown of (1), then we have

(10a) − ( λ + μ ) ∇ ( ∇ ⋅ u ) − μ ∇ 2 u + ∇ p = f , in Ω × ( 0 , J ] ,

(10b) K − 1 ϕ + ∇ p = 0 , in Ω × ( 0 , J ] ,

(10c) c 0 p t + ∇ ⋅ u t + ∇ ⋅ ϕ = g , in Ω × ( 0 , J ] .

To proceed, we will adopt the classical symbols for Sobolev spaces. For example, given a measurable subset E in R d , for a negative integer s , we denote by H s ( E ) Sobolev spaces with the seminorm ∣ ⋅ ∣ s , E and norm ‖ ⋅ ‖ s , E , and ( ⋅ , ⋅ ) 0 , E and ⟨ ⋅ , ⋅ ⟩ 0 , ∂ E will denote the L 2 ( E ) and L 2 ( ∂ E ) inner product, respectively. The subscript E will be omitted, if E = Ω . Furthermore, we recall that

H ( div , Ω ) ≔ { φ ∈ ( L 2 ( Ω ) ) d ∣ ∇ ⋅ φ ∈ L 2 ( Ω ) } ,

equipped with the norm ‖ φ ‖ H ( div ) = ( ‖ φ ‖ 2 + ‖ ∇ ⋅ φ ‖ 2 ) 1 2 is also a Hilbert space. In particular, let X be a Banach space endowed with norm ‖ ⋅ ‖ X , and we denote by C 0 ( ( 0 , J ) ; X ) the set of functions v : ( 0 , J ) → X , which are continuous in t ∈ ( 0 , J ) . Also, we set

(11a) L ∞ ( ( 0 , J ) ; X ) ≔ { v ∈ X : esssup 0 ≤ t ≤ J ‖ v ‖ X 2 < ∞ } ,

(11b) L 2 ( ( 0 , J ) ; X ) ≔ v ∈ X : ∫ 0 J ‖ v ‖ X 2 d t < ∞ ,

(11c) C s ( ( 0 , J ) ; X ) ≔ v ∈ X : ∂ i v ∂ t i ∈ C 0 ( ( 0 , J ) ; X ) , 0 ≤ i ≤ s .

Let

V ≔ ( H 0 , Γ u D 1 ( Ω ) ) d = { v ∈ ( H 1 ( Ω ) ) d ∣ v ∣ u D = 0 } , W ≔ H 0 , Γ p N ( div , Ω ) = { φ ∈ H ( div , Ω ) ∣ φ ⋅ n ∣ Γ p N = 0 } , Q ≔ L 2 ( Ω ) ,

and set

( u , v ) ≔ 2 μ ( ε ( u ) : ε ( v ) ) + λ ( ∇ ⋅ u , ∇ ⋅ v ) , z ( ϕ , φ ) ≔ ( K − 1 ϕ , φ ) .

Then, we will concentrate the following three-filed variational formulation for the Biot’s problem: Find ( u , p ) ∈ C 1 ( ( 0 , J ] ; V × Q ) and ϕ ∈ C 0 ( ( 0 , J ] ; W ) such that

(12a) a ( u , v ) − ( ∇ ⋅ v , p ) = ( f , v ) , ∀ v ∈ V ,

(12b) z ( ϕ , φ ) − ( p , ∇ ⋅ φ ) = 0 , ∀ φ ∈ W ,

(12c) c 0 ( p t , q ) + ( ∇ ⋅ u t , q ) + ( ∇ ⋅ ϕ , q ) = ( g , q ) , ∀ q ∈ Q .

Obviously, a ( ⋅ , ⋅ ) is continuous and symmetric on V . Meanwhile, it follows from (7) that

(13) ‖ v ‖ 1 ≤ C ‖ ε ( v ) ‖ 0 , ∀ v ∈ V ,

which means that a ( ⋅ , ⋅ ) is coercive on V [31].

Here and below, C without or with subscripts will be used to represent a positive constant that is independent of the grid size h and λ , and its value may vary by location.

3 Mixed and nonconforming spaces

Let T h be a quasi-uniformity partition of Ω into the simplicial element for d = 2 and the parallelepiped element for d = 3 . We denote by h the maximum diameter of all element of Ω ¯ = ⋃ T ∈ T h T ¯ , i.e., h ≔ max T ∈ T h h T , by P k ( T ) the set of polynomials of degree ≤ k , and by Q k ( T ) the space of polynomials of degree ≤ k in each coordinate direction. Let Γ be the set of edges ( d = 2 ) or faces ( d = 3 ) of partition T h . For simplicity, we shall use the unified symbol e ∈ Γ to represent an edge ( d = 2 ) or face ( d = 3 ). We denote by ℰ h i = Γ \ ∂ Ω the set of interior edges (faces), and by ℰ h ∂ = Γ \ ℰ h i the set of boundary edges (faces).

Let e be an interior edge (face) shared by element T + and T − . Given a discontinuous vector function u , we denote jumps by

(14) ⟦ u ⟧ = u + ⋅ n + + u − ⋅ n − , on e ∈ ℰ h i and ⟦ u ⟧ = u ⋅ n , on e ∈ ℰ h ∂ ,

where the unit normal vectors n + and n − on e point exterior to T + and T − , respectively. Following [33], we also make a weak restriction on the partition T h (Figure 1), as follows:

Figure 1

The coarsest mesh with h = 1 ⁄ 4 .

Assumption 3.1

Every element of partition T h has at least one vertex in the interior of the domain.

3.1 Discretization of the fluid flux and pressure

We consider using RT0 and the piecewise constant to approximate the fluid flux and pressure, respectively. By imposing the boundary conditions, we have

(15a) W h ≔ { φ ∈ W ∣ φ ∣ T = ( P 0 ( T ) ) d + x P 0 ( T ) , ∀ T ∈ T h } ,

(15b) Q h ≔ { q ∈ Q ∣ q ∣ T ∈ P 0 ( T ) , ∀ T ∈ T h } .

Then, it follows from [34] that:

Lemma 3.1

There exists a positive constant C 1 such that

(16) sup φ h ∈ W h ( ∇ ⋅ φ h , p h ) ‖ φ h ‖ H ( div ) ≥ C 1 ‖ p h ‖ 0 , ∀ p h ∈ Q h .

Meanwhile, there exist the interpolation I h : H ( div ; Ω ) ∩ ( H 1 ( Ω ) ) d → W h and the projection P h : L 2 ( Ω ) ∩ H 1 ( Ω ) → Q h , such that

(17a) ( ∇ ⋅ ( φ − I h φ ) , q h ) = 0 , ∀ q h ∈ Q h ,

(17b) ∥ φ − I h φ ∥ 0 ≤ C h m ‖ φ ‖ m , 0 ≤ m ≤ 1 ,

(17c) ( p − P h p , q h ) = 0 , ∀ q h ∈ Q h ,

(17d) ∥ p − P h p ∥ 0 ≤ C h m ‖ p ‖ m , 0 ≤ m ≤ 1 ,

(17e) ∇ ⋅ I h φ = P h ∇ ⋅ φ .

Alternatively, W h is also chosen as BDM1, so that Lemma 3.1 and (17) hold [34].

3.2 Approximation of the displacement

This subsection will introduce the approximation space of the displacement in the case of d = 2 and d = 3 . For the case d = 2 , the nonconforming linear element is applied to approximate the displacement. That is, the Crouzeix-Raviart element for one of displacement components and linear conforming element for other components:

V 1 , h ≔ v ∈ L 2 ( Ω ) ∣ v ∣ T ∈ P 1 ( T ) , ∀ T ∈ T h and ∫ e ⟦ v ⟧ = 0 , ∀ e ∈ ℰ h i , v = 0 at the midpoints of the edges lying on Γ u D , V 2 , h ≔ { v ∈ H 0 , Γ u D 1 ( Ω ) ∣ v ∣ T ∈ P 1 ( T ) , ∀ T ∈ T h } .

For the case d = 3 , we define the conforming spaces

V ˜ 1 , h ≔ { v ∈ H 0 , Γ u D 1 ( Ω ) ∣ v ^ ∣ T ∈ Q 1 ( T ^ ) , ∀ T ∈ T h } , V ˜ 2 , h ≔ { v ∈ H 0 , Γ u D 1 ( Ω ) ∣ v ^ ∣ T ∈ Q 2 ( T ^ ) , ∀ T ∈ T h } ,

and the nonconforming space

V ˜ 3 , h ≔ v ∈ L 2 ( Ω ) ∣ v ∣ T ∈ R Q 1 ( T ) , ∀ T ∈ T h and ∣ e ∣ − 1 ∫ e [ [ v ] ] = 0 , ∀ e ∈ ℰ h i , ∣ e ∣ − 1 ∫ e v = 0 , ∀ e ∈ ℰ h ∂ ∩ Γ u D ,

where we denote by G T the tri-linear mapping from T ˆ to T [35], and where the local space is chosen as follows:

(18) R Q 1 ( T ) ≔ { v ^ ∘ G T − 1 ∣ v ^ ∈ span { 1 , x , y , z , x 2 − y 2 , x 2 − z 2 } } .

Therefore, the approximation space for displacement is defined as follows:

(19) V h ≔ ( V 1 , h , V 2 , h ) or V h ≔ ( V ˜ 1 , h , V ˜ 2 , h , V ˜ 3 , h ) .

Also, we define the following norms

(20) ∣ v ∣ 1 , h 2 ≔ ∑ T ∈ T h ∣ ∇ v ∣ 0 , T 2 and ‖ v ‖ 1 , h 2 ≔ ‖ v ‖ 0 2 + ∣ v ∣ 1 , h 2 , ∀ v ∈ V + V h ,

due to V h ⊈ V . Then, it follows from Lemma 4.3 of [33] and Theorem 3.1 of [35] that:

Lemma 3.2

Suppose that Assumption 3.1 is valid. Then there exists a positive constant C 2 , such that

(21) sup v h ∈ V h ∑ T ∈ T h ( ∇ ⋅ v h , p h ) 0 , T ‖ v h ‖ 1 , h ≥ C 2 ‖ p h ‖ 0 , ∀ p h ∈ Q h .

As a result, there exists the projection Π h : ( H 1 ( Ω ) ⋂ H 2 ( Ω ) ) d → V h [36], so that

(22a) ( ∇ ⋅ ( v − Π h v ) , q h ) T = 0 , ∀ q h ∈ Q h ,

(22b) ∥ v − Π h v ∥ s , T ≤ C h T m − s ‖ v ‖ m , T , 0 ≤ s < m ≤ 2 ,

(22c) ∇ ⋅ ( Π h v ∣ T ) = P h ∇ ⋅ v ∣ T .

Notice that V h can be combined in other forms, such as V h ≔ ( V 2 , h , V 1 , h ) for d = 2 , V h ≔ ( V ˜ 1 , h , V ˜ 3 , h , V ˜ 2 , h ) for d = 3 , and so on. For such reason, a comparison is also provided for the different V h in numerical experiments.

4 Fully discrete scheme

We shall devote this section to present the fully discrete scheme for the problem (12) at time t n ≔ n Δ t , n = 1 , 2 , … , N , as following: to find ( u h n , ϕ h n , p h n ) ∈ V h × W h × Q h such that

(23a) a h ( u h n , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p h n ) 0 , T = ( f n , v h ) , ∀ v h ∈ V h ,

(23b) z ( ϕ h n , φ h ) − ( p h n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h ,

(23c) c 0 p h n − p h n − 1 Δ t , q h + ∑ T ∈ T h ∇ ⋅ u h n − u h n − 1 Δ t , q h 0 , T + ( ∇ ⋅ ϕ h n , q h ) = ( g n , q h ) , ∀ q h ∈ Q h ,

where Δ t = J ⁄ N , and we have denoted s n ( x ) ≔ s ( x , t n ) , and where a h ( ⋅ , ⋅ ) is given by

(24) a h ( u , v ) ≔ 2 μ ∑ T ∈ T h ( ε ( u ) : ε ( v ) ) 0 , T + λ ∑ T ∈ T h ∇ ⋅ u , ∇ ⋅ ( v ) 0 , T , ∀ u , v ∈ V h .

Also, we choose the discrete initial solutions

(25) u h 0 = u ¯ h 0 and p h 0 = p ¯ h 0 ,

whose definition will be given in (32).

Remark 4.1

The serious drawback of our method lies in Assumption 3.1, so that the partition of Ω has certain limitation. Fortunately, under Assumption 3.1 the discrete Korn’s inequality

(26) ‖ v ‖ 1 , h 2 ≤ C ∑ T ∈ T h ( ε ( v ) : ε ( v ) ) 0 , T , ∀ v ∈ V h ,

has been established in Lemma 4.5 of [33] for d = 2 and Theorem 3.2 of [35] for d = 3 , so that the scheme (23) does therefore not need to include any parameterized stabilization term (e.g., jump terms [24,37]) that provides the stability.

4.1 Existence and uniqueness of solution

The goal of this subsection is, at each time step t n , to investigate the existence and uniqueness of solution to (23). To this end, we shall rewrite (23) as follows:

(27a) a h ( u h n , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p h n ) 0 , T = ( f n , v h ) , ∀ v h ∈ V h ,

(27b) Δ t z ( ϕ h n , φ h ) − Δ t ( p h n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h ,

(27c) c 0 ( p h n , q h ) + ∑ T ∈ T h ( ∇ ⋅ u h n , q h ) 0 , T + Δ t ( ∇ ⋅ ϕ h n , q h ) = ( g ˜ , q h ) , ∀ q h ∈ Q h ,

where (27b) and (27c) are obtained by multiplying (23b) and (23c) with Δ t , respectively, and we have set g ˜ ∣ T ≔ ( Δ t g n + c 0 p h n − 1 + ∇ ⋅ u h n − 1 ) ∣ T .

Theorem 4.1

The fully discrete scheme (23) admits, at each time step t n , a unique solution ( u h n , ϕ h n , p h n ) ∈ V h × W h × Q h .

Proof

It is trivial to see that the number of linear equations (27) equals the number of unknowns. Thus, it suffices to prove that corresponding homogeneous equations:

(28a) a h ( u h n , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p h n ) 0 , T = 0 , ∀ v h ∈ V h ,

(28b) Δ t z ( ϕ h n , φ h ) − Δ t ( p h n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h ,

(28c) c 0 ( p h n , q h ) + ∑ T ∈ T h ( ∇ ⋅ u h n , q h ) 0 , T + ( ∇ ⋅ ϕ h n , q h ) = 0 , ∀ q h ∈ Q h ,

have a unique set of zero solution at each time step t n . To that end, taking v h = u h n , φ h = ϕ h n and q h = p h n in above equations, and then adding them together, we have

(29) a h ( u h n , u h n ) + Δ t z ( ϕ h n , ϕ h n ) + c 0 ( p h n , p h n ) = 0 .

Next, we discuss the constrained specific storage coefficient in two cases: c 0 > 0 and c 0 = 0 . For the case of c 0 > 0 , it follows from (3) and (26) that

(30) ‖ u h n ‖ 1 , h 2 + ‖ ϕ h n ‖ 0 2 + ‖ p h n ‖ 0 2 = 0 ,

which, combining boundary conditions, means u h n = 0 , ϕ h n = 0 and p h n = 0 . For the case of c 0 = 0 , we still obtain u h n = 0 and ϕ h n = 0 , from which we derive that

(31) ∑ T ∈ T h ( ∇ ⋅ v h , p h n ) 0 , T = 0 , ∀ v h ∈ V h and ( p h n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h .

Using Lemmas 3.1 and 3.2, we can also obtain p h n = 0 .

To sum up, we end this proof.□

Notice that our proof covers both the case of an incompressible Newtonian fluid c 0 = 0 and the case of c 0 being positive.

5 Error estimates

The present section will be devoted to provide a priori error estimates. To this end, we start with the construction of an elliptic projection, just to make it easier for the error estimate (independent of λ ) of the scheme (23).

5.1 Elliptic projection

Given that ( u , ϕ , p ) ∈ C 0 ( ( 0 , J ) ; V × W × Q ) , we define the following elliptic projection: for t > 0 , to find ( u ¯ h , ϕ ¯ h , p ¯ h ) ∈ C 0 ( ( 0 , J ) ; V h × W h × Q h ) such that

(32a) a h ( u ¯ h , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p ¯ h ) 0 , T = a ˜ ( u , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p ) 0 , T − ∑ T ∈ T h ⟨ σ ˜ ⋅ n T , v h ⟩ 0 , ∂ T \ Γ u D , ∀ v h ∈ V h ,

(32b) z ( ϕ ¯ h , φ h ) − ( p ¯ h , ∇ ⋅ φ h ) = z ( ϕ , φ h ) − ( p , ∇ ⋅ φ h ) , ∀ φ h ∈ W h ,

(32c) ( ∇ ⋅ ϕ ¯ h , q h ) = ( ∇ ⋅ ϕ , q h ) , ∀ q h ∈ Q h ,

where a ˜ is given by

(33) a ˜ ( u , v h ) ≔ 2 μ ∑ T ∈ T h ( ε ( u ) : ε ( v h ) ) 0 , T + λ ∑ T ∈ T h ( ∇ ⋅ u , ∇ ⋅ v h ) 0 , T , ∀ v h ∈ V h .

Notice that, by decoupling the above projection, (3) and Lemma 3.1 imply that the mixed Poisson equations (32b) and (32c) have a unique solution ( ϕ ¯ h , p ¯ h ) ∈ W h × Q h , and then u ¯ h can be determined by (26) in the linear elasticity equation (32a). Therefore, (32) is well posed. Furthermore, we can establish the following estimates.

Lemma 5.1

Let ( u ¯ h , ϕ ¯ h , p ¯ h ) ∈ C 0 ( ( 0 , J ) ; V h × W h × Q h ) be the solution of elliptic projection (32). Suppose that ( u , ϕ , p ) ∈ C 0 ( ( 0 , J ) ; [ H 2 ( Ω ) ] d × [ H 1 ( Ω ) ] d × H 1 ( Ω ) ) and Assumption 3.1 is valid, then we have

(34a) ‖ u − u ¯ h ‖ 1 , h ≤ C h ( ‖ u ‖ 2 + ‖ λ ∇ ⋅ u ‖ 1 + ‖ ϕ ‖ 1 + ‖ p ‖ 1 ) ,

(34b) ‖ ϕ − ϕ ¯ h ‖ 0 ≤ C h ‖ ϕ ‖ 1 ,

(34c) ‖ p − p ¯ h ‖ 0 ≤ C h ( ‖ ϕ ‖ 1 + ‖ p ‖ 1 ) .

Proof

Let’s split, for t > 0 , errors as follows:

u − u ¯ h = u − Π h u + Π h u − u ¯ h , ϕ − ϕ ¯ h = ϕ − I h ϕ + I h ϕ − ϕ ¯ h , p − p ¯ h = p − P h p + P h p − p ¯ h .

Subtracting z ( I h ϕ , φ h ) and adding ( ∇ ⋅ φ h , P h p ) from both sides of (32b) leads to

(35) z ( ϕ ¯ h − I h ϕ , φ h ) − ( ∇ ⋅ φ h , p ¯ h − P h p ) = z ( ϕ − I h ϕ , φ h ) − ( ∇ ⋅ φ h , p − P h p ) .

On the other hand, it follows from (17e) and (32c) that

(36) ( ∇ ⋅ ϕ ¯ h , q h ) = ( ∇ ⋅ ϕ , q h ) = ( P h ∇ ⋅ ϕ , q h ) = ( ∇ ⋅ I h ϕ , q h ) , ∀ q h ∈ Q h ,

from which we can obtain

(37) ∇ ⋅ ( ϕ ¯ h − I h ϕ ) = 0 .

Then, taking φ h = ϕ ¯ h − I h ϕ in (35), we have

(38) z ( ϕ ¯ h − I h ϕ , ϕ ¯ h − I h ϕ ) = z ( ϕ − I h ϕ , ϕ ¯ h − I h ϕ ) ,

which means

(39) ‖ ϕ ¯ h − I h ϕ ‖ 0 ≤ C ‖ ϕ − I h ϕ ‖ 0 .

Thus, by (17b) and the triangle inequality, (34b) can be proved.

Next, according to Lemma 3.1, there exists φ ˜ h ∈ W h such that

(40) ∇ ⋅ φ ˜ h = − ( p ¯ h − P h p ) and ‖ φ ˜ h ‖ H ( div ) ≤ C ‖ p ¯ h − P h p ‖ 0 .

Similarly, setting φ h = φ ˜ h in (35) yields

(41) ( p ¯ h − P h p , p ¯ h − P h p ) = z ( ϕ − ϕ ¯ h , φ ˜ h ) + ( p ¯ h − P h p , p − P h p ) ,

which implies

(42) ‖ p ¯ h − P h p ‖ 0 ≤ C ( ‖ ϕ − ϕ ¯ h ‖ 0 + ‖ p − P h p ‖ 0 ) .

Hence, (34c) follows directly from (17d), (34b) and the triangle inequality.

Finally, subtracting a h ( Π h u , v h ) and adding ∑ T ∈ T h ( ∇ ⋅ v h , P h p ) 0 , T from both sides of (32a), we have

(43) a h ( u ¯ h − Π h u , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p ¯ h − P h p ) 0 , T = 2 μ ∑ T ∈ T h ( ε ( u − Π h u ) : ε ( v h ) ) 0 , T + λ ∑ T ∈ T h ( ∇ ⋅ u , ∇ ⋅ v h ) 0 , T − λ ∑ T ∈ T h ( P h ∇ ⋅ Π h u , ∇ ⋅ v h ) 0 , T − ∑ T ∈ T h ( ∇ ⋅ v h , p − P h p ) 0 , T − ∑ T ∈ T h ⟨ σ ˜ ⋅ n T , v h ⟩ 0 , ∂ T \ Γ u D .

Then, it follows from (22c) that

(44) a h ( u ¯ h − Π h u , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p ¯ h − P h p ) 0 , T = 2 μ ∑ T ∈ T h ( ε ( u − Π h u ) : ε ( v h ) ) 0 , T + λ ∑ T ∈ T h ( ∇ ⋅ u − P h ∇ ⋅ u , ∇ ⋅ v h ) 0 , T − ∑ T ∈ T h ( ∇ ⋅ v h , p − P h p ) 0 , T − ∑ T ∈ T h ⟨ σ ˜ ⋅ n T , v h ⟩ 0 , ∂ T \ Γ u D .

Meanwhile, for such a typical consistency error, we have

(45) ∑ T ∈ T h ⟨ σ ˜ ⋅ n T , v h ⟩ 0 , ∂ T \ Γ u D = ∑ e ∈ Γ \ Γ u D ⟨ σ ˜ − Π 0 0 ( σ ˜ ) , ⟦ v h ⟧ ⟩ 0 , e ≤ C h 1 2 ‖ σ ˜ ‖ 1 h 1 2 ‖ v h ‖ 1 , h .

where Π 0 0 is the L 2 ( e ) – projection onto the space of piecewise constant on e [38]. To proceed, taking v h = u ¯ h − Π h u in (44), and then employing (26) and the Cauchy-Schwartz inequality, we obtain

(46) ‖ u ¯ h − Π h u ‖ 1 , h ≤ C ( ‖ u − Π h u ‖ 1 , h + λ ‖ ∇ ⋅ u − P h ∇ ⋅ u ‖ 0 + ‖ p − p ¯ h ‖ 0 + h ‖ σ ˜ ‖ 1 ) .

Therefore, by combining (22b), (34c), and the triangle inequality, we complete the proof.□

Remark 5.1

Given that smooth time derivatives of u , ϕ , and p , we have the projection u ¯ h t , ϕ ¯ h t and p ¯ h t of u t , ϕ t and p t , respectively, so that

(47a) ‖ u t − u ¯ h t ‖ 1 , h ≤ C h ( ‖ u t ‖ 2 + ‖ λ ∇ ⋅ u t ‖ 1 + ‖ ϕ t ‖ 1 + ‖ p t ‖ 1 ) ,

(47b) ‖ ϕ t − ϕ ¯ h t ‖ 0 ≤ C h ‖ ϕ t ‖ 1 ,

(47c) ‖ p t − p ¯ h t ‖ 0 ≤ C h ( ‖ ϕ t ‖ 1 + ‖ p t ‖ 1 ) .

More generally, we also define the projection u ¯ h t t , ϕ ¯ h t t and p ¯ h t t of u t t , ϕ t t , and p t t , respectively.

Remark 5.2

As stated in [15, Theorem 3.3], sup 0 ≤ t ≤ J ‖ λ ∇ ⋅ u ‖ 1 and sup 0 ≤ t ≤ J ‖ λ ∇ ⋅ ∂ t u ‖ 1 can be bounded by suitable norms of f and g , and initial conditions of u and p . Therefore, the following error analysis will be independent of λ .

5.2 A priori error estimates

For convenience, for t > 0 , we set

(48) η u ≔ u − u ¯ h , η ϕ ≔ ϕ − ϕ ¯ h , η p ≔ p − p ¯ h , ξ u ≔ u ¯ h − u h , ξ ϕ ≔ ϕ ¯ h − ϕ h , ξ p ≔ p ¯ h − p h ,

giving

u − u h = η u + ξ u , ϕ − ϕ h = η ϕ + ξ ϕ , p − p h = η p + ξ p .

Now, we have the following a priori error estimates.

Theorem 5.2

Let ( u , p ) ∈ C 1 ( ( 0 , J ) ; V × Q ) and ϕ ∈ C 0 ( ( 0 , J ) ; W ) be the solution of (12) and ( u h n , ϕ h n , p h n ) ∈ V h × W h × Q h be the solution of (23). Given that the regularity assumptions

u t ∈ L 2 ( ( 0 , J ) ; ( H 2 ( Ω ) ) d ) , u t t ∈ L 2 ( ( 0 , J ) ; ( H 1 ( Ω ) ) d ) , p t , λ ∇ ⋅ u t , ϕ t ∈ L 2 ( ( 0 , J ) ; H 1 ( Ω ) ) , p t t ∈ L 2 ( ( 0 , J ) ; L 2 ( Ω ) ) ,

then we have

max 1 ≤ n ≤ N ∥ u n − u h n ∥ 1 , h 2 + Δ t ∑ n = 1 N ( ‖ p n − p h n ‖ 0 2 + ‖ ϕ n − ϕ h n ‖ 0 2 ) ≤ C ( h 2 + ( Δ t ) 2 ) .

Proof

First, we consider the case of c 0 = 0 . Note that the following Taylor expansion

(49) u n − u n − 1 Δ t = u t n + 1 Δ t ∫ t n − 1 t n ( s − t n − 1 ) u t t ( s ) d s ,

so, at time step t = t n , (12) can been rewritten as follows:

(50a) a ˜ ( u n , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , p n ) 0 , T = ( f n , v h ) + ∑ T ∈ T h ⟨ σ ˜ n ⋅ n T , v h ⟩ 0 , ∂ T \ Γ u , D , ∀ v h ∈ V + V h ,

(50b) z ( ϕ n , φ h ) − ( p n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h ,

(50c) ∑ T ∈ T h ∇ ⋅ u n − u n − 1 Δ t , q h 0 , T + ( ∇ ⋅ ϕ n , q h ) = ( g n , q h ) + 1 Δ t ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ u t t ( s ) d s , q h , ∀ q h ∈ Q h .

Then, by combining (23), (32), (50), and the notation of (48), we obtain

(51a) a h ( ξ u n , v h ) − ∑ T ∈ T h ( ∇ ⋅ v h , ξ p n ) 0 , T = 0 , ∀ v h ∈ V + V h ,

(51b) z ( ξ ϕ n , φ h ) − ( ξ p n , ∇ ⋅ φ h ) = 0 , ∀ φ h ∈ W h ,

(51c) ∑ T ∈ T h ∇ ⋅ ξ u n − ξ u n − 1 Δ t , q h 0 , T + ( ∇ ⋅ ξ ϕ n , q h ) = − ∑ T ∈ T h ∇ ⋅ η u n − η u n − 1 Δ t , q h 0 , T + 1 Δ t ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ u t t ( s ) d s , q h , ∀ q h ∈ Q h .

To proceed, taking v h = ξ u n − ξ u n − 1 Δ t , φ h = ξ ϕ n , and q h = ξ p n in (51), respectively, summing these equations, and then multiplying both sides of the resulting by Δ t , we also have

(52) a h ( ξ u n , ξ u n − ξ u n − 1 ) + Δ t z ( ξ ϕ n , ξ ϕ n ) = − ∑ T ∈ T h ( ∇ ⋅ ( η u n − η u n − 1 ) , ξ p n ) 0 , T + ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ u t t ( s ) d s , ξ p n .

Let ‖ ξ u n ‖ a h 2 ≔ a h ( ξ u n , ξ u n ) . Since the fact that ξ u 0 = 0 (25) and

(53) a h ( ξ u n , ξ u n − ξ u n − 1 ) ≥ 1 2 ( ‖ ξ u n ‖ a h 2 − ‖ ξ u n − 1 ‖ a h 2 ) ,

summing from 1 to N for (52) yields

(54) 1 2 ‖ ξ u N ‖ a h 2 + k max − 1 Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 ≤ ∑ i = 1 2 ℳ i ,

where

ℳ 1 ≔ − ∑ n = 1 N ∑ T ∈ T h ( ∇ ⋅ ( η u n − η u n − 1 ) , ξ p n ) 0 , T , ℳ 2 ≔ ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ u t t ( s ) d s , ξ p n .

Now, we will bound the qualities ℳ 1 and ℳ 2 . Since the following Taylor expansion:

(55) η u n − η u n − 1 = Δ t η u t n + ∫ t n − 1 t n ( s − t n − 1 ) η u t t ( s ) d s ,

and

(56) ∑ n = 1 N ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ η u t t ( s ) d s 0 ≤ ( Δ t ) 3 2 ∫ 0 J ‖ ∇ ⋅ η u t t ( s ) ‖ 0 2 d s 1 2 ,

using the Cauchy-Schwarz and Young’s inequalities leads to

(57) ∣ ℳ 1 ∣ ≤ Δ t ∑ n = 1 N ∑ T ∈ T h ∣ ( ∇ ⋅ η u t n , ξ p n ) 0 , T ∣ + ∑ n = 1 N ∑ T ∈ T h ∫ t n − 1 t n ( s − t n − 1 ) ∇ ⋅ η u t t ( s ) d s , ξ p n 0 , T ≤ ε 1 2 Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 + C Δ t ∑ n = 1 N ‖ ∇ ⋅ η u t n ‖ 0 2 + ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ η u t t ( s ) ‖ 0 2 d s ,

where we have utilized Lemma 3.1 and (51b) to obtain

(58) ‖ ξ p n ‖ 0 ≤ C sup φ h ∈ W h ( ∇ ⋅ φ h , ξ p n ) ‖ φ h ‖ H ( div ) = C sup φ h ∈ W h z ( ξ ϕ n , φ h ) ‖ φ h ‖ H ( div ) ≤ C ‖ ξ ϕ n ‖ 0 .

In the same way, we also have

(59) ∣ ℳ 2 ∣ ≤ ε 1 2 Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 + C ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ u t t ( s ) ‖ 0 2 d s .

By combining the aforementioned estimates for ℳ 1 − ℳ 2 with (54), we obtain

(60) 1 2 ‖ ξ u N ‖ a h 2 + ( k max − 1 − ε 1 ) Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 ≤ C Δ t ∑ n = 1 N ‖ ∇ ⋅ η u t n ‖ 0 2 + ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ η u t t ( s ) ‖ 0 2 d s + ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ u t t ( s ) ‖ 0 2 d s ,

which, together with the suitable ε 1 and (26), implies that

(61) ‖ ξ u N ‖ 1 , h 2 + Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 ≤ C Δ t ∑ n = 1 N ‖ ∇ ⋅ η u t n ‖ 0 2 + ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ η u t t ( s ) ‖ 0 2 d s + ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ u t t ( s ) ‖ 0 2 d s .

At the same time, from Remark 5.1, it follows that

(62) ‖ ∇ ⋅ η u t n ‖ 0 ≤ ‖ η u t n ‖ 1 , h ≤ ‖ u t n − u ¯ h t n − ( u t ( s ) − u ¯ h t ( s ) ) ‖ 1 , h + ‖ u t ( s ) − u ¯ h t ( s ) ‖ 1 , h ≤ C ‖ u t n − u t ( s ) ‖ 1 + ‖ u t ( s ) − u ¯ h t ( s ) ‖ 1 , h , ∀ s ∈ [ t n − 1 , t n ] .

Since

(63) u t n − u t ( s ) = ∫ s t n u t t ( r ) d r ,

it holds

(64) ‖ u t n − u t ( s ) ‖ 1 2 ≤ Δ t ∫ t n − 1 t n ‖ u t t ( s ) ‖ 1 2 d s .

Using again Remark 5.1 yields

(65) ‖ u t ( s ) − u ¯ h t ( s ) ‖ 1 , h 2 ≤ C h 2 ( ‖ u t ( s ) ‖ 2 2 + ‖ λ ∇ ⋅ u t ( s ) ‖ 1 2 + ‖ ϕ t ( s ) ‖ 1 + ‖ p t ( s ) ‖ 1 2 ) ,

As a consequence, we have

(66) Δ t ∑ n = 1 N ‖ ∇ ⋅ η u t n ‖ 0 2 ≤ C ( Δ t ) 2 ∫ 0 J ‖ u t t ( s ) ‖ 1 2 d s + C h 2 ∫ 0 J ‖ u t ( s ) ‖ 2 2 d s + ∫ 0 J ( ‖ λ ∇ ⋅ u t ( s ) ‖ 1 2 + ‖ ϕ t ( s ) ‖ 1 + ‖ p t ( s ) ‖ 1 2 ) d s .

It is also trivial to find that

(67) ( Δ t ) 2 ∫ 0 J ‖ ∇ ⋅ η u t t ( s ) ‖ 0 2 d s ≤ C ( Δ t ) 2 ∫ 0 J ‖ u t t ( s ) ‖ 1 2 d s .

Therefore, it follows from (61), (66), and (67) that

(68) max 1 ≤ n ≤ N ‖ ξ u n ‖ 1 , h 2 + Δ t ∑ n = 1 N ‖ ξ ϕ n ‖ 0 2 ≤ C ( h 2 + ( Δ t ) 2 ) ,

which, together with (58), implies

(69) max 1 ≤ n ≤ N ‖ ξ u n ‖ 1 , h 2 + Δ t ∑ n = 1 N ( ‖ ξ ϕ n ‖ 0 2 + ‖ ξ p n ‖ 0 2 ) ≤ C ( h 2 + ( Δ t ) 2 ) .

On the other hand, for c 0 > 0 , we can easily obtain

(70) max 1 ≤ n ≤ N ( ‖ ξ u n ‖ 1 , h 2 + ‖ ξ p n ‖ 0 2 ) + Δ t ∑ n = 1 N ( ‖ ξ ϕ n ‖ 0 2 + ‖ ξ p n ‖ 0 2 ) ≤ C ( h 2 + ( Δ t ) 2 ) ,

which means that (69) holds.

In conclusion, by combining Lemma 5.1 with (69), and then using the triangle inequality, we complete this proof.□

Remark 5.3

One can find that only suboptimal L 2 -error estimation for the fluid flux can be obtained in the above proof, when using BDM1 for the fluid flux. However, through numerical experiments, we find that optimal L 2 -error estimations can still be obtained in this case, and it also keeps the locking-free performance.

6 Numerical tests

In this section, we shall provide numerical tests to confirm the theoretical results presented in Section 5, i.e., to assess the accuracy and locking-free performance of our method. All the computations are performed on uniform meshes that just satisfy Assumption 3.1; see Figure 1.

6.1 Convergence and alleviating the volumetric locking

Following [15], we consider, on the domain Ω = [ 0 , 1 ] × [ 0 , 1 ] , the following smooth solution:

u 1 ( x , y , t ) = e − t sin ( 2 π y ) ( − 1 + cos ( 2 π x ) ) + 1 μ + λ sin ( π x ) sin ( π y ) , u 2 ( x , y , t ) = e − t sin ( 2 π x ) ( 1 − cos ( 2 π y ) ) + 1 μ + λ sin ( π x ) sin ( π y ) , p ( x , y , t ) = e − t sin ( π x ) sin ( π y ) ,

from which the corresponding body force f and the source/sink term g can be readily obtained from (10). One can find that, at any time t , the solution satisfies

∇ ⋅ u = π e − t sin ( π ( x + y ) ) ⁄ ( μ + λ ) → 0 , as λ → ∞ ,

which is designed to emphasize that our method is free of volumetric locking for a nearly incompressible porous media ( λ → ∞ ). To this end, we choose three different λ = 1 , 1 × 1 0 4 , 1 × 1 0 8 and parameters c 0 = 0 , α = 1 , K = 1 , μ = 1 , and define the following error quantities:

‖ u − u h ‖ L ∞ ( H 1 ) ≔ max 1 ≤ n ≤ N ∥ u n − u h n ∥ 1 , h , ‖ ϕ − ϕ h ‖ L 2 ( L 2 ) ≔ Δ t ∑ n = 1 N ‖ ϕ n − ϕ h n ‖ 0 , ‖ p − p h ‖ L 2 ( L 2 ) ≔ Δ t ∑ n = 1 N ‖ p n − p h n ‖ 0 .

In Tables 1–6, we list the numerical results, including the errors and convergence rates for displacements, fluid flux, and pressure. We observe from Table 1 that the errors ‖ u − u h ‖ L ∞ ( H 1 ) , ‖ ϕ − ϕ h ‖ L 2 ( L 2 ) , and ‖ p − p h ‖ L 2 ( L 2 ) are of the order O ( h ) with ( V 1 , h , V 2 , h ) -RT0-P0 for λ = 1 , while Tables 2 and 3 show that the errors are almost the same with Table 1 for the two larger λ , even for λ = 1 × 1 0 8 . On the other hand, Tables 2 and 4 display the subtle difference of ( V 1 , h , V 2 , h ) -RT0-P0 and ( V 2 , h , V 1 , h ) -RT0-P0, but they yield the same rates. In particular, Tables 5 and 6 illustrate that the error ‖ ϕ − ϕ h ‖ L 2 ( L 2 ) are of the order O ( h 2 ) with ( V 1 , h , V 2 , h ) -BDM1-P0 for the two larger λ , except the errors ‖ u − u h ‖ L ∞ ( H 1 ) and ‖ p − p h ‖ L 2 ( L 2 ) of the order O ( h ) . This is a interesting discovery, but rarely mentioned in the previous literature. Therefore, our method has the optimal rates for the aforementioned variables and is free of volumetric locking for nearly incompressible porous skeleton.

Table 1

Convergence rates with ( V 1 , h , V 2 , h ) -RT0-P0 for λ = 1

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.365346 × 1 0 00	—	3.324706 × 1 0 ‒ 1	—	8.725866 × 1 0 ‒ 2	—
1/8	1/20	2.478010 × 1 0 00	0.82	1.657478 × 1 0 ‒ 1	1.00	4.295737 × 1 0 ‒ 2	1.02
1/16	1/40	1.297575 × 1 0 00	0.94	8.281841 × 1 0 ‒ 2	1.00	2.149890 × 1 0 ‒ 2	1.00
1/32	1/80	6.608179 × 1 0 ‒ 1	0.97	4.140140 × 1 0 ‒ 2	1.00	1.075622 × 1 0 ‒ 2	1.00
1/64	1/160	3.329846 × 1 0 ‒ 1	0.99	2.069981 × 1 0 ‒ 2	1.00	5.379108 × 1 0 ‒ 3	1.00

Table 2

Convergence rates with ( V 1 , h , V 2 , h ) -RT0-P0 for λ = 1 × 1 0 4

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.598895 × 1 0 00	—	3.320902 × 1 0 ‒ 1	—	8.739410 × 1 0 ‒ 2	—
1/8	1/20	2.543699 × 1 0 00	0.85	1.655511 × 1 0 ‒ 1	1.00	4.296939 × 1 0 ‒ 2	1.02
1/16	1/40	1.324786 × 1 0 00	0.94	8.278514 × 1 0 ‒ 2	1.00	2.150007 × 1 0 ‒ 2	1.00
1/32	1/80	6.738277 × 1 0 ‒ 1	0.98	4.139474 × 1 0 ‒ 2	1.00	1.075619 × 1 0 ‒ 2	1.00
1/64	1/160	3.394360 × 1 0 ‒ 1	0.99	2.069765 × 1 0 ‒ 2	1.00	5.379006 × 1 0 ‒ 3	1.00

Table 3

Convergence rates with ( V 1 , h , V 2 , h ) -RT0-P0 for λ = 1 × 1 0 8

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.599073 × 1 0 00	—	3.320918 × 1 0 ‒ 1	—	8.739421 × 1 0 ‒ 2	—
1/8	1/20	2.543739 × 1 0 00	0.85	1.655512 × 1 0 ‒ 1	1.00	4.296940 × 1 0 ‒ 2	1.02
1/16	1/40	1.324801 × 1 0 00	0.94	8.278516 × 1 0 ‒ 2	1.00	2.150007 × 1 0 ‒ 2	1.00
1/32	1/80	6.738341 × 1 0 ‒ 1	0.98	4.139475 × 1 0 ‒ 2	1.00	1.075619 × 1 0 ‒ 2	1.00
1/64	1/160	3.394391 × 1 0 ‒ 1	0.99	2.069765 × 1 0 ‒ 2	1.00	5.379006 × 1 0 ‒ 3	1.00

Table 4

Convergence rates with ( V 2 , h , V 1 , h ) -RT0-P0 for λ = 1 × 1 0 4

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.594397 × 1 0 00	—	3.320933 × 1 0 ‒ 1	—	8.739432 × 1 0 ‒ 2	—
1/8	1/20	2.540725 × 1 0 00	0.85	1.655514 × 1 0 ‒ 1	1.00	4.296941 × 1 0 ‒ 2	1.02
1/16	1/40	1.323119 × 1 0 00	0.94	8.278518 × 1 0 ‒ 2	1.00	2.150007 × 1 0 ‒ 2	1.00
1/32	1/80	6.729600 × 1 0 ‒ 1	0.98	4.139475 × 1 0 ‒ 2	1.00	1.075619 × 1 0 ‒ 2	1.00
1/64	1/160	3.389963 × 1 0 ‒ 1	0.99	2.069765 × 1 0 ‒ 2	1.00	5.379006 × 1 0 ‒ 3	1.00

Table 5

Convergence rates with ( V 1 , h , V 2 , h ) -BDM1-P0 for λ = 1 × 1 0 4

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.598895 × 1 0 00	—	1.289337 × 1 0 ‒ 1	—	8.968894 × 1 0 ‒ 2	—
1/8	1/20	2.543699 × 1 0 00	0.85	3.162615 × 1 0 ‒ 2	2.03	4.329664 × 1 0 ‒ 2	1.05
1/16	1/40	1.324786 × 1 0 00	0.94	7.946882 × 1 0 ‒ 3	1.99	2.154265 × 1 0 ‒ 2	1.01
1/32	1/80	6.738277 × 1 0 ‒ 1	0.98	1.991767 × 1 0 ‒ 3	2.00	1.076157 × 1 0 ‒ 2	1.00
1/64	1/160	3.394360 × 1 0 ‒ 1	0.99	4.983715 × 1 0 ‒ 4	2.00	5.379680 × 1 0 ‒ 3	1.00

Table 6

Convergence rates with ( V 1 , h , V 2 , h ) -BDM1-P0 for λ = 1 × 1 0 8

h	Δ t	‖ u ‒ u h ‖ L ∞ ( H 1 )	Rate	‖ ϕ ‒ ϕ h ‖ L 2 ( L 2 )	Rate	‖ p ‒ p h ‖ L 2 ( L 2 )	Rate
1/4	1/10	4.599073 × 1 0 00	—	1.289388 × 1 0 ‒ 1	—	8.968914 × 1 0 ‒ 2	—
1/8	1/20	2.543739 × 1 0 00	0.85	3.162772 × 1 0 ‒ 2	2.03	4.329668 × 1 0 ‒ 2	1.05
1/16	1/40	1.324801 × 1 0 00	0.95	7.947305 × 1 0 ‒ 3	1.99	2.154265 × 1 0 ‒ 2	1.01
1/32	1/80	6.738341 × 1 0 ‒ 1	0.98	1.991875 × 1 0 ‒ 3	2.00	1.076157 × 1 0 ‒ 2	1.00
1/64	1/160	3.394391 × 1 0 ‒ 1	0.99	4.983989 × 1 0 ‒ 4	2.00	5.379681 × 1 0 ‒ 3	1.00

6.2 Overcoming the nonphysical pressure oscillations

To show the performance of our method to overcome NPOs, we test the benchmark cantilever bracket problem [15]. For the porous fluid, the no-flow flux boundary conditions are applied along all sides of the unit area [ 0 , 1 ] × [ 0 , 1 ] ; for the porous skeleton, the deformation is fixed along the left edge, a downward traction force is imposed along the top edge, and the right and bottom edges are traction free. We assume that the initial conditions u 0 = 0 and p 0 = 0 . Also, we choose K = 1 × 1 0 − 7 , 1 × 1 0 − 9 , 1 × 1 0 − 12 , and parameters: c 0 = 0 , α = 0.93 , E = 1 0 5 , ν = 0.4 .

In Figures 2, 3, 4, we plot the numerical pressure with ( V 1 , h , V 2 , h ) -RT0-P0, ( V 1 , h , V 2 , h ) -BDM1-P0, and ( P 1 ) 2 -RT0-P0 for with respect different K at the first step (i.e., Δ t = 0.001 ). It is observed that ( P 1 ) 2 -RT0-P0 encounters NPOs that almost concentrate in the left half of the domain, while there is no spurious oscillations for our method whether K = 1 × 1 0 − 9 or K = 1 × 1 0 − 12 . Therefore, such method is free of NPOs.

$Figure 2 The numerical pressure with ( V 1 , h , V 2 , h ) \left({V}_{1,h},{V}_{2,h}) -RT0-P0 for K = 1 × 1 0 ‒ 7 {\bf{K}}=1\times 1{0}^{‒7} .$

Figure 2

The numerical pressure with ( V 1 , h , V 2 , h ) -RT0-P0 for K = 1 × 1 0 ‒ 7 .

$Figure 3 The numerical pressure with ( V 1 , h , V 2 , h ) \left({V}_{1,h},{V}_{2,h}) -RT0-P0 (1eft) and ( P 1 ) 2 {\left({P}_{1})}^{2} -RT0-P0 (right) for K = 1 × 1 0 ‒ 9 {\bf{K}}=1\times 1{0}^{‒9} .$

Figure 3

The numerical pressure with ( V 1 , h , V 2 , h ) -RT0-P0 (1eft) and ( P 1 ) 2 -RT0-P0 (right) for K = 1 × 1 0 ‒ 9 .

$Figure 4 The numerical pressure with ( V 1 , h , V 2 , h ) \left({V}_{1,h},{V}_{2,h}) -BDM1-P0 (1eft) and ( P 1 ) 2 {\left({P}_{1})}^{2} -RT0-P0 (right) for K = 1 × 1 0 ‒ 12 {\bf{K}}=1\times 1{0}^{‒12} .$

Figure 4

The numerical pressure with ( V 1 , h , V 2 , h ) -BDM1-P0 (1eft) and ( P 1 ) 2 -RT0-P0 (right) for K = 1 × 1 0 ‒ 12 .

7 Conclusion

In this article, we have proposed and analyzed two novel nonconforming finite element methods for the three-field Biot’s consolidation model in poroelasticity. These methods are crafted to satisfy the discrete Korn’s inequality inherently, eschewing the need for additional stabilization terms. They are characterized by optimal convergence rates for the approximations of displacement, fluid flux, and pressure. Remarkably, they also effectively mitigate volumetric locking and suppress nonphysical pressure oscillations, thereby enhancing the robustness of numerical solutions in poroelasticity.

Acknowledgements

The authors are very grateful to Editors and anonymous reviewers for their valuable suggestions on an earlier version.

Funding information: The authors extend their appreciation to the supported by the Scientific Research Foundation of Chengdu University of Information Technology KYQN202324; the Key Laboratory of Numerical Simulation of Sichuan Provincial Universities KLNS-2023SZFZ002; the Sichuan National Applied Mathematics Center-Chengdu University of Information Technology 2024ZX001 the Scientific Research Fund of Hunan Provincial Education Department 22A0483.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved it submission.
Conflict of interest: The authors declare no potential conflict of interest.
Ethical approval: This article does not contain any studies with human participants or animals performed by any of the authors.
Data availability statement: Data on the results of the study may be obtained from the corresponding author upon reasonable request.

References

[1] K. Terzaghi, Theoretical Soil Mechanics, Chapman, 1959, DOI: https://doi.org/10.1002/9780470172766. 10.1002/9780470172766Suche in Google Scholar

[2] M. A. Biot, General theory of three dimensional consolidation, J. Appl. Phys. 12 (1941), no. 2, 155–164, DOI: https://doi.org/10.1063/1.1712886. 10.1063/1.1712886Suche in Google Scholar

[3] A. Anandarajah, Computational Methods in Elasticity and Plasticity, Springer, New York, 2010. 10.1007/978-1-4419-6379-6Suche in Google Scholar

[4] D. Boffi, M. Botti, and D. A. DiPietro, A nonconforming high-order method for the Biot problem on general meshes, SIAM J. Sci. Comput. 38 (2016), no. 3, A1508–A1537, DOI: https://doi.org/10.1137/15M1025505. 10.1137/15M1025505Suche in Google Scholar

[5] Y. Chen, Y. Luo, and M. Feng, Analysis of a discontinuous Galerkin method for the Biot’s consolidation problem, Appl. Math. Comput. 219 (2013), no. 17, 9043–9056, DOI: https://doi.org/10.1016/j.amc.2013.03.104. 10.1016/j.amc.2013.03.104Suche in Google Scholar

[6] M. A. Murad, V. Thomée, and A. F. Loula, Asymptotic behavior of semidiscrete finite-element approximations of Biot’s consolidation problem, SIAM J. Numer. Anal. 33 (1996), no. 3, 1065–1083, DOI: https://doi.org/10.1137/0733052. 10.1137/0733052Suche in Google Scholar

[7] L. Berger, R. Bordas, D. Kay, and S. Tavener, Stabilized lowest-order finite element approximation for linear three-field poroelasticity, SIAM J. Sci. Comput. 37 (2015), no.5, A2222–A2245, DOI: https://doi.org/10.1137/15M1009822. 10.1137/15M1009822Suche in Google Scholar

[8] J. Kraus, P. L. Lederer, M. Lymbery, and J. Schöberl, Uniformly well-posed hybridized discontinuous Galerkin/hybrid mixed discretizations for Biot’s consolidation model, Comput. Methods Appl. Mech. Engrg. 384 (2021), no. 23, 113991, DOI: https://doi.org/10.1016/j.cma.2021.113991. 10.1016/j.cma.2021.113991Suche in Google Scholar

[9] S. Kumar, R. Oyarzúa, R. Ruiz-Baier, and R. Sandilya, Conservative discontinuous finite volume and mixed schemes for a new four-field formulation in poroelasticity, ESAIM Math. Model. Numer. Anal. 54 (2020), no. 1, 273–299, DOI: https://doi.org/10.1051/m2an/2019063. 10.1051/m2an/2019063Suche in Google Scholar

[10] X. Feng, Z. Ge, and Y. Li, Analysis of a multiphysics finite element method for a poroelasticity model, IMA J. Numer. Anal. 38 (2018), no. 1, 330–359, DOI: https://doi.org/10.1093/imanum/drx003. 10.1093/imanum/drx003Suche in Google Scholar

[11] Q. Hong, J. Kraus, M. Lymbery, and F. Philo, Conservative discretizations and parameter-robust preconditioners for Biot and multiple-network flux-based poroelasticity models, Numer. Linear Algebra Appl. 26 (2019), no. 4, e2242, DOI: https://doi.org/10.1002/nla.2242. 10.1002/nla.2242Suche in Google Scholar

[12] S. Y. Yi, A coupling of nonconforming and mixed finite element methods for Biot’s consolidation model, Numer. Methods Partial Differential Equations 29 (2013), no. 5, 1749–1777, DOI: https://doi.org/10.1002/num.21775. 10.1002/num.21775Suche in Google Scholar

[13] J. Guo, and M. Feng, A robust and mass conservative virtual element method for linear three-field poroelasticity, J. Sci. Comput. 92 (2022), no. 95, 1–32, DOI: https://doi.org/10.1007/s10915-022-01960-2. 10.1007/s10915-022-01960-2Suche in Google Scholar

[14] J. J. Lee, K. A. Mardal, and R. Winther, Parameter-robust discretization and preconditioning of Biot’s consolidation model, SIAM J. Sci. Comput. 39 (2017), no. 1, A1–A24, DOI: https://doi.org/10.1137/15M1029473. 10.1137/15M1029473Suche in Google Scholar

[15] S. Y. Yi, A study of two modes of locking in poroelasticity, SIAM J. Numer. Anal. 55 (2017), no. 4, 1915–1936, DOI: https://doi.org/10.1137/16M1056109. 10.1137/16M1056109Suche in Google Scholar

[16] C. Rodrigo, F. J. Gaspar, X. Hu, and L. T. Zikatanov, Stability and monotonicity for some discretizations of the Biot’s consolidation model, Comput. Methods Appl. Mech. Engrg. 298 (2016), 183–204, DOI: https://doi.org/10.1016/j.cma.2015.09.019. 10.1016/j.cma.2015.09.019Suche in Google Scholar

[17] M. A. Murad and F. D. Loula, On stability and convergence of finite element approximations of Biot’s consolidation problem, Internat. J. Numer. Methods Engrg. 37 (1994), no. 4, 645–667, DOI: https://doi.org/10.1002/nme.1620370407. 10.1002/nme.1620370407Suche in Google Scholar

[18] P. J. Phillips and M. F. Wheeler, Overcoming the problem of locking in linear elasticity and poroelasticity: an heuristic approach, Comput. Geosci. 13 (2009), 5–12, DOI: https://doi.org/10.1007/s10596-008-9114-x. 10.1007/s10596-008-9114-xSuche in Google Scholar

[19] P. Hansbo and M. G. Larson, Discontinuous Galerkin methods for incompressible and nearly incompressible elasticity by Nitscheas method, Comput. Methods Appl. Mech. Engrg. 191 (2002), no. 17, 1895–1908, DOI: https://doi.org/10.1016/S0045-7825(01)00358-9. 10.1016/S0045-7825(01)00358-9Suche in Google Scholar

[20] J. J. Lee, Robust error analysis of coupled mixed methods for Biot’s consolidation model, J. Sci. Comput. 69 (2016), no. 2, 610–632, DOI: https://doi.org/10.1007/s10915-016-0210-0. 10.1007/s10915-016-0210-0Suche in Google Scholar

[21] Q. Hong and J. Kraus, Parameter-robust stability of classical three-field formulation of Biot’s consolidation model, Electron. T. Numer. Ana. 48 (2018), 202–226, DOI: https://doi.org/10.1553/etna_vol48s202.10.1553/etna_vol48s202Suche in Google Scholar

[22] G. Aguilar, F. Gaspar, F. Lisbona, and C. Rodrigo, Numerical stabilization of Biot’s consolidation model by a perturbation on the flow equation, Internat. J. Numer. Methods Engrg. 75 (2008), no. 11, 1282–1300, DOI: https://doi.org/10.1002/nme.2295. 10.1002/nme.2295Suche in Google Scholar

[23] G. Fu, A high-order HDG method for the Biot’s consolidation model, Comput. Math. Appl. 77 (2019), no. 1, 237–252, DOI: https://doi.org/10.1016/j.camwa.2018.09.029. 10.1016/j.camwa.2018.09.029Suche in Google Scholar

[24] X. Hu, C. Rodrigo, F. J. Gaspar, and L. T. Zikatanov, A nonconforming finite element method for the Biot’s consolidation model in poroelasticity, J. Comput. Appl. Math. 310 (2017), 143–154, DOI: https://doi.org/10.1016/j.cam.2016.06.003. 10.1016/j.cam.2016.06.003Suche in Google Scholar

[25] S. Y. Yi, Convergence analysis of a new mixed finite element method for Biot’s consolidation model, Numer. Methods Partial Differential Equations 30 (2014), no. 4, 1189–1210, DOI: https://doi.org/10.1002/num.21865. 10.1002/num.21865Suche in Google Scholar

[26] T. Bærland, J. J. Lee, K. A. Mardal, and R. Winther, Weakly imposed symmetry and robust preconditioners for Biot’s consolidation model, Comput. Methods Appl. Math. 17 (2017), no. 3, 377–396, DOI: https://doi.org/10.1515/cmam-2017-0016. 10.1515/cmam-2017-0016Suche in Google Scholar

[27] A. Khan, and C. E. Powell, Parameter-robust stochastic Galerkin mixed approximation for linear poroelasticity with uncertain inputs, SIAM J. Sci. Comput. 43 (2021), no. 4, B855–B883, DOI: https://doi.org/10.1137/20M1324296. 10.1137/20M1324296Suche in Google Scholar

[28] M. Botti, D. A. Di Pietro, O. L. Maítre, and P. Sochala, Numerical approximation of poroelasticity with random coefficients using polynomial chaos and hybrid high-order methods, Comput. Methods Appl. Mech. Engrg. 25 (2020), no. 361, 112736, DOI: https://doi.org/10.1016/j.cma.2019.112736. 10.1016/j.cma.2019.112736Suche in Google Scholar

[29] S. C. Brenner and L. Sung, Linear finite element methods for planar linear elasticity, Math. Comp. 59 (1992), no. 200, 321–338, DOI: https://doi.org/10.1090/S0025-5718-1992-1140646-2. 10.1090/S0025-5718-1992-1140646-2Suche in Google Scholar

[30] P. J. Phillips and M. F. Wheeler, A coupling of mixed and discontinuous Galerkin finite element methods for poroelasticity, Comput. Geosci. 12 (2008), no. 4, 417–435, DOI: https://doi.org/10.1007/s10596-008-9082-1. 10.1007/s10596-008-9082-1Suche in Google Scholar

[31] S. C. Brenner and L. R. Scott, The Mathematical Theory of Finite Element Methods, Springer, 2008, DOI: https://doi.org/10.1007/978-0-387-75934-0. 10.1007/978-0-387-75934-0Suche in Google Scholar

[32] R. E. Showalter, Diffusion in poro-elastic media, J. Math. Anal. Appl. 251 (2000), no. 1, 310–340, DOI: https://doi.org/10.1006/jmaa.2000.7048. 10.1006/jmaa.2000.7048Suche in Google Scholar

[33] R. Kouhia and R. Stenberg, A linear nonconforming finite element method for nearly incompressible elasticity and Stokes flow, Comput. Methods Appl. Mech. Engrg. 124 (1995), no. 3, 195–212, DOI: https://doi.org/10.1016/0045-7825(95)00829-P. 10.1016/0045-7825(95)00829-PSuche in Google Scholar

[34] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, 1991, DOI: https://doi.org/10.1002/nme.1620191115. 10.1007/978-1-4612-3172-1Suche in Google Scholar

[35] M. Zhang and S. Zhang, A 3D conforming-nonconforming mixed finite element for solving symmetric stress Stokes equations, Int. J. Numer. Anal. Model. 14 (2017), no. 4, 730–743. Suche in Google Scholar

[36] B. Rivière, Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations, Society for Industrial and Applied Mathematics, 2008, DOI: https://doi.org/10.1137/1.9780898717440. 10.1137/1.9780898717440Suche in Google Scholar

[37] P. Hansbo and M. G. Larson, Discontinuous Galerkin and the Crouzeix-Raviart element: application to elasticity, ESAIM Math. Model. Numer. Anal. 37 (2003), no. 1, 63–72, DOI: https://doi.org/10.1051/m2an:2003020. 10.1051/m2an:2003020Suche in Google Scholar

[38] S. C. Brenner, Forty years of the Crouzeix-Raviart element, Numer. Methods Partial Differential Equations 31 (2015), no. 2, 367–396, DOI: https://doi.org/10.1002/num.21892. 10.1002/num.21892Suche in Google Scholar

Received: 2023-12-31

Revised: 2024-06-18

Accepted: 2024-09-22

Published Online: 2024-12-11

This work is licensed under the Creative Commons Attribution 4.0 International License.

Artikel in diesem Heft

https://doi.org/10.1515/dema-2024-0073

Schlagwörter für diesen Artikel

Biot’s consolidation model; nonconforming finite element method; poroelasticity locking; locking-free

Creative Commons

BY 4.0