Convergence Analysis of H(div)-Conforming Finite Element Methods for a Nonlinear Poroelasticity Problem

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


Introduction
e 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][3][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][9][10], mixed FEMs [11,12], discontinuous Galerkin (DG) methods [13][14][15], hybrid high-order methods [16], and weak Galerkin methods [17][18][19][20]. We also refer the interested reader to [21][22][23][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][26][27][28] and the references therein). Compared with linear model in poroelasticity, the numerical methods on nonlinear problem are still very rare. e goal of this paper is to develop H(div)-conforming finite element methods for a nonlinear poroelasticity model. e nonlinearity arises from the dependence of the hydraulic conductivity κ on the dilation ∇ · u, that is, κ � κ(∇ · u). 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][36][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][39][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). e main advantage of our method is that we present a unified treatment for both flow variables and the displacement. In addition, we relax H 1 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. e 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.

e Model Problem.
Let Ω ⊂ R 2 be a bounded convex polygonal domain, with Lipschitz boundary zΩ. We consider the following nonlinear poroelastic model in Ω over a time interval (0, T] (see [29]): in which u(x, t) is the displacement of the solid phase, p(x, t) is the fluid pressure, σ(x, t) is the total stress tensor, c 0 > 0 is the storage coefficient, 0 < α ≤ 1 is the Biot-Willis constant, and λ and μ are Lamé constants; κ(·) is the hydraulic conductivity that satisfies assumptions (6) and (7) below. Denote by Γ d (resp. Γ t ) the Dirichlet (resp. traction) boundary for the elastic variables. Similarly, denote by Γ p (resp. Γ f ) the pressure Dirichlet (resp. fluid normal flux) boundary. We assume that zΩ � Γ d ∪Γ t with |Γ d | > 0, and zΩ � Γ p ∪Γ f . e boundary conditions and initial conditions for the above system read as where n denotes the outward unit normal vector.

Weak Formulation.
We start by introducing some notation. Let H s (D) (s ≥ 0) be the standard Sobolev space, equipped with norm ‖ · ‖ s,D and seminorm | · | s,D . If s � 0, H 0 (D) is understood as L 2 (D). Denote (w, v) D � D wv dx and 〈w, v〉 zD � zD wvd F . When D � Ω, we will omit the index Ω. For the space H s (D) � (H s (D)) 2 , its norm is still denoted by ‖ · ‖ s,D . A subspace of H 1 (Ω) with vanishing trace on Γ d is given by Integrating by parts, we deduce that the mixed weak formulation for (4) reads as follows: find (p, q, u ∈ P × Q × V) for every t ∈ (0, T] such that . Additionally, we assume that there exist positive k a and k b such that and κ − 1 is Lipschitz continuous: 2 Discrete Dynamics in Nature and Society To deal with functions of time and space, we introduce the standard Bochner space L p (0, T; H s (Ω)), which consists of all functions u:

Error Estimates for the Semidiscrete Scheme
In this section, we aim to present the semidiscrete numerical method which focuses on discretizing the spatial variables.
the set of boundary edges on Γ d , and E t h the set of boundary edges on Γ t . erefore, the set of all edges For every edge e ∈ E h , we fix a unit normal n e such that for edges on the boundary zΩ, n e is the outward unit normal.
For e ∈ E I h shared by elements K + and K − , let v be a scalar or vector piecewise smooth function and set v ± � v| e∩zK ± , and we define the average and jump by On a boundary edge, we set Consider the following finite element spaces: where BDM k (k ≥ 1) is the H(div)-conforming space introduced by Brezzi, Douglas, and Marini and P k (K) denotes the space of polynomials of degree less than or equal to k on K.
be the BDM interpolation and P h : P ⟶ P h be the L 2 projection. It is well known that they satisfy the following properties [45,46]: Here and throughout the paper, we utilize C to represent a positive generic constant that is independent of h and Δt but may take different values at different occurrences.

Semidiscrete H(div)-Conforming Methods with DG Formulation.
e expression of a h (·, ·) 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 v ∈ V h on each element K, integrating by parts, and then summing over all elements in T h , we arrive at where n T is outward normal of T.
We first have Here in the third line, we have used σn e � 0 in E t h , and in the last line, we have used the fact that [v · n e ] � 0 on e ∈ E I h ∪E d h and the regularity of the exact solutions. In addition, it is easy to check that Substituting (21) and (22) Let τ e be the unit tangential direction to the edge e so that n e and τ e form a right-hand side coordinate system; it follows that for any w, we have w � w · n e n e + w · τ e τ e . (24) erefore, � ϵ(u)n e · n e n e + ϵ(u)n e · τ e τ e · v · n e n e + v · τ e τ e � ϵ(u)n e · n e v · n e + ϵ(u)n e · τ e v · τ e .
Combining (18) and (19) and the regularity of the exact solutions implies that Here in the third line, we have used the fact that [v · n e ] � 0 on e ∈ E I h ∪E d h . erefore, (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) satisfy us, the corresponding semidiscrete H(div)-conforming FEMs for (4) can be designed as follows: find for any (w, z, v) ∈ P h × Q h × V h , with the initial conditions given by

Error Estimates.
We begin by defining the mesh-de- en, from [43], we can obtain the following result which shows the continuity and coercivity of the bilinear a h (·, ·).

Lemma 1.
ere exists a constant C a > 0 such that Moreover, for sufficiently large penalty parameter β, it holds that 4 Discrete Dynamics in Nature and Society We are now in a position to state the following error estimate, which is the main result of this section. (4) and (31), respectively. en, the following holds: where ‖x‖ L ∞ (0,T;E h ) � sup 0≤s≤T ‖x(s)‖ h and C depends on the regularity of u, u t , p, p t , and q.
Proof. Subtracting (30) from (31) yields We then split the error Since the estimates for ξ p , ξ q , and ξ u can be derived by the interpolation error bounds in (16), (18), and (19), it remains to estimate θ p , θ q , and θ u . To this end, using (15) and (17), we can rewrite (37) as Setting w � θ p , z � θ q , and v � z t θ u in the above equations, we infer that Summing equation (39), integrating in time from 0 to t ≤ T, and noting that θ p (0) � 0, θ u (0) � 0, and we obtain with For the first term A 1 , using (6) and Cauchy-Schwarz and Young's inequalities, we find that where ε 1 is an arbitrarily small number.

Discrete Dynamics in Nature and Society
For the second term A 2 , with the aid of (6) and (7) and Cauchy-Schwarz and Young's inequalities, we deduce that where ϵ 2 is an arbitrarily small number and C q � max 0≤s≤t ‖Π h q(s)‖ ∞ . To estimate the third term A 3 , integrating by parts and noting that θ u (0) � 0, we first obtain which combined with (34) and Young's inequality implies that where ε 3 is an arbitrarily small number.
Combining the bounds in (41)-(46) and using (6) and (35), we have Now, we choose appropriate ε i (i � 1, 2, 3) such that Consequently, the above inequality still holds if we replace its lefthand side term by C min (‖θ u (t)‖ 2 h + ‖θ p (t)‖ 2 0 + t 0 ‖θ q (s)‖ 2 0 ds). 6 Discrete Dynamics in Nature and Society en, dividing both sides of the above inequality by C min and using Gronwall's inequality, we have Noting that the above estimate holds for all 0 ≤ t ≤ T and using appropriate approximation properties stated in (16), (18), and (19), we obtain is can be stated by the following equivalent formulation: is together with the standard interpolation error estimates for ξ p , ξ q , and ξ u and the triangle inequality yields the desired estimate (36). e fully discrete mixed element method with back Euler time stepping reads as follows: at each time t � t n , find (p n h , q n h , u n h ) ∈ P h × Q h × V h such that c 0 z t p n h , w + α ∇ · z t u n h , w + ∇ · q n h , w � g n , w ,

Error
Estimates. e standard Taylor expansion leads to where z t f n � z t f(t n ). Now, we prove the following error estimate which is the main result of this section. (4) and (51), respectively. en, the following finite element error estimate holds: where C depends on the regularity of u, u t , u tt , p, p t , p tt , and q.
Proof. Since (30) holds for the exact solution u at any time t � t n , this together with Taylor expansion (53) yields c 0 z t p n , w + α ∇ · z t u n , w + ∇ · q n , w � g n , w + c 0 Δt Discrete Dynamics in Nature and Society 7 Subtracting (51) from (55), we obtain then we further obtain Similarly, the second term B 2 can be estimated by where ε 4 is an arbitrarily small number. Similarly, we can bound B 4 by Δt θ n q � � � � � � � � � � where ε 5 is an arbitrarily small number and C m � max 1≤n≤m ‖Π h q n ‖ ∞ .
Discrete Dynamics in Nature and Society