Abstract

In this paper, we introduce and analyze H(div)-conforming finite element methods for a nonlinear model in poroelasticity. More precisely, the flow variables are discretized by H(div)-conforming mixed finite elements, while the elastic displacement is approximated by the H(div)-conforming finite element with the interior penalty discontinuous Galerkin formulation. Optimal a priori error estimates are derived for both semidiscrete and fully discrete schemes.

1. Introduction

The Biot poroelasticity model [1] describes the phenomena of coupled mechanics and flow in porous media, and it is of increasing importance today in a divergence of science and engineering applications. Applications of poroelasticity have been made to areas that include carbon sequestration in environment engineering, predictive ability in earthquake engineering, surface subsidence in field phenomena, logging technologies in reservoir engineering, and pathological condition in biomechanics.

In the context of linear poroelasticity systems, a large number of numerical methods have been developed in recent years. Among them, finite element methods (FEMs) are commonly used approaches for such coupled system. In [2–4], Taylor–Hood elements were applied for the approximation of the displacement field and fluid pressure. Recently, some techniques based on mixed formulations have been developed. A method coupling standard conforming FEM for the displacement with mixed FEMs for flow variables has been provided in [5, 6]. However, it is well known that a conforming FEM approximation of displacement may give rise to locking or nonphysical oscillations [7]. Some remedies of locking include nonconforming FEMs [8–10], mixed FEMs [11, 12], discontinuous Galerkin (DG) methods [13–15], hybrid high-order methods [16], and weak Galerkin methods [17–20]. We also refer the interested reader to [21–24] for further study of locking-free FEMs on linear poroelasticity. Moreover, some preconditioners that are robust with respect to the physical and discretization parameters have been developed in this direction (refer to [25–28] and the references therein).

Compared with linear model in poroelasticity, the numerical methods on nonlinear problem are still very rare. The goal of this paper is to develop H(div)-conforming finite element methods for a nonlinear poroelasticity model. The nonlinearity arises from the dependence of the hydraulic conductivity on the dilation , that is, . Such model has been studied in [29], where the existence and uniqueness of the weak solution were proved. Moreover, the authors in [29] introduce a FEM approximation and give the optimal order error estimate. We also find that similar problem has been treated in [30], where a mixed method combining conforming and mixed FEMs was proposed and analyzed. Some other works of numerical methods for nonlinear poroelasticity model can be found in [31, 32]. In the current work, we will propose and analyze H(div)-conforming methods to solve the nonlinear problem studied in [29]. More precisely, the flow variables are discretized by H(div)-conforming mixed FEMs, while the displacement is approximated by the H(div)-conforming FEM with interior penalty DG method. Combining H(div)-conforming finite elements with DG methods was initially proposed in [33, 34] (see also [35–37]); they mainly intended to solve Stokes and Navier–Stokes equations in fluid mechanics. Later, this method was extended to more complex Darcy–Stokes interface problems [38–40], Brinkman problems [41], a magnetic induction model [42], and linear Biot model in poroelasticity [43]. We also refer the interested reader to [44] on robust numerical simulations of H(div)-conforming FEMs for Stokes, coupled Darcy–Stokes, and Brinkman problems. Our work is to extend H(div)-conforming FEMs of linear Biot model to the nonlinear case. With the help of the framework presented in [30], we will give a detailed a priori error analysis of the proposed numerical scheme. It is worth mentioning that we shall address some difficulties arising from the inherent nonlinearity (see the terms (43), (44), (68), and (69) below). The main advantage of our method is that we present a unified treatment for both flow variables and the displacement. In addition, we relax conformity by using H(div) space to approximate the displacement, and thus the normal components are continuous and the tangential components are treated by interior penalty DG approach; this implies that our method is locally conservative.

The rest of this paper is organized as follows. In Section 2, we describe the nonlinear poroelasticity model and present its mixed weak formulation. We propose a spatial semidiscrete scheme which is based on H(div)-conforming finite elements in Section 3, where a priori error estimate is also provided. In Section 4, a fully discrete scheme with the backward Euler time stepping is proposed and analyzed. Some conclusions are made in Section 5.

2. Preliminaries and Model Problem

2.1. The Model Problem

Let be a bounded convex polygonal domain, with Lipschitz boundary . We consider the following nonlinear poroelastic model in over a time interval (see [29]):wherein which is the displacement of the solid phase, is the fluid pressure, is the total stress tensor, is the storage coefficient, is the Biot–Willis constant, and and are Lamé constants; is the hydraulic conductivity that satisfies assumptions (6) and (7) below.

Denote by (resp. ) the Dirichlet (resp. traction) boundary for the elastic variables. Similarly, denote by (resp. ) the pressure Dirichlet (resp. fluid normal flux) boundary. We assume that with , and . The boundary conditions and initial conditions for the above system read aswhere denotes the outward unit normal vector.

Setting , we can reformulate equation (1) by

2.2. Weak Formulation

We start by introducing some notation. Let be the standard Sobolev space, equipped with norm and seminorm . If , is understood as . Denote and . When , we will omit the index . For the space , its norm is still denoted by . A subspace of with vanishing trace on is given by . Additionally, we define with its graph norm . Two subspaces are then introduced by and . For convenience, we set , and .

Integrating by parts, we deduce that the mixed weak formulation for (4) reads as follows: find for every such thatfor any , where . Additionally, we assume that there exist positive and such thatand is Lipschitz continuous:

To deal with functions of time and space, we introduce the standard Bochner space , which consists of all functions with normfor . When , the norm is defined as

3. Error Estimates for the Semidiscrete Scheme

In this section, we aim to present the semidiscrete numerical method which focuses on discretizing the spatial variables. Let be a shape-regular decomposition of into triangles . denotes the diameter of and . Denote by the set of interior edges of elements in , the set of boundary edges on , and the set of boundary edges on . Therefore, the set of all edges . The length of the edge is denoted by . For every edge , we fix a unit normal such that for edges on the boundary , is the outward unit normal.

For shared by elements and , let be a scalar or vector piecewise smooth function and set , and we define the average and jump by

On a boundary edge, we set

Consider the following finite element spaces:where is the H(div)-conforming space introduced by Brezzi, Douglas, and Marini and denotes the space of polynomials of degree less than or equal to on .

Let be the BDM interpolation and be the projection. It is well known that they satisfy the following properties [45, 46]:

Here and throughout the paper, we utilize to represent a positive generic constant that is independent of and but may take different values at different occurrences.

3.1. Semidiscrete H(div)-Conforming Methods with DG Formulation

The expression of below is similar to the statement in Section 3.1 in the work [43]. For completeness, we derive its scheme as follows.

Multiplying the third equation in (4) by any on each element , integrating by parts, and then summing over all elements in , we arrive atwhere is outward normal of .

We first have

Here in the third line, we have used in , and in the last line, we have used the fact that on and the regularity of the exact solutions.

In addition, it is easy to check that

Substituting (21) and (22) into (20) yields

Let be the unit tangential direction to the edge so that and form a right-hand side coordinate system; it follows that for any , we have

Therefore,

Combining (18) and (19) and the regularity of the exact solutions implies that

Here in the third line, we have used the fact that on . Therefore, (23) is reduced to

Adding some stabilized terms as in [47], we propose the following DG method for the third equation in (4):with

Consequently, after integrating by parts, we find that the exact solutions of (4) satisfyfor any .

Thus, the corresponding semidiscrete H(div)-conforming FEMs for (4) can be designed as follows: find such thatfor any , with the initial conditions given by

3.2. Error Estimates

We begin by defining the mesh-dependent norm on :

Then, from [43], we can obtain the following result which shows the continuity and coercivity of the bilinear .

Lemma 1. There exists a constant such that

Moreover, for sufficiently large penalty parameter , it holds that

We are now in a position to state the following error estimate, which is the main result of this section.

Theorem 1. Let and be the solutions of (4) and (31), respectively. Then, the following holds:where and depends on the regularity of , and .

Proof. Subtracting (30) from (31) yieldsfor any .
We then split the error into with and . Similarly, with and ; with and . Since the estimates for , , and can be derived by the interpolation error bounds in (16), (18), and (19), it remains to estimate , , and . To this end, using (15) and (17), we can rewrite (37) asfor any .
Setting , , and in the above equations, we infer thatSumming equation (39), integrating in time from to , and noting that , , andwe obtainwithFor the first term , using (6) and Cauchy–Schwarz and Young’s inequalities, we find thatwhere is an arbitrarily small number.
For the second term , with the aid of (6) and (7) and Cauchy–Schwarz and Young’s inequalities, we deduce thatwhere is an arbitrarily small number and .
To estimate the third term , integrating by parts and noting that , we first obtainwhich combined with (34) and Young’s inequality implies thatwhere is an arbitrarily small number.
Combining the bounds in (41)–(46) and using (6) and (35), we haveNow, we choose appropriate such that . Consequently, the above inequality still holds if we replace its left-hand side term by . Then, dividing both sides of the above inequality by and using Gronwall’s inequality, we haveNoting that the above estimate holds for all and using appropriate approximation properties stated in (16), (18), and (19), we obtainThis can be stated by the following equivalent formulation:This together with the standard interpolation error estimates for , , and and the triangle inequality yields the desired estimate (36).