Locking free staggered DG method for the Biot system of poroelasticity on general polygonal meshes
aa r X i v : . [ m a t h . NA ] N ov Locking free staggered DG method for the Biot system ofporoelasticity on general polygonal meshes
Lina Zhao ∗ Eric Chung † Eun-Jae Park ‡ November 11, 2020
Abstract
In this paper we propose and analyze a staggered discontinuous Galerkin method for a five-fieldformulation of the Biot system of poroelasticity on general polygonal meshes. Elasticity is equippedwith stress-displacement-rotation formulation with weak stress symmetry for arbitrary polynomialorders, which extends the piecewise constant approximation developed in (L. Zhao and E.-J. Park,SIAM J. Sci. Comput. 42:A2158–A2181,2020). The proposed method is locking free and canhandle highly distorted grids possibly including hanging nodes, which is desirable for practicalapplications. We prove the convergence estimates for the semi-discrete scheme and fully discretescheme for all the variables in their natural norms. In particular, the stability and convergenceanalysis do not need a uniformly positive storativity coefficient. Moreover, to reduce the size ofthe global system, we propose a five-field formulation based fixed stress splitting scheme, wherethe linear convergence of the scheme is proved. Several numerical experiments are carried out toconfirm the optimal convergence rates and the locking-free property of the proposed method.
Keywords:
Staggered DG, General polygonal mesh, Locking-free, Fixed stress splitting, Weaksymmetry, Biot system, Poroelasticity
The Biot system of poroelasticity [3, 41] models fluid flow within deformable porous media, and hasbeen extensively studied in the literature due to its broad range of applications such as geosciences,carbon sequestration and biomedical applications. The system is a fully coupled system, consisting ofan equilibrium equation for the solid and a mass balance equation for the fluid. A great amount ofeffort has been devoted to developing efficient numerical schemes for the Biot system. The two-fielddisplacement-pressure formulation is studied by numerous numerical methods such as finite differencemethod [20], finite volume method [35], finite element methods [31, 32, 33, 40], coupling of HHOmethod with SWIP method [4], weak Galerkin method [22] and HDG method [19]. The three-fielddisplacement-pressure-Darcy velocity formulation with various couplings of continuous and discontin-uous Galerkin methods, and mixed finite element methods is studied in [37, 38, 45, 44, 36, 29, 47, 27].Very recently, a stabilized hybrid mixed finite element method uniformly stable with respect to thephysical parameters is proposed in [34]. Alternatively, fully-mixed formulations of the Biot system hasalso drawn great attention [24, 46, 26, 1]. Typical difficulty encountered in the design of numericalschemes for the Biot system is to avoid locking and to get estimates independent of the storativitycoefficient.Therefore, the purpose of this paper is to develop and analyze a staggered discontinuous Galerkin(DG) method on general polygonal meshes for the Biot system of poroelasticity based on a five-fieldformulation. As a new generation of numerical schemes, staggered DG method offers many salientfeatures such as the easy treatment of general polygonal elements possibly including hanging nodes,the local and global mass conservation and superconvergence. All these features make it a goodcandidate for practical applications, thus it has been extensively studied for a wide range of partial ∗ Department of Mathematics,The Chinese University of Hong Kong, Hong Kong SAR, China. ([email protected]). † Department of Mathematics, The Chinese University of Hong Kong, Hong Kong SAR, China. ([email protected]). ‡ School of Mathematics and Computing (Computational Science and Engineering), Yonsei University, Seoul 03722, Korea.([email protected]). ocking free SDG for poroelasticity c and are valid for c = 0. Moreover, like the linear elastic-ity [53], the error bounds are robust for nearly incompressible materials. The numerical experimentsincluding cantilever bracket benchmark problem presented verify the optimal convergence rates andlocking free property of the proposed method.The rest of the paper is organized as follows. In the next section, we describe the model problemand the associated weak formulation. Then in section 3, we derive the discrete formulation for theBiot system of poroelasticity and prove the unique solvability. Convergence analysis for semi-discreteand fully discrete methods is presented in section 4. In section 5, a fixed stress splitting scheme isintroduced and analyzed. Several numerical experiments are carried out in section 6 to confirm theproposed theories. Finally, a conclusion is given. In this section we introduce the poroelasticity system and its fully mixed weak formulation, which laysfoundation for the derivation of our staggered DG scheme. Let Ω ⊂ R be a simply connected boundeddomain, occupied by a poroelastic media saturated with fluid. Then the stress-strain constitutiverelationship for the poroelastic body is given by Aσ e = ǫ ( u ) , (2.1)where ǫ ( u ) = ( ∇ u + ∇ u T ) and σ e = 2 µǫ ( u ) + λ ∇ · uI . Here A is a symmetric, bounded and uniformlypositive definite linear operator representing the compliance tensor, σ e is the elastic strain, u is thesolid displacement. ocking free SDG for poroelasticity Aσ = 12 µ ( σ − λ µ + 2 λ tr ( σ ) I ) , where I is the 2 × µ > λ ≥ σ e = 2 µǫ ( u ) + λ ∇ · uI . The poroelastic stress, which includes the effect of the fluid pressure p , is given as σ = σ e − αpI, (2.2)where 0 < α ≤ f representing the body forces and a source term q , the quasi-static Biot systemthat governs the fluid flow within the poroelastic media is given as follows: −∇ · σ = f in Ω × [0 , T ] ,K − z + ∇ p = 0 in Ω × [0 , T ] ,∂∂t ( c p + α ∇ · u ) + ∇ · z = q in Ω × [0 , T ] , (2.3)where z is the Darcy velocity, c ≥ T > K isa symmetric and positive definite tensor representing the permeability of the porous media divided bythe fluid viscosity, i.e., there exist two strictly positive real numbers K and K satisfying for almostevery x ∈ Ω and all ξ ∈ R such that | ξ | = 10 < K ≤ K ( x ) ξ · ξ ≤ K . Further, we introduce the global anisotropy ratio ̺ B = K K . (2.4)For the sake of simplicity, we assume that the system satisfies homogeneous Dirichlet boundary con-ditions for u and p , i.e., u = 0 , p = 0 on ∂ Ω and the initial condition p ( x,
0) = p . Now we are readyto derive the weak formulation, we have from (2.1) and (2.2) ∇ · u = tr ( ǫ ( u )) = tr ( Aσ e ) = trA ( σ + αpI ) , which can be substituted in the third equation of (2.3) to get ∂ t ( c p + αtrA ( σ + αpI )) + ∇ · z = q. In the weakly symmetric stress formulation, we allow for σ to be non-symmetric and introduce theLagrange multiplier γ = rot( u ) /
2, rot( u ) = − ∂u ∂y + ∂u ∂x . The constitutive equation (2.1) can berewritten as A ( σ + αpI ) = ∇ u − γχ, where χ = (cid:18) −
11 0 (cid:19) . By the preceding arguments, we can propose the following mixed weak formulation for the Biot systemof poroelasticity: Find ( σ, u, γ, z, p ) ∈ L (Ω) × × H (Ω) × L (Ω) × L (Ω) × H (Ω) such that( A ( σ + αpI ) , ψ ) − ( ∇ u, ψ ) + ( γ, as ( ψ )) = 0 ∀ ψ ∈ L (Ω) × , (2.5)( σ, ∇ v ) = ( f, v ) ∀ v ∈ H (Ω) , (2.6)( as ( σ ) , η ) = 0 ∀ η ∈ L (Ω) , (2.7)( K − z, ξ ) + ( ∇ p, ξ ) = 0 ∀ ξ ∈ L (Ω) , (2.8) c ( ∂ t p, w ) + α ( ∂ t A ( σ + αpI ) , wI ) − ( z, ∇ w ) = ( q, w ) ∀ w ∈ H (Ω) , (2.9) ocking free SDG for poroelasticity trAψ, w ) = ( Aψ, wI ) and as ( ψ ) = − ψ + ψ for ψ ∈ L (Ω) × .For the initial condition we can set z = − K ∇ p which naturally satisfies (2.8). We can alsodetermine ( σ , u , γ ) by solving the elasticity problem (2.5)-(2.7) with p given as initial data. Werefer to the initial data obtained by this procedure as compatible initial data.Before closing this section, we also introduce some notations that will be used later. Let D ⊂ R . By( · , · ) D , we denote the scalar product in L ( D ) : ( p, q ) D := R D p q dx . When D coincides with Ω, the sub-script Ω will be dropped. We use the same notation for the scalar product in L ( D ) and in L ( D ) × .More precisely, ( ξ, w ) D := P i =1 ( ξ i , w i ) for ξ, w ∈ L ( D ) and ( ψ, ζ ) D := P i =1 P j =1 ( ψ i,j , ζ i,j ) D for ψ, ζ ∈ L ( D ) × . The associated norm is denoted by k · k ,D . Given an integer m ≥ n ≥ W m,n ( D ) and W m,n ( D ) denote the usual Sobolev space provided the norm and semi-norm k v k W m,n ( D ) = { P | ℓ |≤ m k D ℓ v k nL n ( D ) } /n , | v | W m,n ( D ) = { P | ℓ | = m k D ℓ v k nL n ( D ) } /n . If n = 2 we usu-ally write H m ( D ) = W m, ( D ) and H m ( D ) = W m, ( D ), k v k H m ( D ) = k v k W m, ( D ) and | v | H m ( D ) = | v | W m, ( D ) . In addition, for an integer l ≥ H l (0 , T ; H m (Ω)) and L ∞ (0 , T ; H m (Ω)) with the norms k ψ k H l (0 ,T ; H m (Ω)) = l X i =0 (cid:16) Z T k d it ψ ( t ) k H m (Ω) dt (cid:17) / , k ψ k L ∞ (0 ,T ; H m (Ω)) = max ≤ t ≤ T k ψ ( t ) k H m (Ω) . In the sequel, we use C to denote a positive constant independent of the meshsize and of the anisotropyratio ̺ B in (2.4), which may have different values at different occurrences. In addition, for any ψ ∈ L (Ω) × , ( Aψ, ψ ) / is a norm equivalent to k ψ k , which will be denoted by k A / ψ k throughoutthe paper. In this section, we present the staggered DG discretization for the Biot system of poroelasticity.Specifically, a fully mixed formulation will be used, and the symmetry of stress is imposed weakly viathe introduction of a lagrange multiplier.To begin with, we introduce the construction of staggered meshes, which lays foundation for theconstruction of finite element spaces suitable for our setting. Let T u be a shape-regular family ofmatching meshes covering Ω exactly, where hanging nodes are allowed. Note that the shape-regularityfor general polygonal meshes can be described as follows: (1) Every element M ∈ T u is star-shapedwith respect to a ball of radius ≥ ρ B h E , where ρ B is a positive constant and h E is the diameter of M ;(2) For every element M ∈ T u and every edge e ⊂ ∂M , it satisfies h e ≥ ρ e h M , where ρ e is a positiveconstant and h e represents the length of edge e . For an element M ∈ T u , we select an interior point ν and connect it to all the vertices of M . The resulting sub-triangulation is denoted as T h . We remarkthat ν is an arbitrary point interior to M and for simplicity we choose it as center point. The unionof all the edges in the primal partition T u is denoted as F u . We use F u to represent the subset of F u , that is the set of edges in F u that does not lie on ∂ Ω. In addition, the edges generated duringthe subdivision process due to the connection of the interior point ν to all the vertices are called dualedges, which is denoted as F p . In addition, we define F = F u ∩ F p and F = F u ∩ F p . For eachtriangle τ ∈ T h , we let h τ be the diameter of τ and h = max { h τ , τ ∈ T h } . Next, we introduce the dualmesh. For each interior edge e ∈ F u , we use D ( e ) to represent the dual mesh, which is the union ofthe two triangles in T h sharing the common edge e . For each edge e ∈ F u \F u , we use D ( e ) to denotethe triangle in T h having the edge e , see Figure 1 for an illustration. For each edge e , we define a unitnormal vector n e as follows: If e is an interior edge, then n e is the unit normal vector of e pointingtowards the outside of Ω. If e is an interior edge, we then fix n e as one of the two possible unit normalvectors on e . When there is no ambiguity, we use n instead of n e to simplify the notation. Let k ≥ τ ∈ T h and e ∈ F , we define P k ( τ ) and P k ( e ) as the spacesof polynomials of degree less than or equal to k on τ and e , respectively.For a scalar, vector, or tensor function v that is double-valued on an interior edge e ∈ F , its jumpand average on e are defined as[ v ] = v | τ − v | τ , { v } = 12 ( v | τ + v | τ ) , ocking free SDG for poroelasticity F u anddashed lines represent edges in F p .where τ and τ are the two triangles belonging to T h which share the common edge e . We set [ v ] = v | τ and { v } = v | τ on boundary edges.Now we are ready to introduce the finite element spaces that will be used for the construction ofour method. First, the locally H (Ω) − conforming SDG space S h is defined by S h := { w ∈ L (Ω) : w | τ ∈ P k ( τ ) , [ w ] | e = 0 ∀ e ∈ F u ; w | ∂ Ω = 0 } , which is equipped by the following norm k w k Z = X τ ∈T h k∇ w k ,τ + X e ∈F p h − e k [ w ] k ,e . For a function v = ( v , v ) ∈ [ S h ] we define k v k h = k v k Z + k v k Z .Next, the locally H (div; Ω)-conforming SDG space Σ h is defined byΣ h := { ξ ∈ L (Ω) : ξ | τ ∈ P k ( τ ) ∀ τ ∈ T h , [ ξ · n ] | e = 0 ∀ e ∈ F p } . Finally, the locally H (Ω) − conforming finite element space is defined as M h := { η ∈ L (Ω) : η | τ ∈ P k ( τ ) ∀ τ ∈ T h , [ η ] | e = 0 ∀ e ∈ F p } . We can derive the staggered DG discretization for the Biot system of poroelasticity by following[51, 53], which reads as follows: Find ( σ h , u h , γ h , z h , p h ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h such that( A ( σ h + αp h I ) , ψ ) − B ∗ h ( u h , ψ ) + ( γ h , as ( ψ )) = 0 ∀ ψ ∈ [Σ h ] , (3.1) B h ( σ h , v ) = ( f, v ) ∀ v ∈ [ S h ] , (3.2)( as ( σ h ) , η ) = 0 ∀ η ∈ M h , (3.3)( K − z h , ξ ) + b ∗ h ( p h , ξ ) = 0 ∀ ξ ∈ Σ h , (3.4) c ( ∂ t p h , w ) + α ( ∂ t A ( σ h + αp h I ) , wI ) − b h ( z h , w ) = ( q, w ) ∀ w ∈ S h , (3.5)where B h ( ψ, v ) = − X e ∈F p ( ψn, [ v ]) e + ( ψ, ∇ v ) ∀ ( ψ, v ) ∈ [Σ h ] × [ S h ] ,B ∗ h ( v, ψ ) = X e ∈F u ( v, [ ψn ]) e − ( v, ∇ · ψ ) ∀ ( ψ, v ) ∈ [Σ h ] × [ S h ] and b h ( ξ, w ) = − X e ∈F p ( ξ · n, [ w ]) e + ( ξ, ∇ w ) ∀ ( ξ, w ) ∈ Σ h × S h ,b ∗ h ( w, ξ ) = X e ∈F u ( w, [ ξ · n ]) e − ( w, ∇ · ξ ) ∀ ( ξ, w ) ∈ Σ h × S h . ocking free SDG for poroelasticity b h ( ξ, w ) = b ∗ h ( w, ξ ) ∀ ( w, ξ ) ∈ S h × Σ h ,B h ( ψ, v ) = B ∗ h ( v, ψ ) ∀ ( ψ, v ) ∈ [Σ h ] × [ S h ] . In addition, the following inf-sup condition holds; cf. [13]. k w k Z ≤ C sup ξ ∈ Σ h b ∗ h ( w, ξ ) k ξ k ∀ w ∈ S h (3.6)and k v k h ≤ C sup ψ ∈ [Σ h ] − B ∗ h ( v, ψ ) k ψ k ∀ v ∈ [ S h ] . (3.7) Remark 3.1. (local conservation on the dual mesh). Taking v = 1 in (3.2) to be identically one overthe dual mesh D ( e ) , ∀ e ∈ F u , we have − ( σ h n D , ∂D ( e ) = ( f, D ( e ) , where n D is the outward unit normal vector of D ( e ). This property is the crux for a posteriori estimatesbased on equilibrated fluxes; cf. [51, 54].Now we prove the inf-sup condition, which is the key to proving the stability and convergenceestimates of the proposed method. The methodology exploited here is different from that of [28]. Indeedthe norm used here bypasses the construction of divergence free functions in proving the convergenceestimates. Lemma 3.1.
The following inf-sup condition holds sup ψ h ∈ [Σ h ] − B ∗ h ( v h , ψ h ) + ( η h , as ( ψ h )) k ψ h k ≥ C ( k v h k h + k η h k ) , ∀ ( v h , η h ) ∈ [ S h ] × M h . (3.8) Proof.
To prove the inf-sup condition (3.8), we apply Fortin interpolation operator (cf. [5]), namely,we need to build an interpolation operator J h : H (Ω) × → [Σ h ] such that it satisfies for any ψ ∈ H (Ω) × − B ∗ h ( v h , ψ − J h ψ ) + ( η h , as ( ψ − J h ψ )) = 0 ∀ ( v h , η h ) ∈ [ S h ] × M h , k J h ψ k ≤ C k ψ k . To this end we will first build ψ ∈ [Σ h ] such that it satisfies B ∗ h ( v h , ψ − ψ ) = 0 ∀ v h ∈ [ S h ] , k ψ k ≤ C k ψ k , (3.9)then we correct ψ by using ψ ∈ [Σ h ] which satisfies B ∗ h ( v h , ψ ) = 0 , ∀ v h ∈ [ S h ] , in addition, thereholds ( as ( ψ − ψ ) , η h ) = ( as ( ψ ) , η h ) ∀ η h ∈ M h , (3.10) k ψ k ≤ C k ψ − ψ k . (3.11)Since B ∗ h ( · , · ) satisfies the inf-sup condition (3.7), which means there exists ψ ∈ [Σ h ] that satisfies B ∗ h ( v h , ψ − ψ ) = 0 ∀ v h ∈ [ S h ] , k ψ k ≤ C k ψ k . Next, we correct ψ in the following way. We define the following spaces for each primal element M ∈ T u : V kh ( M ) = { ς | M ∈ H ( M ) : ς | τ ∈ P k ( τ ) , τ ⊂ M } ,Q kh ( M ) = { q | M ∈ L ( M ) : q | τ ∈ P k ( τ ) , τ ⊂ M, Z M q dx = 0 } . ocking free SDG for poroelasticity M ∈ T u , we consider the stable approximation V k +1 h ( M ) × Q kh ( M ) for theStokes element such that curl V k +1 h ( M ) ⊂ Σ h ( M ) , where Σ h ( M ) denotes Σ h restricted to the primal element M . We then solve for ( ψ ∗ h | M , p ∗ h | M ) ∈ V k +1 h ( M ) × Q kh ( M ) such that(2 µǫ ( ψ ∗ h ) , ǫ ( φ h )) M + ( p ∗ h , ∇ · φ h ) M = ( f, φ h ) M ∀ φ h ∈ V k +1 h ( M ) , ( q h , ∇ · ψ ∗ h ) M = ( as ( ψ − ψ ) , q h ) M ∀ q h ∈ Q kh ( M ) . Since ψ ∗ h = 0 on ∂M it is easy to check that the second equation also satisfies for a constant function.Here we remark that the existing stable pair that fits into the above spaces can be Taylor-Hood element(cf. [7]).Let ψ | M = − curl ψ ∗ h | M , then we can check that( as ( ψ − ψ ) , η h ) M = ( as ( ψ ) , η h ) ∀ η h ∈ M h ( M ) , k ψ k ,M ≤ C k ψ − ψ k ,M , which satisfies (3.10) by summing over all the primal elements. Here M h ( M ) denotes the space M h restricted to each primal element M ∈ T u . In addition, it also holds B ∗ h ( v h , ψ ) = 0 ∀ v h ∈ [ S h ] . Therefore, we have built ψ and ψ which satisfy (3.9)-(3.11). Hence, by the well-known argumentsfrom [5] we can infer that the existence of ψ and ψ leads to the inf-sup condition. Remark 3.2.
In the above proof, we employ Taylor-Hood stable pair, in fact we can also exploitScott-Vogelius element. In which case, we need to be careful about the singular nodes, indeed ifsingular nodes exists then one needs to enforce additional restriction on those nodes (cf. [21]), whichis naturally satisfied for our choice of space M h .We also emphasize that we choose locally H -conforming space for γ h to enforce the symmetry ofstress, which is different from that of [28]. In fact, we can also exploit fully discontinuous space for γ h , which then leads to strongly symmetric stress. In this case, the final system will be singular forrectangular mesh if one chooses the center point as ν . One can perturb the position of ν to avoidsingularity. We want to develop a numerical scheme which is stable regardless of the position of ν ,thus we choose locally H -conforming space for γ h .In order to prove the well-posedness of (3.1)-(3.5), we will make use of the following theorem; cf.[43]. Theorem 3.1.
Let the linear, symmetric and monotone operators N be given for the real vector space E to its dual algebraic dual E ∗ , and let E ′ b be the Hilbert space which is the dual of E with the seminorm | x | b = ( N x ( x )) , ∀ x ∈ E. Let
M ⊂ E × E ′ b be a relation with domain D = { x ∈ E : M ( x ) = 0 } . Assume M is monotone and Rg ( N + M ) = E ′ b . Then for each x ∈ D and for each F ∈ W , (0 , T ; E ′ b ) , there is a solution x of ddt ( N x ( t )) + M ( x ( t )) ∋ F ( t ) , ∀ < t < T (3.12) with N x ∈ W , ∞ (0 , T ; E ′ b ) , x ( t ) ∈ D, ∀ ≤ t ≤ T, and N x (0) = N x . Theorem 3.2.
There exists a unique solution to (3.1) - (3.5) . ocking free SDG for poroelasticity Proof.
In order to fit (3.1)-(3.5) in the form of Theorem 3.1, we consider a slightly modified formu-lation, with (3.1) differentiated in time and the new variables ˙ u h and ˙ γ h representing d t u h and d t γ h ,respectively: ( ∂ t A ( σ h + αp h I ) , ψ ) − B ∗ h ( ˙ u h , ψ ) + ( ˙ γ h , as ( ψ )) = 0 ∀ ψ ∈ [Σ h ] . (3.13)We now introduce the operators( A σσ σ h , ψ ) = ( Aσ h , ψ ) , ( A σp σ h , w ) = ( Aσ h , αwI ) , ( A σu σ h , v ) = B h ( σ h , v ) , ( A σγ σ h , η ) = ( as ( σ h ) , η ) , ( A zz z h , ξ ) = ( K − z h , ξ ) , ( A zp z h , w ) = b h ( z h , w ) , ( A pp p h , w ) = c ( p h , w ) + α ( Aαp h I, wI ) . We have a system in the form of (3.12), where˙ x = σ h ˙ u h ˙ γ h z h p h , N = A σσ A Tσp A σp A pp , M = − A Tσu A Tσγ A σu A σγ A zz A Tzp − A zp and F = (0 , f, , , q ) T .The dual space E ′ b is L (Ω) × × × × × L (Ω), and the condition F ∈ W , (0 , T ; E ′ b ) inTheorem 3.1 allows for non-zero source terms only in the equations with time derivatives, whichmeans f = 0. We can reduce our problem to a system with f = 0 by solving for each t ∈ (0 , T ] anelasticity problem with a source term f , cf. [42].( Aσ fh , ψ ) − B ∗ h ( ˙ u fh , ψ ) + ( ˙ γ fh , as ( ψ )) = 0 ∀ ψ ∈ [Σ h ] ,B h ( σ fh , v ) = ( f, v ) ∀ v ∈ [ S h ] , ( as ( σ fh ) , η ) = 0 ∀ η ∈ M h . Subtracting the above equations from (3.1)-(3.5), we can obtain( A (( σ h − ˙ σ fh ) + αp h I ) , ψ ) − B ∗ h ( u h − ˙ u fh , ψ ) + ( γ h − ˙ γ fh , as ( ψ )) = ( A ( σ fh − ∂ t σ fh ) , ψ ) , (3.14) B h ( σ h − σ fh , v ) = 0 , (3.15)( as ( σ h − σ fh ) , η ) = 0 , (3.16)( K − z h , ξ ) + b ∗ h ( p h , ξ ) = 0 , (3.17) c ( ∂ t p h , w ) + α ( ∂ t A (( σ h − σ fh ) + αp h I ) , wI ) − b h ( z h , w ) = ( q, w ) − ( α∂ t Aσ fh , wI ) , (3.18)which results in a problem with a modified right hand side F = ( A σσ ( σ fh − ∂ t σ fh ) , , , , q − A σp ∂ t σ fh ) T .The range condition Rg ( N + M ) = E ′ b can be verified by showing that the square finite dimensionalhomogeneous system: Find (ˆ σ h , ˆ u h , ˆ γ h , ˆ z h , ˆ p h ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h such that( A (ˆ σ h + α ˆ p h I ) , ψ ) − B ∗ h (ˆ u h , ψ ) + (ˆ γ h , as ( ψ )) = 0 ∀ ψ ∈ [Σ h ] , (3.19) B h (ˆ σ h , v ) = 0 ∀ v ∈ [ S h ] , (3.20)( as (ˆ σ h ) , η ) = 0 ∀ η ∈ M h , (3.21)( K − ˆ z h , ξ ) + b ∗ h (ˆ p h , ξ ) = 0 ∀ ξ ∈ Σ h , (3.22) c ( ∂ t ˆ p h , w ) + α ( ∂ t A (ˆ σ h + α ˆ p h I ) , wI ) − b h (ˆ z h , w ) = 0 ∀ w ∈ S h , (3.23)has only zero solution. Taking ( ψ, v, η, ξ, w ) = (ˆ σ h , ˆ u h , ˆ γ h , ˆ z h , ˆ p h ) in (3.19)-(3.23) and summing up theresulting equations yield k A (ˆ σ h + α ˆ p h I ) k + c k ˆ p h k + k K − ˆ z h k = 0 , ocking free SDG for poroelasticity σ h + α ˆ p h I = 0 and ˆ z h = 0. By the inf-sup condition (3.6), we have k ˆ p h k Z ≤ C k ˆ z h k , thereby ˆ p h = 0 by using the discrete Poincar´e-Friedrichs inequality (cf. [6]), consequently, ˆ σ h = 0.The inf-sup condition (3.8) and (3.19) lead to k ˆ u h k h + k ˆ γ h k ≤ C k A (ˆ σ h + α ˆ p h I ) k , which yields ˆ u h = 0 and ˆ γ h = 0.Finally, we need compatible initial data ˙ x ∈ D , i.e., M ˙ x ∈ E ′ b . Let us first consider the initial data x = ( σ h , u h , γ h , z h , p h ) for the discrete formulation (3.1)-(3.5). We take x to be the elliptic projectionof the initial data ˜ x = ( σ , u , γ , z , p ) for the weak formulation (2.5)-(2.9), which is constructedfrom p by the procedure described at the end of section 2. With the reduction to a problem with f = 0, the construction satisfies ( N + M )˜ x ∈ E ′ b . Then by( N + M ) x = ( N + M )˜ x , we have M x = ( N + M )˜ x − N x ∈ E ′ b . For the initial data of the differentiated problem (3.13),(3.2)-(3.5), we simply take ˙ x = ( σ h , , , z h , p h ), which also satisfies M ˙ x ∈ E ′ b . Here u h and γ h are not needed for the differentiated problem, but will be used to recover the solution of the originalproblem.Now all conditions of Theorem 3.1 are satisfied and we conclude the existence of a solution of (3.13),(3.2)-(3.5) with σ h (0) = σ h and p h (0) = p h . By taking t → z h and p h satisfy(3.4) at t = 0, we also infer that z h (0) = z h .Next, we recover the solution of the original problem. Let us define u h ( t ) = u h + Z t ˙ u h ds, γ h ( t ) = γ h + Z t ˙ γ h ds ∀ t ∈ (0 , T ] . By construction, u h (0) = u h and γ h (0) = γ h . Integrating (3.13) from 0 to any t ∈ (0 , T ] and using that σ h , u h , γ h satisfy (3.1) at t = 0, we conclude that (3.1) holds for all t . This completes the existenceproof. Uniqueness follows from the stability bound given in Theorem 4.1 in the next section.Before closing this section, we introduce some interpolation operators which will play an importantrole for later analysis. There exists the interpolation operator I h : H (Ω) → S h satisfying for any ω ∈ H (Ω) ( I h ω − ω, ψ ) e = 0 ∀ ψ ∈ P k ( e ) , e ∈ F u , ( I h ω − ω, ψ ) τ = 0 ∀ ψ ∈ P k − ( τ ) , τ ∈ T h . (3.24)Similarly, there also exists an operator Π h : [ H δ (Ω)] → Σ h , δ > / ϕ ∈ [ H δ (Ω)] ((Π h ϕ − ϕ ) · n, φ ) e = 0 ∀ φ ∈ P k ( e ) , e ∈ F p , (Π h ϕ − ϕ, φ ) τ = 0 ∀ φ ∈ P k − ( τ ) , τ ∈ T h . (3.25)The definition of the interpolation operators implies that b h ( z − Π h z, w ) = 0 ∀ w ∈ S h , (3.26) b ∗ h ( p − I h p, v ) = 0 ∀ v ∈ Σ h . (3.27)In addition for Π h σ = (Π h σ , Π h σ ) ∈ [Σ h ] and I h u = ( I h u , I h u ) ∈ [ S h ] , we also have B h ( σ − Π h σ, v ) = 0 ∀ v ∈ [ S h ] , (3.28) B ∗ h ( u − I h u, ψ ) = 0 ∀ ψ ∈ [Σ h ] . (3.29) ocking free SDG for poroelasticity I h and Π h are polynomial preserving operators, which satisfy the followinginterpolation error estimates (cf. [17, 13]) k ϕ − Π h ϕ k ≤ Ch k +1 k ϕ k H k +1 (Ω) ∀ ϕ ∈ H k +1 (Ω) , k ω − I h ω k ≤ Ch k +1 k ω k H k +1 (Ω) ∀ ω ∈ H k +1 (Ω) . (3.30)In addition, we take π h to be the standard nodal interpolation operator, which satisfies k π h ζ − ζ k ≤ Ch k +1 k ζ k H k +1 (Ω) ∀ ζ ∈ H k +1 (Ω) . (3.31) In this section, we prove the convergence estimates for the semi-discrete scheme and fully discretescheme with backward Euler in time. In particular we will show that the stability and convergenceerror estimates are independent of the anisotropy ratio ̺ B and of the storativity coefficient c , and arevalid for c = 0. In this subsection we prove the stability and convergence estimates for the semi-discrete scheme. Tobegin, we consider the following stability estimate.
Theorem 4.1. (stability). Let ( σ h , u h , γ h , z h , p h ) be the numerical solution of (3.1) - (3.5) , then wehave k σ h k L ∞ (0 ,T ; L (Ω)) + k p h k L ∞ (0 ,T ; L (Ω)) + k u h k L ∞ (0 ,T ; L (Ω)) + k K − z h k L ∞ (0 ,T ; L (Ω)) + k γ h k L ∞ (0 ,T ; L (Ω)) + k K − z h k L (0 ,T ; L (Ω)) + k u h k L (0 ,T ; L (Ω)) + k p h k L (0 ,T ; L (Ω)) + k σ h k L (0 ,T ; L (Ω)) + k γ h k L (0 ,T ; L (Ω)) ≤ C (cid:16) k f k L ∞ (0 ,T ; L (Ω)) + k f k H (0 ,T ; L (Ω)) + k q k H (0 ,T ; L (Ω)) + k q k L ∞ (0 ,T ; L (Ω)) + k p k H (Ω) (cid:17) . Proof.
Differentiating (3.1) in time, taking ψ = σ h , v = ∂ t u h , η = ∂ t γ h , ξ = z h , w = p h in (3.1)-(3.5)and summing up the resulting equations yield12 ddt (cid:16) k A ( σ h + αp h I ) k + c k p h k (cid:17) + k K − z h k = ( f, ∂ t u h ) + ( q, p h ) . Integrating over time implies12 (cid:16) k A ( σ h + αp h I )( t ) k + c k p h ( t ) k (cid:17) + Z t k K − z h k ds = − Z t ( ∂ t f, u h ) ds + Z t ( q, p h ) ds + 12 (cid:16) k A ( σ h + αp h I )(0) k + c k p h k (cid:17) + ( f, u h )( t ) − ( f, u h )(0) . The Cauchy-Schwarz inequality and Young’s inequality lead to k A ( σ h + αp h I )( t ) k + c k p h ( t ) k + 2 Z t k K − z h k ds ≤ ǫ (cid:16) k u h ( t ) k + Z t ( k p h k + k u h k ) ds (cid:17) + 12 ǫ (cid:16) k f ( t ) k + Z t ( k q k + k ∂ t f k ) ds (cid:17) + k A ( σ h + αp h I )(0) k + c k p h k + k u h k + k f (0) k . (4.1)By (3.1) and the inf-sup condition (3.8), we have k u h k h + k γ h k ≤ C sup ψ ∈ [Σ h ] − B ∗ h ( u h , ψ ) + ( as ( ψ ) , γ h ) k ψ k ≤ C k A ( σ h + αp h I ) k . ocking free SDG for poroelasticity k u h k ≤ C k A ( σ h + αp h I ) k . (4.2)Taking ψ = σ h , v = u h and η = γ h in (3.1)-(3.3), and summing up the resulting equations, we deducethat k σ h k ≤ C ( k p h k + ǫ k u h k + 1 ǫ k f k ) . (4.3)Thus, we can infer from (4.2) and (4.3) that k u h k ≤ C k A ( σ h + αp h I ) k ≤ C ( k σ h k + k p h k ) ≤ C (cid:16) k p h k + ǫ k u h k + 1 ǫ k f k (cid:17) . Choosing ǫ small enough yields k u h k ≤ C ( k p h k + 1 ǫ k f k ) . (4.4)On the other hand, we have from (3.4) and the inf-sup condition (3.6) k p h k Z ≤ C sup ξ ∈ Σ h b ∗ h ( p h , ξ ) k ξ k = C sup ξ ∈ Σ h ( K − z h , ξ ) k ξ k ≤ C k K − z h k . (4.5)Combining the discrete Poincar´e-Friedrichs inequalities, (4.4) and (4.5) lead to Z t ( k u h k + k p h k ) ds ≤ C Z t ( k K − z h k + k f k ) ds. (4.6)Therefore, choosing ǫ small enough, we can infer from (4.1), (4.5) and (4.6) that k A ( σ h + αp h I )( t ) k + c k p h ( t ) k + k u h ( t ) k + k γ h ( t ) k + Z t (cid:16) k K − z h k + k u h k + k p h k + k σ h k + k γ h k (cid:17) ds ≤ C (cid:16) k f ( t ) k + Z t ( k f k + k q k + k ∂ t f k ) ds + k p h k + k u h k + k f (0) k + k σ h k (cid:17) . (4.7)Next, differentiating (3.1)-(3.4) in time and combining them with (3.5) with the choices ( ψ, v, η, ξ, w ) =( ∂ t σ h , ∂ t u h , ∂ t γ h , z h , ∂ t p h ), we can obtain2 Z t (cid:16) k ∂ t A ( σ h + αp h I ) k + c k ∂ t p h k (cid:17) ds + k K − z h ( t ) k ≤ ǫ (cid:16) k p h ( t ) k + Z t k ∂ t u h k (cid:17) + 1 ǫ (cid:16) k q ( t ) k + Z t k ∂ t f k ds (cid:17) + Z t ( k p h k + k ∂ t q k ) ds + k K − z h k + k p h k + k q (0) k . (4.8)Then by inf-sup condition (3.8) and (3.1) differentiated in time, we have k ∂ t u h k h + k ∂ t γ h k ≤ C k ∂ t A ( σ h + αp h I ) k . (4.9)Choosing ǫ in (4.8) small enough yields Z t (cid:16) k ∂ t A ( σ h + αp h I ) k + c k ∂ t p h k (cid:17) ds + k K − z h ( t ) k + k p h ( t ) k ≤ C (cid:16) Z t ( k p h k + k ∂ t q k + k ∂ t f k ) ds + k q ( t ) k + k K − z h k + k p h k + k q (0) k (cid:17) , ocking free SDG for poroelasticity k σ h ( t ) k + k p h ( t ) k + k u h ( t ) k + k K − z h ( t ) k + k γ h ( t ) k + Z t (cid:16) k K − z h k + k u h k + k p h k + k σ h k + k γ h k (cid:17) ds ≤ C (cid:16) k f ( t ) k + Z t (cid:16) k f k + k q k + k ∂ t f k + k ∂ t q k (cid:17) ds + k q ( t ) k + k p h k + k u h k + k f (0) k + k σ h k + k K − z h k + k q (0) k (cid:17) . Since the discrete initial data ( σ h , u h , γ h , z h , p h ) is the elliptic projection of the initial data ( σ , u , γ , z , p ),we have k σ h k + k u h k + k γ h k + k z h k + k p h k ≤ C ( k p k H (Ω) + k f (0) k ) . (4.10)Therefore, the proof is completed.Next, we will prove the convergence estimate. To simplify the notation, we let σ − σ h = ( σ − Π h σ ) + (Π h σ − σ h ) := C σ + D σ ,u − u h = ( u − I h u ) + ( I h u − u h ) := C u + D u γ − γ h = ( γ − π h γ ) + ( π h γ − γ h ) := C γ + D γ ,z − z h = ( z − Π h z ) + (Π h z − z h ) := C z + D z ,p − p h = ( p − I h p ) + ( I h p − p h ) := C p + D p . Theorem 4.2.
Let ( σ h , u h , γ h , z h , p h ) be the numerical solution of (3.1) - (3.5) and assume that thesolution of (2.5) - (2.9) is sufficiently smooth, then we have the following error estimate k σ − σ h k L ∞ (0 ,T ; L (Ω)) + k p − p h k L ∞ (0 ,T ; L (Ω)) + k u − u h k L ∞ (0 ,T ; L (Ω)) + k γ − γ h k L ∞ (0 ,T ; L (Ω)) + k K − ( z − z h ) k L ∞ (0 ,T ; L (Ω)) + k σ − σ h k L (0 ,T ; L (Ω)) + k p − p h k L (0 ,T ; L (Ω)) + k u − u h k L (0 ,T ; L (Ω)) + k σ − σ h k L (0 ,T ; L (Ω)) + k K − ( z − z h ) k L (0 ,T ; L (Ω)) ≤ Ch k +1) (cid:16) k u k L (0 ,T ; H k +1 (Ω)) + k γ k H (0 ,T ; H k +1 (Ω)) + k p k H (0 ,T ; H k +1 (Ω)) + k σ k H (0 ,T ; H k +1 (Ω)) + k K − z k H (0 ,T ; H k +1 (Ω)) + k γ k L ∞ (0 ,T ; H k +1 (Ω)) + k K − z k L ∞ (0 ,T ; H k +1 (Ω)) + k σ k L ∞ (0 ,T ; H k +1 (Ω)) + k p k L ∞ (0 ,T ; H k +1 (Ω)) + k u k L ∞ (0 ,T ; H k +1 (Ω)) (cid:17) . Proof.
We obtain the following error equations by replacing ( σ, u, γ, z, p ) by ( σ h , u h , γ h , z h , p h ) andperforming integration by parts on (3.1)-(3.5)( A ( D σ + αD p I ) , ψ ) − B ∗ h ( D u , ψ ) + ( D γ , as ( ψ )) = − ( A ( C σ + αC p I ) , ψ )+ B ∗ h ( C u , ψ ) − ( C γ , as ( ψ )) , (4.11) B h ( D σ , v ) = − B h ( C σ , v ) , (4.12)( as ( D σ ) , η ) = − ( as ( C σ ) , η ) , (4.13)( K − D z , ξ ) + b ∗ h ( D p , ξ ) = − ( K − C z , ξ ) − b ∗ h ( C p , ξ ) , (4.14) c ( ∂ t D p , w ) + α ( ∂ t A ( D σ + αD p I ) , wI ) − b h ( D z , w ) = − c ( ∂ t C p , w ) − α ( ∂ t A ( C σ + αD p I ) , wI ) + b h ( C z , w ) (4.15)for ( ψ, v, η, ξ, w ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h .Differentiating (4.11) with time, taking ψ = D σ , v = ∂ t D u , η = ∂ t D γ , ξ = D z , w = D p , and summingup the resulting equations yield12 ∂ t (cid:16) k A ( D σ + αD p I ) k + c k D p k (cid:17) + k K − D z k = − ( ∂ t A ( C σ + αC p I ) , D σ )+ B ∗ h ( ∂ t C u , D σ ) − ( ∂ t C γ , as ( D σ )) − B h ( C σ , ∂ t D u ) + ( as ( C σ ) , ∂ t D γ ) − ( K − C z , D z ) − b ∗ h ( C p , D z ) − c ( ∂ t C p , D p ) − α ( ∂ t A ( C σ + αC p I ) , D p I ) + b h ( C z , D p ) . (4.16) ocking free SDG for poroelasticity B ∗ h ( ∂ t C u , D σ ) = 0 , − B h ( C σ , ∂ t D u ) = 0 , b ∗ h ( C p , D z ) = 0 , b h ( C z , D p ) = 0 . In addition, the Cauchy-Schwarz inequality yields | − ( ∂ t A ( C σ + αC p I ) , D σ ) | ≤ k ∂ t A ( C σ + αC p I k k A D σ k , | ( ∂ t C γ , as ( D σ )) | ≤ k ∂ t C γ k k as ( D σ ) k , | ( K − C z , D z ) | ≤ k K − C z k k K − D z k , | c ( ∂ t C p , D p ) | ≤ c k ∂ t C p k k D p k , | ( ∂ t A ( C σ + αC p I ) , D p I ) | ≤ k ∂ t A ( C σ + αC p I ) k k D p k . Therefore, we have from (4.16) that12 ∂ t (cid:16) k A ( D σ + D p I ) k + c k D p k (cid:17) + k K − D z k ≤ k ∂ t A ( C σ + αC p I ) k k A D σ k + k ∂ t C γ k k D σ k + ( as ( C σ ) , ∂ t D γ ) + k K − C z k k K − D z k + c k ∂ t C p k k D p k + k ∂ t A ( C σ + αC p I ) k k D p k . Integrating in time from 0 to an arbitrary t ∈ (0 , T ] results in12 (cid:16) k A ( D σ + D p I )( t ) k + c k D p ( t ) k (cid:17) + Z t k K − D z k ds ≤ (cid:16) k A ( D σ + D p I )(0) k + c k D p (0) k (cid:17) − Z t ( ∂ t as ( C σ ) , D γ ) ds + 1 ǫ k C σ ( t ) k + ǫ k D γ ( t ) k − ( as ( C σ ) , D γ )(0)+ ǫ Z t (cid:16) k D p k + k K − D z k + k D σ k (cid:17) ds + 1 ǫ Z t (cid:16) k ∂ t C γ k + k K − C z k + k ∂ t C p k + k ∂ t A ( C σ + αC p I ) k (cid:17) ds. (4.17)It follows from the inf-sup condition (3.8) and (4.11) that k D u k h + k D γ k ≤ C sup ψ ∈ [Σ h ] − B ∗ h ( D u , ψ ) + ( as ( ψ ) , D γ ) k ψ k = C sup ψ ∈ [Σ h ] − ( A ( D σ + αD p I ) , ψ ) − ( A ( C σ + αC p I ) , ψ ) − ( C γ , as ( ψ )) k ψ k ≤ C (cid:16) k A ( D σ + αD p I ) k + k A ( C σ + αC p I ) k + k C γ k (cid:17) . (4.18)Taking ψ = D σ , v = D u , η = D γ in (4.11)-(4.13), and summing up the resulting equations yield( A ( D σ + αD p I ) , D σ ) = − ( A ( C σ + αC p I ) , D σ ) − ( C γ , as ( D σ )) + ( as ( C σ ) , D γ ) , which together with Young’s inequality implies k A D σ k ≤ C (cid:16) k D p k + k C σ k + k C p k + k C γ k + 1 ǫ k C σ k + ǫ k D γ k (cid:17) . Taking ǫ small enough and using (4.18), we can obtain k D σ k ≤ C (cid:16) k C p k + k D p k + k C σ k + k C γ k (cid:17) . (4.19)It follows from (3.6) and (4.14) that k D p k Z ≤ C sup ξ ∈ Σ h b ∗ h ( D p , ξ ) k ξ k = C sup ξ ∈ Σ h − ( K − D z , ξ ) − ( K − C z , ξ ) k ξ k ≤ C ( k K − D z k + k K − C z k ) , (4.20) ocking free SDG for poroelasticity Z t k D p k ds ≤ Z t k D p k Z ds ≤ Z t ( k K − D z k + k K − C z k ) ds. An application of (4.19) and (4.20) leads to Z t k D σ k ds ≤ C Z t (cid:16) k K − D z k + k K − C z k + k C σ k + k C γ k + k C p k (cid:17) ds. Besides, we can infer from (4.18) and (4.19) that k D u k h + k D γ k ≤ C (cid:16) k K − D z k + k K − C z k + k C σ k + k C γ k + k C p k (cid:17) . Choosing ǫ and ǫ in (4.17) small enough, we can infer from the preceding arguments k A ( D σ + D p I )( t ) k + c k D p ( t ) k + Z t (cid:16) k K − D z k + k D σ k + k D p k + k D γ k + k D u k (cid:17) ds ≤ (cid:16) k A ( D σ + D p I )(0) k + c k D p (0) k (cid:17) + k C σ (0) k + k D γ (0) k + k C σ ( t ) k + Z t (cid:16) k C γ k + k K − C z k + k ∂ t C p k + k C σ k + k ∂ t C σ k + k C p k + k ∂ t C γ k (cid:17) ds. (4.21)Differentiating (4.11)-(4.14) in time and setting ( ψ, v, η, ξ, w ) = ( ∂ t D σ , ∂ t D u , ∂ t D γ , D z , ∂ t D p ), we canobtain k ∂ t A ( D σ + αD p I ) k + c k ∂ t D p k + 12 ∂ t k K − D z k = − ( ∂ t A ( C σ + αC p I ) , ∂ t D σ + α∂ t D p I ) − ( ∂ t C γ , as ( ∂ t D σ )) + ( ∂ t as ( C σ ) , ∂ t D γ ) − ( ∂ t K − C z , D z ) − c ( ∂ t C p , ∂ t D p ) . Integrating over time and using as ( ∂ t D σ ) = as ( ∂ t ( D σ + αD p I )) yield Z t k ∂ t A ( D σ + αD p I ) k ds + c Z t k ∂ t D p k ds + 12 k K − D z ( t ) k ≤ ǫ Z t k ∂ t A ( D σ + αD p I ) k ds + 14 ǫ (cid:16) Z t k ∂ t A ( C σ + αC p I ) k ds + Z t k ∂ t C γ k ds (cid:17) + 14 ǫ Z t k ∂ t C σ k ds + ǫ Z t k ∂ t D γ k ds + Z t k K − ∂ t C z k ds + Z t k K − D z k ds + c Z t k ∂ t C p k ds + c Z t k ∂ t D p k ds + 12 k K − D z (0) k . (4.22)By inf-sup condition (3.8) and (4.11) differentiated in time, we can obtain k ∂ t D u k h + k ∂ t D γ k ≤ C k ∂ t A ( D σ + αD p I ) k . Choosing ǫ and ǫ in (4.22) small enough leads to Z t k ∂ t A ( D σ + αD p I ) k ds + c Z t k ∂ t D p k ds + k K − D z ( t ) k ≤ C (cid:16) Z t k ∂ t A ( C σ + αC p I ) k ds + Z t k ∂ t C γ k ds + Z t k ∂ t C σ k ds + Z t k K − ∂ t C z k ds + c Z t k ∂ t C p k ds + Z t k K − D z k ds + k K − D z (0) k (cid:17) , ocking free SDG for poroelasticity k D σ ( t ) k + k D p ( t ) k + k D u ( t ) k + k D γ ( t ) k + k K − D z ( t ) k + Z t (cid:16) k K − D z k + k D σ k + k D p k + k D γ k + k D u k (cid:17) ds ≤ C (cid:16) k A ( D σ + D p I )(0) k + c k D p (0) k + Z t (cid:16) k C γ k + k C z k + k C σ k + k C p k + k ∂ t C γ k + k ∂ t C p k + k ∂ t A ( C σ + αC p I ) k + k ∂ t C σ k (cid:17) ds + k K − D z (0) k + k C σ (0) k + k D γ (0) k + k K − ∂ t C z ( t ) k + k C σ ( t ) k + k C p ( t ) k + k C γ ( t ) k (cid:17) . The initial data satisfies k D σ (0) k + k D γ (0) k + k D p (0) k + k K − D z (0) k ≤ C (cid:16) k C σ (0) k + k C γ (0) k + k C p (0) k + k K − C z (0) k (cid:17) . (4.23)Therefore, the proof is completed by combining the preceding arguments and the interpolation errorestimates (3.30)-(3.31). Remark 4.1. (robustness of error estimates for nearly incompressible materials). We can observefrom the above proof that we employ the coercivity of A , while the coercivity of A on L (Ω) × is notuniform in λ . Indeed, c λ k ψ k ≤ k A ψ k ∀ ψ ∈ L (Ω) × holds with a constant c λ > c λ → λ → + ∞ . It means that our error estimates obtainedusing the coercivity of A may have error bounds growing unboundedly as λ → + ∞ . In order to get auniform bound, we can follow the framework given in [28, 53], which is omitted for simplicity. In fact,proceeding analogously to [28, 53] we can check that our proposed method has uniformly boundedestimate. In this subsection we analyze the convergence estimates for the fully discrete scheme. To this end weintroduce a partition of the time interval [0 , T ] into subintervals [ t n − , t n ] , ≤ n ≤ N ( N is an integer)and denote the time step size by ∆ t = TN . Using backward Euler scheme in time, we get the fullydiscrete staggered DG method as follows: Find ( σ h , u h , γ h , z h , p h ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h suchthat ( A ( σ nh + αp nh I ) , ψ ) − B ∗ h ( u nh , ψ ) + ( γ nh , as ( ψ )) = 0 , (4.24) B h ( σ nh , v ) = ( f n , v ) , (4.25)( as ( σ nh ) , η ) = 0 , (4.26)( K − z nh , ξ ) + b ∗ h ( p nh , ξ ) = 0 , (4.27) c ( p nh − p n − h ∆ t , w ) + α ( A ( σ nh + αp nh I ) − A ( σ n − h + αp n − h I )∆ t , wI ) − b h ( z nh , w ) = ( q n , w ) (4.28)for ( ψ, v, η, ξ, w ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h .For a Sobolev space H on Ω with a norm k · k H , we define the discrete in time norms by k ψ k l (0 ,T ; H ) := (cid:16) N X n =1 ∆ t k ψ k H (cid:17) / , k ψ k l ∞ (0 ,T ; H ) := max ≤ n ≤ N k ψ k H ocking free SDG for poroelasticity σ n − σ nh = ( σ n − Π h σ n ) + (Π h σ n − σ nh ) := C nσ + D nσ ,u n − u nh = ( u n − I h u n ) + ( I h u n − u nh ) := C nu + D nu γ n − γ nh = ( γ n − π h γ n ) + ( π h γ n − γ nh ) := C nγ + D nγ ,z n − z nh = ( z n − Π h z n ) + (Π h z n − z nh ) := C nz + D nz ,p n − p nh = ( p n − I h p n ) + ( I h p n − p nh ) := C np + D np . Theorem 4.3. (existence and uniqueness). There exists a unique solution to the fully discrete scheme (4.24) - (4.28) .Proof. Since (4.24)-(4.24) is a square linear system, uniqueness implies existence. It suffices to showthe uniqueness. To this end, suppose f n , q n and the ( n − A ( σ nh + αp nh I ) , ψ ) − B ∗ h ( u nh , ψ ) + ( γ nh , as ( ψ )) = 0 ∀ ψ ∈ [Σ h ] , (4.29) B h ( σ nh , v ) = 0 ∀ v ∈ [ S h ] , (4.30)( as ( σ nh ) , η ) = 0 ∀ η ∈ M h , (4.31)( K − z nh , ξ ) + b ∗ h ( p nh , ξ ) = 0 ∀ ξ ∈ Σ h , (4.32) c ( p nh , w ) + α ( A ( σ nh + αp nh I ) , wI ) − ∆ tb h ( z nh , w ) = 0 ∀ w ∈ S h . (4.33)Taking ψ = σ nh , v = u nh , η = γ nh , ξ = z nh , w = p nh in (4.29)-(4.33) and summing up the resulting equationsyield k A ( σ nh + αp nh I ) k + c k p nh k + ∆ t k K − z nh k = 0 , thereby σ nh + αp nh I = 0 and z nh = 0.We can infer from inf-sup condition (3.6) and (4.32) that k p nh k Z ≤ C k K − z nh k , which implies p nh = 0. Hence σ nh = 0 because of σ nh + αp nh I = 0.It follows from (3.8) and (4.29) that k u nh k h + k γ nh k ≤ C k A ( σ nh + αp nh I ) k , which gives u nh = 0 and γ nh = 0. This completes the proof. Theorem 4.4.
Let ( σ h , u h , γ h , z h , p h ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h be the numerical solution of thefully discrete scheme (4.24) - (4.28) and assume that the solution of (2.5) - (2.9) is sufficiently smooth,then we have the following convergence estimate k K − ( z − z h ) k l (0 ,T ; L (Ω)) + k σ − σ h k l (0 ,T ; L (Ω)) + k z − z h k l (0 ,T ; L (Ω)) + k u − u h k l (0 ,T ; L (Ω)) + k σ − σ h k l (0 ,T ; L (Ω)) + k p − p h k l ∞ (0 ,T ; L (Ω)) + k K − ( z − z h ) k l ∞ (0 ,T ; L (Ω)) + k σ − σ h k l ∞ (0 ,T ; L (Ω)) + k γ − γ h k l ∞ (0 ,T ; L (Ω)) + k u − u h k l ∞ (0 ,T ; L (Ω)) ≤ C (cid:16) h k +1) (cid:16) k p k H (0 ,T ; H k +1 (Ω)) + k σ k H (0 ,T ; H k +1 (Ω)) + k K − z k H (0 ,T ; H k +1 (Ω)) + k γ k H (0 ,T ; H k +1 (Ω)) + k σ k L ∞ (0 ,T ; H k +1 (Ω)) + k p k L ∞ (0 ,T ; H k +1 (Ω)) + k u k L (0 ,T ; H k +1 (Ω)) + k u k L ∞ (0 ,T ; H k +1 (Ω)) + k γ k L ∞ (0 ,T ; H k +1 (Ω)) + k K − z k L ∞ (0 ,T ; H k +1 (Ω)) (cid:17) + ∆ t (cid:16) k p k H (0 ,T ; L (Ω)) + k σ k H (0 ,T ; L (Ω)) (cid:17)(cid:17) . ocking free SDG for poroelasticity Proof.
We have the following error equations by performing integration by parts on (4.24)-(4.28)( A ( σ n − σ nh + α ( p n − p nh ) I ) , ψ ) − B ∗ h ( u n − u nh , ψ ) + ( γ n − γ nh , as ( ψ )) = 0 , (4.34) B h ( σ n − σ nh , v ) = 0 , (4.35)( as ( σ n − σ nh ) , η ) = 0 , (4.36)( K − ( z n − z nh ) , ξ ) + b ∗ h ( p n − p nh , ξ ) = 0 , (4.37) α ∆ t ( A ( σ n − σ nh + αp n − p nh I ) − A ( σ n − − σ n − h + α ( p n − − p n − h I )) , wI )+ c ∆ t ( p n − p nh − ( p n − − p n − h ) , w ) − b h ( z n − z nh , w ) = c ( p n − p n − ∆ t − p t (: , t n ) , w )+ α ∆ t ( A ( σ n − σ n − + α ( p n − p n − ) I ) − ∆ tA ( σ t (: , t n ) + αp t (: , t n )) , wI ) (4.38)for ( ψ, v, η, ξ, w ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h .Consider (4.34) at iterations n and n − ψ = D nσ , v = D nu , η = D nγ − D n − γ , ξ = D nz , w = D np in the difference equation and (4.35)-(4.38), and summing upthe resulting equations, we have( K − ( z n − z nh ) , D nz ) + c ∆ t ( p n − p nh − ( p n − − p n − h ) , D np )+ 1∆ t ( A ( σ n − σ nh + α ( p n − p nh ) I ) − A ( σ n − − σ n − h + α ( p n − − p n − h ) I ) , D nσ + αD np I )= c ∆ t ( p n − p n − − ∆ tp t (: , t n ) , D np ) + 1∆ t ( A ( σ n − σ n − + α ( p n − p n − ) I ) − ∆ tA ( σ t (: , t n ) + αp t (: , t n )) , αD np I ) − t ( C nγ − C n − γ , as ( D nσ )) + 1∆ t ( D nγ − D n − γ , as ( C nσ )) . Using the inequality ( a − b, a ) ≥ | a | −| b | and the Cauchy-Schwarz inequality, we can obtain k K − D nz k + c t (cid:16) k D np k − k D n − p k + k D np − D n − p k (cid:17) + 12∆ t (cid:16) k A ( D nσ + αD np ) k − k A ( D n − σ + αD n − p I ) k (cid:17) ≤ C (cid:16) k K − C nz k k K − D nz k + c ∆ t k C np − C n − p kk D np k + 1∆ t ( D nγ − D n − γ , as ( C nσ ))+ 1∆ t k A ( C nσ + αC np I ) − A ( C n − σ + αC n − p I ) k (cid:16) k D nσ k + k D np k (cid:17) + 1∆ t k A ( σ n − σ n − + α ( p n − p n − ) I ) − ∆ tA ( σ t (: , t n ) + αp t (: , t n ) I ) k k D np k + 1∆ t k C nγ − C n − γ k k D nσ k + c ∆ t k p n − p n − − ∆ tp t (: , t n ) k k D np k (cid:17) . (4.39)It follows from (3.6), (4.37) and the discrete Poincar´e-Friedrichs inequality that k D np k ≤ C k D np k Z ≤ C sup ξ ∈ Σ h b ∗ h ( p nh − I h p n , ξ ) k ξ k ≤ C k K − ( z n − z nh ) k . (4.40)Similarly, we can infer from (3.8) and (4.34) that k D nu k h + k D nγ k ≤ C ( k A ( σ n − σ nh + α ( p n − p nh ) I ) k + k C nγ k ) . (4.41)Besides, taking ( ψ, v, η ) = ( D nσ , D nu , D nγ ) in (4.34)-(4.36) and summing up the resulting equations yield( A ( σ n − σ nh + α ( p n − p nh ) I ) , D nσ ) + ( C nγ , as ( D nσ )) − ( as ( C nσ ) , D nγ ) = 0 , then an application of Young’s inequality yields k D nσ k ≤ C (cid:16) k C np k + k D np k + k C nσ k + k C nγ k + ǫ k D nγ k + 1 ǫ k C nσ k (cid:17) . (4.42) ocking free SDG for poroelasticity ǫ small enough and combining with (4.40), (4.41) give k D nσ k ≤ C (cid:16) k C np k + k D np k + k C nσ k + k C nγ k (cid:17) ≤ C (cid:16) k p n − I h p n k + k K − ( z n − z nh ) k + k Π h σ n − σ n k + k γ n − π h γ n k (cid:17) . (4.43)On the other hand, we have from Taylor’s expansion p n − p n − − ∆ tp t (: , t n ) = − Z t n t n − ( t − t n − ) p tt (: , t n ) ds, where the right hand side can be estimated by the Cauchy-Schwarz inequality k t Z t n t n − ( t − t n − ) p tt (: , t n ) ds k ≤ ∆ t Z t n t n − k p tt k ds. Similarly, we have k t A ( σ n − σ n − − ∆ tσ t (: , t n )) k ≤ C ∆ t Z t n t n − k σ tt k ds, k t A (( p n − p n − ) I − ∆ tp t (: , t n ) I ) k ≤ C ∆ t Z t n t n − k p tt k ds. Therefore, we can infer from (4.39) and Young’s inequalities that k K − D nz k + c t (cid:16) k D np k − k D n − p k + k D np − D n − p k (cid:17) + 12∆ t (cid:16) k A ( D nσ + αD np I ) k − k A ( D n − σ + αD n − p I ) k (cid:17) ≤ C (cid:16) k K − C nz k + k C nσ k + k C nγ k + k C np k + 1(∆ t ) ( k C np − C n − p k + k A ( C nσ − C n − σ ) k + k C nγ − C n − γ k )+ ∆ t Z t n t n − (cid:16) k p tt k + k σ tt k (cid:17) ds + 1∆ t ( D nγ − D n − γ , as ( C nσ )) (cid:17) . (4.44)Changing n to j in (4.39) and summing over j = 1 , · · · , n , we can obtain2∆ t n X j =1 k K − D jz k + c (cid:16) k D np k − k D p k + n X j =1 k D np − D n − p k (cid:17) + k A ( D nσ + αD np I ) k ≤ C (cid:16) n X j =1 ∆ t ( k K − C jz k + k C jσ k + k C jγ k + k C jp k ) + n X j =1 t (cid:16) k C jp − C j − p k + k C jσ − C j − σ k + k C jγ − C j − γ k (cid:17) + ∆ t n X j =1 Z t j t j − ( k p tt k + k σ tt k ) ds + k A ( D σ + αD p I ) k + n X j =1 ( D jγ − D j − γ , as ( C jσ )) (cid:17) . (4.45)It follows from the interpolation error estimates (3.30)-(3.31) that k C jz k ≤ Ch k +1 k z j k H k +1 (Ω) , k C jγ k ≤ Ch k +1 k γ j k H k +1 (Ω) , k C jσ k ≤ Ch k +1 k σ j k H k +1 (Ω) , k C jp k ≤ Ch k +1 k p j k H k +1 (Ω) . The Cauchy-Schwarz inequality and the interpolation error estimates (3.30) imply1∆ t k C jp − C j − p k = 1∆ t k ( p j − I h p j ) − ( p j − − I h p j − ) k = 1∆ t k Z t j t j − ( p t − I h p t ) k ≤ Ch k +1) Z t j t j − k p t k H k +1 (Ω) , t k C jγ − C j − γ k ≤ Ch k +1) Z t j t j − k γ t k H k +1 (Ω) . (4.46) ocking free SDG for poroelasticity n X j =1 ( D jγ − D j − γ , as ( C jσ )) = D nγ C nσ − D γ C σ − n X j =1 ( D j − γ , C jσ − C j − σ ) , where the last term can be estimated by k C jσ − C j − σ k ≤ C ∆ th k +1) Z t j t j − k σ t k H k +1 (Ω) . Thus, the preceding arguments lead to∆ t n X j =1 k K − D jz k + c (cid:16) k D np k − k D p k + n X j =1 k D np − D n − p k (cid:17) + k A ( D nσ + αD np I ) k ≤ C (cid:16) h k +1) ( k γ k H (0 ,T ; H k +1 (Ω)) + k p k H (0 ,T ; H k +1 (Ω)) + k z k L (0 ,T ; H k +1 (Ω)) + k σ k H (0 ,T ; H k +1 (Ω)) + k σ k L ∞ (0 ,T ; H k +1 (Ω)) ) + ∆ t n X j =1 Z t j t j − ( k p tt k + k σ tt k ) ds + k A ( D σ + αD p I ) k + k C σ k + k D γ k (cid:17) . Considering (4.34)-(4.37) at iterations n and n −
1, and taking the difference, setting ψ = D nσ − D n − σ , v = D nu − D n − u , η = D nγ − D n − γ , ξ = D nz , w = D np − D n − p and summing up the resultingequations yield1∆ t k A ( D nσ + αD np I ) − A ( D n − σ + αD n − p I ) k + ( K − ( D nz − D n − z ) , D nz ) + c ∆ t ( p n − p nh − ( p n − − p n − h ) , D np − D n − p )= c ( p n − p n − ∆ t − p t (: , t n ) , D np − D n − p ) − ( K − ( C nz − C n − z ) , D nz )+ 1∆ t ( A ( σ n − σ n − + α ( p n − p n − ) I ) − ∆ tA ( σ t (: , t n ) + αp t (: , t n ) I ) , α ( D np − D n − p ) I ) − t ( A ( C nσ + αC np I ) − A ( C n − σ + αC n − p I ) , D nσ + αD np I − ( D n − σ + αD n − p I )) − t ( C nγ − C n − γ ) , as ( D nσ − D n − σ )) + 1∆ t ( as ( C nσ − C n − σ ) , D nγ − D n − γ ) . We can infer from inf-sup condition (3.6) and (3.8) that k I h p n − p nh − ( I h p n − − p n − h ) k ≤ C k K − ( z n − z nh − ( z n − − z n − h )) k , k π h γ n − γ nh − ( π h γ n − − γ n − h ) k ≤ C k A ( D nσ + αD np I ) − A ( D n − σ + αD n − p I ) k . Therefore, we have1∆ t k A ( D nσ + αD np I ) − A ( D n − σ + αD n − p I ) k + 12 (cid:16) k K − D nz k − k K − D n − z k + k K − D nz − K − D n − z k (cid:17) + c ∆ t k D np − D n − p k ≤ C (cid:16) ∆ t Z t n t n − ( k p tt k + k σ tt k ) ds + 1∆ t (cid:16) k C nσ − C n − σ k + k C np − C n − p k + k C nγ − C n − γ k + k K − ( C nz − C n − z ) k (cid:17) + ∆ t k D nz k (cid:17) . ocking free SDG for poroelasticity n to j and summing over j = 1 , · · · , n , we can obtain n X j =1 t k A ( D jσ + αD jp I ) − A ( D j − σ + αD j − p I ) k + k K − D nz k + n X j =1 c ∆ t k D jp − D j − p k ≤ C (cid:16) n X j =1 t (cid:16) k C jσ − C j − σ k + k C jp − C j − p k + k C jγ − C j − γ k + k K − ( C jz − C j − z ) k (cid:17) + k K − ( I h z − z h ) k + ∆ t Z t n ( k p tt k + k σ tt k ) ds + n X j =1 ∆ t k K − D jz k (cid:17) . Proceeding analogously to (4.46), we can get1∆ t k C jσ − C j − σ k ≤ Ch k +1) Z t j t j − k σ t k H k +1 (Ω) , t k C jz − C j − z k ≤ Ch k +1) Z t j t j − k K − z t k H k +1 (Ω) . Combining the above estimates with the interpolation error estimates (3.30)-(3.31) completes the proof.
The global system (4.24)-(4.28) consists of five variables which is relatively large. In order to reducethe computational costs, we propose the following fixed stress splitting scheme inspired by [30]. Thisincludes two steps: First, the flow problem is solved independently. Second, the mechanics problemis solved using updated pressure and flux. For fixed n, i ∈ N , the detailed splitting scheme reads asfollows: Step 1 : Given ( z n,i − h , p n,i − h , σ n,i − h ), find ( z n,ih , p n,ih ) such that( K − z n,ih , ξ ) + b ∗ h ( p n,ih , ξ ) = 0 ∀ ξ ∈ Σ h , (5.1)(( c + β ) p n,ih ∆ t , w ) + α ( A ( αp n,ih I )∆ t , wI ) − b h ( z n,ih , w ) = ( q n , w ) + ( c p n − h ∆ t , w )+ ( β p n,i − h ∆ t , w ) − α ( A ( σ n,i − h ) − A ( σ n − h + αp n − h I )∆ t , wI ) ∀ w ∈ S h . (5.2) Step 2 : Given p n,ih , find ( σ n,ih , u n,ih , γ n,ih ) such that( A ( σ n,ih ) , ψ ) − B ∗ h ( u n,ih , ψ ) + ( γ n,ih , as ( ψ )) = − ( A ( αp n,ih I ) , ψ ) ∀ ψ ∈ [Σ h ] , (5.3) B h ( σ n,ih , v ) = ( f n , v ) ∀ v ∈ [ S h ] , (5.4)( as ( σ n,ih ) , η ) = 0 ∀ η ∈ M h . (5.5)The initial guess for the iterations is chosen to be the solution at the last time step, i.e., ( σ n, h , p n, h ) =( σ n − h , p n − h ). Remark 5.1.
Since the mass matrix is block diagonal, we can apply local elimination for (5.1)-(5.2)and (5.3)-(5.5), respectively. In this way, we can get a reduced system which solely depends on thedisplacement and pressure.
Theorem 5.1. (linear convergence of fixed stress splitting). Let ( σ nh , u nh , γ nh , z nh , p nh ) be the solution of (4.24) - (4.28) and let ( σ n,ih , u n,ih , γ n,ih , z n,ih , p n,ih ) be the solution of (5.1) - (5.5) . Then for all β satisfying β ≥ α µ + λ ) . ocking free SDG for poroelasticity We have k p nh − p n,ih k ≤ β c + β + α µ + λ + ∆ tC − p C − inf K k p nh − p n,i − h k , ∀ i ≥ , (5.6) where C p is the Poincar´e constant and C inf is the inf-sup constant in (3.6) .Proof. Let e iu = − u n,ih + u nh , e iσ = − σ n,ih + σ nh , e iγ = − γ n,ih + γ nh , e ip = − p n,ih + p nh and e iz = − z n,ih + z nh .Subtracting (5.1)-(5.5) from (4.24)-(4.28), we obtain( A ( e iσ + αe ip I ) , ψ ) − B ∗ h ( e iu , ψ ) + ( e iγ , as ( ψ )) = 0 ,B h ( e iσ , v ) = 0 ,as ( e iσ , η ) = 0 , ( K − e iz , ξ ) + b ∗ h ( e ip , ξ ) = 0 ,c ∆ t ( e ip , w ) − β ∆ t ( p n,ih − p n,i − h , w ) + α ∆ t ( A ( e i − σ + αe ip I ) , wI ) − b h ( e iz , w ) = 0for ( ψ, v, η, ξ, w ) ∈ [Σ h ] × [ S h ] × M h × Σ h × S h .Taking ψ = e i − σ , v = e iu , η = e iγ , ξ = e iz , w = e ip in the above equations, and summing up theresulting equations, we can infer that c k e ip k + β ( e ip − e i − p , e ip ) + α ( A ( e ip I ) , e ip I ) + ( A ( e i − σ ) , e iσ ) + ∆ t k K − e iz k = 0 . (5.7)Consider (5.3)-(5.4) at iterations i and i − ψ = σ n,ih − σ n,i − h andsumming up the resulting equations yield( A ( σ n,ih − σ n,i − h ) , σ n,ih − σ n,i − h ) = − α ( A (( p n,ih − p n,i − h ) I ) , σ n,ih − σ n,i − h ) , which leads to k A ( σ n,ih − σ n,i − h ) k ≤ α k A (( p n,ih − p n,i − h ) I ) k k A ( σ n,ih − σ n,i − h ) k ≤ (cid:16) α k A (( p n,ih − p n,i − h ) I ) k + 12 k A ( σ n,ih − σ n,i − h ) k (cid:17) . Thus, we have k A / ( e iσ − e i − σ ) k ≤ α k A / (( p n,ih − p n,i − h ) I ) k ≤ α k A / (( e ip − e i − p ) I ) k ≤ α ( λ + µ ) / k e ip − e i − p k . (5.8)We can infer from (5.7), the equality ( a − b, a ) = a − b +( a − b ) and the identity ( A ( e i − σ ) , e iσ ) = (cid:16) k A ( e iσ + e i − σ ) k − k A ( e iσ − e i − σ ) k (cid:17) that c k e ip k + β k e ip k − k e i − p k + k e ip − e i − p k ) + α ( A ( e ip I ) , e ip I )+ 14 (cid:16) k A ( e iσ + e i − σ ) k − k A ( e iσ − e i − σ ) k (cid:17) + ∆ t k K − e iz k = 0 . Therefore, an appeal to (5.8) yields c k e ip k + β k e ip k − k e i − p k ) + ( β − α λ + µ ) ) k e ip − e i − p k + α ( A ( e ip I ) , e ip I )+ 14 k A ( e iσ + e i − σ ) k + ∆ t k K − e iz k ≤ . It follows from the inf-sup condition (3.6) and the discrete Poincar´e-Friedrichs inequality that k e ip k ≤ C p C inf k K − e iz k , ocking free SDG for poroelasticity C p is the Poincar´e constant and C inf is the inf-sup constant in (3.6). Thereby, if we require β ≥ α λ + µ ) , then there holds c k e ip k + β k e ip k − k e i − p k ) + α ( A ( e ip I ) , e ip I )+ ∆ tC − p C − inf K k e ip k + 14 k A ( e iσ + e i − σ ) k ≤ , (5.9)which leads to the desired estimate (5.6). In this section we present several numerical experiments to verify the theoretical convergence ratesand illustrate the behavior of the proposed method. We also present one benchmark example showingthe locking-free property of the method. In the simulation presented below, we use the polynomialorder k = 1. In this example, we test the convergence of our method. To this end, we set Ω = (0 , , where thedisplacement and pressure are given by u = (cid:18) sin(2 πt ) sin(2 πx ) sin(2 πy ) x cos( t ) (cid:19) , p = exp( t )(sin( πx ) cos( πy ) + 10) . The corresponding f and q can be calculated by (2.1)-(2.3). We set α = 1 , λ = 1 , µ = 1 and thesimulation time is T = 0 .
01 with time steps ∆ t = h . We remark that ∆ t should be chosen smallenough so that the time discretization error will not influence the convergence rates. For this example,we consider two values of K , i.e., K = I and K = (cid:18)
00 1 (cid:19) (6.1)We exploit square meshes and trapezoidal meshes for this example, see Figure 2 for an illustration.The convergence history against the number of degrees of freedom for rectangular meshes and trape-zoidal meshes with various values of c is depicted in Figure 3. Optimal convergence rates matchingthe theoretical results can be obtained for various values of c . In addition, with different values of c ,the accuracy for all the errors remain almost the same and the method is still valid for c = 0. Fur-ther, the accuracy for L errors of σ, u, γ, z, p on rectangular meshes and trapezoidal meshes is slightlydifferent. Figure 4 shows the convergence history on rectangular meshes for anisotropic K defined in(6.1) and similar performances can be observed, which illustrates the robustness of our method to theanisotropy ratio. The computational domain is the unit square. We impose no-flow boundary condition along all sides.The deformation is fixed at the left edge, and a downward traction is applied along the top. Thebottom and right sides are traction-free, see Figure 5 for an illustration of the boundary condition.The time step for this example is ∆ t = 0 . E = 10 , ν = 0 . , α = 0 . , c = 0 , K = 10 − , ocking free SDG for poroelasticity Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Figure 3: Example 6.1. Top plots: convergence history on square meshes for c = 1 (left), c = 10 − (middle) and c = 0 (right) for K = 1. Bottom plots: convergence history on trapezoidal meshes for c = 1 (left), c = 10 − (middle) and c = 0 (right) for K = 1.where E is the Young’s modulus and ν is the Poisson ratio, and the Lam´e parameters are computedby λ = Eν (1 + ν )(1 − ν ) , µ = E ν ) . Figure 6 shows that our proposed method yields a smooth pressure field, in contrast to the non-physicalcheckboard pattern that one obtains with continuous elasticity elements at the early time steps, see[39]. In addition, Figure 6 also indicates that the pressure solution along different x -lines at time T = 0 . The computational domain is taken to be Ω = (0 , and we set α = 1 , λ = 1 , µ = 1 , K = I . Thesimulation time is T = 0 .
01 with time steps ∆ t = h . The exact displacement and pressure are givenas u = (cid:18) sin(2 πt ) sin(2 πy )exp( − x/ϑ )sin(2 πt ) cos(2 πy )exp( − x/ϑ ) (cid:19) , p = exp( t ) exp( − x/ϑ ) x (1 − x ) y (1 − y ) . ocking free SDG for poroelasticity Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Degrees of freedom -5 -4 -3 -2 -1 E rr o r Convergence history h || || - h || || - h || ||p -p h || ||z-z h || Figure 4: Example 6.1. Convergence history on square meshes for c = 1 (left), c = 10 − (middle)and c = 0 (right) for anisotropic K .Figure 5: Example 6.2. Profile of boundary condition.The used meshes are of Shishkin type (cf. [2]), see the example in Figure 8. For a parameter N ≥ δ ∈ (0 ,
1) and generating a grid of points( x j , y j ) by x j = ( j δN , ≤ j ≤ N , j ∈ N ,δ + ( j − N ) − δ ) N , N < j ≤ N, j ∈ N ,y i = iN ≤ i ≤ N, i ∈ N . The transition point parameter δ is chosen as δ = min { , ϑ | ln( ϑ ) |} . Here we choose ϑ = 10 − andthe corresponding exact velocity and pressure at t = 0 .
01 are displayed in Figures 7 and 8, wherethe exponential boundary layer near x = 0 is clearly visible. The layer has a width of O ( ϑ ) and ispresent in the velocity and pressure solution. The convergence history against the number of degreesof freedom for various values of c is shown in Figure 9 on square meshes and anisotropic meshes. Wecan observe that optimal convergence rates matching the theoretical results can be obtained, whichindicates that our method can work well on anisotropic meshes and the method works for c = 0.Further, the errors on anisotropic meshes are much smaller than that of rectangular meshes, whichillustrates the superior performances of anisotropic meshes in particular for layered problem. In this example, we consider a more practical case motivated by the example given in [25], where threesedimentary layers overlying an impermeable bedrock in a basin. The sediment stack totals 420m atthe deepest point of the basin ( x = 0m) but thins to 120m above the step ( x > ocking free SDG for poroelasticity y axis -0.6-0.4-0.200.20.40.60.811.2 P r e ss u r e x=0.26x=0.33x=0.4x=0.46 Figure 6: Example 6.2. Pressure field at T = 0 .
001 (left) and pressure at different x -lines, T = 0 . u (left) and u (right).region A to E are described in Table 1, in addition, we let f = q = 0. The values for the parametersare given in Table 2. The final simulation time is 10 years, where we take the simulation time to be(0 : 360 : 360 × In this paper we analyzed a staggered DG method for a five-field formulation of the Biot system ofporoelasticity on general polygonal meshes. The method is proved to converge in optimal rates for allthe variables in their natural norms. In addition, the error estimates are independent of the storativitycoefficient c and are valid even for c = 0, which is also confirmed by our numerical simulation.Several numerical experiments including cantilever bracket benchmark problem are presented, whichfurther confirms the locking-free property of our method. In addition, the accuracy of our method isslightly influenced by the shape of the grid. Numerical results illustrate that our method is a goodoption for solving problems on anisotropic meshes. Moreover, our method can handle nonmatchingmeshes straightforwardly, which is advantageous in solving problem with local features. In the futurewe will extend this method to solve nonlinear poroelasticity. Acknowledgment
The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund(Project numbers 14304719 and 14302018) and CUHK Faculty of Science Direct Grant 2019-20 andNSFC/RGC Joint Research Scheme (Project number HKUST620/15). The research of Eun-Jae Parkwas supported by the National Research Foundation of Korea (NRF) grant funded by the Ministry ofScience and ICT (NRF-2015R1A5A1009350 and NRF-2019R1A2C2090021). ocking free SDG for poroelasticity ϑ = 10 − , N = 8 used inthe simulations. Degrees of freedom -3 -2 -1 E rr o r Convergence history h || ,Rect|| - h || ,Rect|| - h || ,Rect||p -p h || ,Rect||z-z h || ,Rect||u-u h || ,Aniso|| - h || ,Aniso|| - h || ,Aniso||p -p h || ,Aniso||z-z h || ,Aniso 10 Degrees of freedom -3 -2 -1 E rr o r Convergence history h || ,Rect|| - h || ,Rect|| - h || ,Rect||p -p h || ,Rect||z-z h || ,Rect||u-u h || ,Aniso|| - h || ,Aniso|| - h || ,Aniso||p -p h || ,Aniso||z-z h || ,Aniso 10 Degrees of freedom -3 -2 -1 E rr o r Convergence history h || ,Rect|| - h || ,Rect|| - h || ,Rect||p -p h || ,Rect||z-z h || ,Rect||u-u h || ,Aniso|| - h || ,Aniso|| - h || ,Aniso||p -p h || ,Aniso||z-z h || ,Aniso Figure 9: Example 6.3. Convergence history with uniform rectangular meshes (solid line) andanisotropic meshes (dashed line). c = 1 (left), c = 10 − (middle), and c = 0 (right). References [1] I. Ambartsumyan, E. Khattatov, and I. Yotov. A coupled multipoint stress – multipoint fluxmixed finite element method for the Biot system of poroelasticity. arXiv:2001.04582 , 2020.[2] T. Apel, A. Linke, and C. Merdon. A nonconforming pressure-robust finite element method forthe Stokes equations on anisotropic meshes. arXiv:2002.12127 , 2020.[3] M. A. Biot. General theory of three-dimensional consolidation.
J. Appl. Phys. , 12(2):155–164,1941.[4] D. Boffi, M. Botti, and D. A. Di Pietro. A nonconforming high-order method for the Biot problemon general meshes.
SIAM J. Sci. Comput. , 38(3):A1508–A1537, 2016.[5] D. Braess.
Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics . Cam-bridge University Press, 2007.[6] S. C. Brenner. Poincar´e-Friedrichs inequalities for piecewise H functions. SIAM J. Numer. Anal. ,41(1):306–324, 2003.[7] F. Brezzi and R. S. Falk. Stability of higher-order Hood-Taylor methods.
SIAM J. Numer. Anal. ,28(3):581–590, 1991.[8] S. Cheung, E. T. Chung, H. H. Kim, and Y. Qian. Staggered discontinuous Galerkin methods forthe incompressible Navier–Stokes equations.
J. Comput. Phys. , 302(1):251–266, 2015.[9] E. T. Chung, B. Cockburn, and G. Fu. The staggered DG method is the limit of a hybridizableDG method.
SIAM J. Numer. Anal. , 52(2):915–932, 2014.[10] E. T. Chung, B. Cockburn, and G. Fu. The staggered DG method is the limit of a hybridizableDG method. Part II: The Stokes flow.
J. Sci. Comput. , 66(2):870–887, 2016. ocking free SDG for poroelasticity × z · n = 0 u = (0 , T B z · n = 0 u = 0C p = 0 u = 0D p = 0 σn = (0 , T E p = H ( t ) u = 0Table 1: Example 6.4. Description of boundary conditions. u = ( u , u ) T .[11] E. T. Chung, J. Du, and C. Lam. Discontinuous Galerkin methods with staggered hybridizationfor linear elastodynamics. Comput. Math. Appl. , 74(6):1198–1214, 2017.[12] E. T. Chung and B. Engquist. Optimal discontinuous Galerkin methods for wave propagation.
SIAM J. Numer. Anal. , 44(5):2131–2158, 2006.[13] E. T. Chung and B. Engquist. Optimal discontinuous Galerkin methods for the acoustic waveequation in higher dimensions.
SIAM J. Numer. Anal. , 47(5):3820–3848, 2009.[14] E. T. Chung, H. H. Kim, and O. B. Widlund. Two-level overlapping Schwarz algorithms for astaggered discontinuous Galerkin method.
SIAM J. Numer. Anal. , 51(1):47–67, 2013.[15] E. T. Chung and C. Lee. A staggered discontinuous Galerkin method for the curl–curl operator.
IMA J. Numer. Anal. , 32(3):1241–1265, 2012.[16] E. T. Chung and W. Qiu. Analysis of an SDG method for the incompressible Navier–Stokesequations.
SIAM J. Numer. Anal. , 55(2):543–569, 2017.[17] P. G. Ciarlet.
The Finite Element Method for Elliptic Problems . Society for Industrial and AppliedMathematics, 2002.[18] J. Du and E. T. Chung. An adaptive staggered discontinuous Galerkin method for the steadystate convection–diffusion equation.
J. Sci. Comput. , 77(3):1490–1518, 2018.[19] G. Fu. A high-order HDG method for the Biot’s consolidation model.
Comput. Math. Appl. ,77(1):237–252, 2019.[20] F. J. Gaspar, F. J. Lisbona, and P. N. Vabishchevich. A finite difference analysis of Biot’sconsolidation model.
Appl. Numer. Math. , 44(4):487–506, 2003.[21] J. Guzm´an and L. R. Scott. The Scott-Vogelius finite elements revisited.
Math. Comp. ,88(316):515–529, 2019.[22] X. Hu, L. Mu, and X. Ye. Weak Galerkin method for the Biot’s consolidation model.
Comput.Math. Appl. , 75(6):2017–2030, 2018.[23] H. H. Kim, E. T. Chung, and C. Lee. A staggered discontinuous Galerkin method for the Stokessystem.
SIAM J. Numer. Anal. , 51(6):3327–3350, 2013. ocking free SDG for poroelasticity c storage coefficient, aquifer layers 10 − m − c storage coefficient, confining layers 10 − m − K permeability, aquifer layers 25 m/day K permeability, confining layers 0.01 m/day α Biot-Willis coefficient 1E Young’s modulus, aquifer layers 8 · PaE Young’s modulus, confining layers 8 · Pa ν Poisson ratio, all regions 0.25 H ( t ) Declining head boundary (6 m/year) · t Table 2: Example 6.4. Values of parameters.Figure 11: Example 6.4. Numerical results at 2 years (left) and 10 years (right). Surface: pressure;Contour: displacement ( | u h | ).[24] J. Korsawe and G. Starke. A least-squares mixed finite element method for Biot’s consolidationproblem in porous media. SIAM J. Numer. Anal. , 43(1):318–339, 2005.[25] S. A. Leake and P. A. Hsieh. Simulation of deformation of sediments from decline of ground-waterlevels in an aquifer underlain by a bedrock step.
U.S. Geological Survey Open File Report: 97-47 ,1997.[26] J. J. Lee. Robust error analysis of coupled mixed methods for Biot’s consolidation model.
J. Sci.Comput. , 69(2):610–632, 2016.[27] J. J. Lee. Robust three-field finite element methods for Biot’s consolidation model in poroelasticity.
BIT Numer. Math , 58(2):347–372, 2018.[28] J. J. Lee and H. H. Kim. Analysis of a staggered discontinuous Galerkin method for linearelasticity.
J. Sci. Comput. , 66(2):625–649, 2016.[29] J. J. Lee, K.-A. Mardal, and R. Winther. Parameter-robust discretization and preconditioning ofBiot’s consolidation model.
SIAM J. Sci. Comput. , 39(1):A1–A24, 2017.[30] A. Mikeli´c and M. F. Wheeler. Convergence of iterative coupling for coupled flow and geome-chanics.
Comput. Geosci. , 17(3):455–461, 2013.[31] M. A. Murad and A. F. D. Loula. Improved accuracy in finite element analysis of Biot’s consoli-dation problem.
Comput. Methods Appl. Mech. Engrg. , 95(3):359–382, 1992.[32] M. A. Murad and A. F. D. Loula. On stability and convergence of finite element approximationsof Biot’s consolidation problem.
Internat. J. Numer. Methods Engrg. , 37(4):645–667, 1994.[33] M. A. Murad, V. Thom´ee, and A. F. D. Loula. Asymptotic behavior of semidiscrete finite-elementapproximations of Biot’s consolidation problem.
SIAM J. Numer. Anal. , 33(3):1065–1083, 1996. ocking free SDG for poroelasticity
Comput. Geosci. , 2020.[35] J. M. Nordbotten. Stable cell-centered finite volume discretization for Biot equations.
SIAM J.Numer. Anal. , 54(2):942–968, 2016.[36] R. Oyarz´ua and R. Ruiz-Baier. Locking-free finite element methods for poroelasticity.
SIAM J.Numer. Anal. , 54(5):2951–2973, 2016.[37] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin finite elementmethods for poroelasticity I: the continuous in time case.
Comput. Geosci. , 11(2):131–144, 2007.[38] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin finite elementmethods for poroelasticity II: the discrete-in-time case.
Comput. Geosci. , 11(2):145–158, 2007.[39] P. J. Phillips and M. F. Wheeler. Overcoming the problem of locking in linear elasticity andporoelasticity: an heuristic approach.
Comput. Geosci. , 13(1):5–12, 2009.[40] C. Rodrigo, F. J. Gaspar, X. Hu, and L.T. Zikatanov. Stability and monotonicity for some dis-cretizations of the Biot’s consolidation model.
Comput. Methods Appl. Mech. Engrg. , 298(1):183–204, 2016.[41] R. E. Showalter. Diffusion in poro-elastic media.
J. Math. Anal. Appl. , 251(1):310–340, 2000.[42] R. E. Showalter. Nonlinear degenerate evolution equations in mixed formulation.
SIAM J. Numer.Anal. , 42(5):2114–2131, 2010.[43] R. E. Showalter.
Monotone operators in Banach space and nonlinear partial differential equations .American Mathematical Society, 2013.[44] M. Wheeler, G. Xue, and I. Yotov. Coupling multipoint flux mixed finite element methods withcontinuous Galerkin methods for poroelasticity.
Comput. Geosci. , 18(1):57–75, 2014.[45] S.-Y. Yi. A coupling of nonconforming and mixed finite element methods for Biot’s consolidationmodel.
Numer. Methods Partial Differential Equations , 29(5):1749–1777, 2013.[46] S.-Y. Yi. Convergence analysis of a new mixed finite element method for Biot’s consolidationmodel.
Numer. Methods Partial Differential Equations , 30(4):1189–1210, 2014.[47] S.-Y. Yi. A study of two modes of locking in poroelasticity.
SIAM J. Numer. Anal. , 55(4):1915–1936, 2017.[48] L. Zhao, E. T. Chung, and M. Lam. A new staggered DG method for the Brinkman problemrobust in the Darcy and Stokes limits.
Comput. Methods Appl. Mech. Engrg. , 364(1), 2020.[49] L. Zhao, E. T. Chung, E.-J. Park, and G. Zhou. Staggered DG method for coupling of the Stokesand Darcy-Forchheimer problems.
SIAM J. Numer. Anal. , 2020, to appear.[50] L. Zhao, D. Kim, E.-J. Park, and E. Chung. Staggered DG method with small edges for Darcyflows in fractured porous media. arXiv:2005.10955 , 2020.[51] L. Zhao and E.-J. Park. A staggered discontinuous Galerkin method of minimal dimension onquadrilateral and polygonal meshes.
SIAM J. Sci. Comput. , 40(4):A2543–A2567, 2018.[52] L. Zhao and E.-J. Park. A lowest-order staggered DG method for the coupled Stokes–Darcyproblem.
IMA J. Numer. Anal. , 40(4):2871–2897, 2020.[53] L. Zhao and E.-J. Park. A staggered cell-centered DG method for linear elasticity on polygonalmeshes.
SIAM J. Sci. Comput. , 42(4):A2158–A2181, 2020.[54] L. Zhao, E.-J. Park, and D. w. Shin. A staggered DG method of minimal dimension for the Stokesequations on general meshes.