1 Introduction
Modeling flow in fractured porous media has drawn great attention in the past decades, being fundamental for addressing many environmental and energy problems, such as water resources management, isolation of radioactive waste and ground water contamination. Given the wide applications of fractured model in practical applications, many advances has been made in the accomplishments of designing efficient numerical methods for fractured porous media. In [39], a mixed finite element method is developed and error estimates are also proved. Later, a mixed finite element method on nonmatching grids is considered in [31]. In [17], a hybridhigh order method is analyzed on fairly general meshes. The error estimates proposed therein show that the method is fully robust with respect to the heterogeneity of the permeability coefficients. In [2], a discontinuous Galerkin approximation for flows in fractured porous media on polytopal grids is analyzed, where optimal convergence estimates in meshdependent energy norm are derived on fairly general meshes possibly including elements with unbounded number of faces. In addition to the aforementioned methods, we also mention other methods that have been developed for fractured porous media, see [38, 1, 41, 8, 35, 9, 12, 3, 19, 20, 13, 32, 11, 34].
Staggered discontinuous Galerkin (DG) methods are initially introduced to solve wave propagation problems [24, 25]
. The salient features of staggered DG method make it desirable for practical applications and the applications to various partial differential equations important for both science and engineering have been considered in
[27, 26, 36, 22, 37, 29, 28, 43]. Recently, staggered DG methods have been successfully designed on fairly general meshes possibly including hanging nodes for Darcy law and the Stokes equations, respectively [44, 45]. It is further developed with essential modifications to solve coupled Stokes and Darcy problem, and Brinkman problem [46, 47]. Staggered DG methods designed therein earn many desirable features, including: 1) It can be flexibly applied to fairly general meshes with possible inclusion of hanging nodes, and the handing nodes can be simply incorporated in the construction of the method; 2) superconvergence can be obtained, which can deliver one order high convergence with proper postprocessing scheme designed; 3) local mass conservations can be preserved, which is highly appreciated in the simulation of multiphase flow. In addition, the mass matrix is block diagonal which is desirable when explicit time stepping schemes are used; 4) no numerical flux or penalty term is needed in contrast to other DG methods.The purpose of this paper is to develop and analyze staggered DG method for the coupled bulkfracture model stemming from the modeling of flows in fractured porous media, allowing more general meshes such as elements with arbitrarily small edges. The flexibility of staggered DG method in handling fairly general meshes, and the preservation of physical properties indeed make it an attractive candidate for such kind of problems. In this paper we propose a discretization which combines a staggered DG approximation for the problem in the bulk domains with a conforming finite element approximation on the fracture. Unlike the strategies employed in [31, 17]
, we impose the coupling conditions by replacing all the terms with respect to the jump and average of flux by the corresponding pressure term, which can compensate for the degrees of freedom for bulk pressure across the fracture. The existence and uniqueness of the resulting system is proved and a rigorous error analysis is carried out. In particular, we prove the convergence estimates under weaker assumption on the polygonal mesh by exploiting some novel strategies. Research in this direction has drawn great attention, see
[7, 10, 15, 2, 16] for works considering general polygonal elements allowing arbitrarily small edges. The primary difficulty arising from a priori error estimates lies in the fact that error estimate for flux is coupled with energy error of fracture pressure, which will naturally lead to suboptimal convergence for error of flux. To overcome this issue, we construct the Ritz projection for fracture pressure so that the term causing suboptimal convergence can vanish. Moreover, we are able to show that the Ritz projection superconverges to numerical approximation of fracture pressure. Then without duality argument we can achieve optimal convergence for error of fracture pressure and bulk pressure, respectively. It is noteworthy that our error estimates are shown to be fully robust with respect to the heterogeneity and anisotropy of the permeability coefficients, which is desirable feature for fractured flow simulation. The theoretical findings are verified by a series of numerical tests. Especially, numerical tests indicate that our method is robust to anisotropy of meshes. We emphasize that our method allows general meshes with arbitrarily small edges, thus it can be easily adapted to solve problem on unfitted grids. In fact, we only need to update the interface elements by connecting the intersection points between background grids and fracture, thereby the resulting grids are again fitted with fracture and thus can be naturally embedded into our current framework. Therefore, this paper focuses on the heart of the novelty on the fitted mesh to make the presentation clear.The rest of this paper is organized as follows. In the next section, we describe the model problem and formulate the staggered DG formulation for the bulk region coupled with standard conforming Galerkin formulation inside the fracture. In addition, some fundamental ingredients are given in order to prove the a priori error estimates. In Section 3, a priori error analysis is derived for bulk flux, bulk pressure and fracture pressure measured in error, where a discrete trace inequality is proved. Then several numerical experiments are given in Section 4 to confirm the theoretical findings, where various tests including elements with small edges and anisotropic meshes are demonstrated. Finally, a conclusion is given.
2 Description of staggered DG method
In this section we first describe the governing equations modeling Darcy flows in fractured porous media. Then staggered DG discretization is derived for the model problem under consideration. Finally, we introduce some technical results that are vital for subsequent sections.
2.1 Model problem
We consider a porous medium saturated by an incompressible fluid that occupies the space region and is crossed by a single fracture . Here, represents the bulk region and can be decomposed as . In addition, we denote by and denote by the boundary of fracture .
denotes a unit normal vector to
with a fixed orientation. The schematic of the bulk and fracture domain is illustrated in Figure 1. Without loss of generality, we assume in the following that the subdomains are numbered so that coincides with the outward normal direction of .In the bulk region, we model the motion of the incompressible fluid by Darcy’s law in mixed form, so that the pressure and the flux satisfy
(2.1)  
(2.2)  
(2.3) 
Here, the boundary pressure, and
the bulk permeability tensor, which is assumed to be a symmetric, piecewise constant. Further, we assume that
is uniformly elliptic so that there exist two strictly positive real numbers and satisfying for almost every and all such thatInside the fracture, we consider the motion of the fluid as governed by Darcy’s law in primal form, so that the fracture pressure satisfies
(2.4)  
where and with and denoting the tangential permeability and thickness of the fracture, respectively. The quantities and are assumed to be piecewise constants. Here, and denote the tangential divergence and gradient operators along , respectively. For the sake of simplicity, we assume , in the analysis.
The above problems are coupled by the following interface conditions
(2.5)  
where we set
Here is a model parameter, and represents the normal permeability of the fracture, which is assumed to be a piecewise constant. As in the bulk domain, we assume that there exists positive constants such that, almost everywhere on ,
Also, and are jump and average operators, respectively, and their precise definitions can be found in the next subsection. The wellposedness of the coupled problem for has been proved in [39].
Remark 2.1 (Neumann boundary conditions).
When the fracture tip is immersed in the domain , the boundary condition at the immersed tip can be modeled as a homogeneous Neumann boundary condition, see [1]. For both bulk and fracture domains, the Neumann boundary condition can be treated as a natural boundary condition. Since the analysis for such boundary condition is parallel to the analysis for Dirichlet boundary conditions, we only consider the latter one for simplicity.
Before closing this subsection, we introduce some notations that will be employed throughout the paper. Let , we adopt the standard notations for the Sobolev spaces and their associated norms , and seminorms for . The space coincides with , for which the norm is denoted as . We use to denote the inner product for and for . If , the subscript will be dropped unless otherwise mentioned. In the sequel, we use to denote a generic positive constant which may have different values at different occurrences.
2.2 Staggered DG method
In this subsection, we begin with introducing the construction of our staggered DG spaces, in line with this we then present the staggered DG method for the model problem (2.1)(2.5). We consider a family of meshes made of disjoint polygonal (primal) elements which are aligned with the fracture so that any element can not be cut by . Note that, since and are disjoint, each element belongs to one of the two subdomains. The union of all the edges excluding the edges lying on the fracture in the decomposition is called primal edges, which is denoted as . Here we use to stand for the subset of , that is the set of edges in that do not lie on . In addition, we use to denote the onedimensional mesh of the fracture . For the construction of staggered DG method, we decompose each element into the union of triangles by connecting the interior point of to all the vertices. Here the interior point is chosen as the center point for simplicity. We rename the union of these subtriangles by to indicate that the triangles sharing common vertex . In addition, the resulting simplicial submeshes are denoted as . Moreover, some additional edges are generated in the subdivision process due to the connection of to all the vertices of the primal element, and these edges are denoted by . For each triangle , we let be the diameter of and . In addition, we define and . The construction for general meshes is illustrated in Figure 2, where the black solid lines are edges in and black dotted lines are edges in .
Finally, we construct the dual mesh. For each interior edge , we use to represent the dual mesh, which is the union of the two triangles in sharing the edge . For each edge , we use to denote the triangle in having the edge , see Figure 2.
For each edge , we define a unit normal vector as follows: If , then is the unit normal vector of pointing towards the outside of . If , an interior edge, we then fix as one of the two possible unit normal vectors on . When there is no ambiguity, we use instead of to simplify the notation.
Typical analysis for polygonal element usually requires the following mesh regularity assumptions (cf. [6, 14]):
 Assumption (A)

Every element in is starshaped with respect to a ball of radius , where is a positive constant and denotes the diameter of .
 Assumption (B)

For every element and every edge , it satisfies , where is a positive constant and denotes the length of edge .
Assumption (A) and (B) can guarantee that the triangulation is shape regular. However, it excludes the elements with arbitrarily small edges, which is interesting from the practical applications. Thus, in this paper, we will show the convergence estimates by only assuming Assumption (A).
Let be the order of approximation. For every and , we define and as the spaces of polynomials of degree less than or equal to on and , respectively. For and belonging to the broken Sobolev space the jump and the jump over are defined respectively as
where , and , are the two triangles in having the edge . Moreover, for , we define . In the above definitions, we assume is pointing from to .
Similarly, we define the average and the average over by
where , and , are the two triangles in having the edge .
Next, we will introduce some finite dimensional spaces. First, we define the following locally conforming space :
Notice that, if , then for each edge and no continuity is imposed across for function . The discrete norm for are defined as follows
where represents the area of triangle . Note that the scaling in the second term used here is different from that of [44], and this modification enables us to show the convergence estimates without Assumption (B). We specify the degrees of freedom for similar to that of [25].
(SD1) For , we have
(SD2) For , we have
(SD3) For , we have for each
Note that in original staggered DG method, the finite dimensional space for pressure is continuous over all the primal edges, in which case (SD3) can be compliant with (SD1). In this paper we consider Darcy flows with fracture where the pressure is discontinuous across the fracture, thereby (SD3) can not be compliant with (SD1). Proceeding analogously to Lemma 2.2 of [25], we can show that any function is uniquely determined by the degrees of freedom (SD1)(SD3), which is omitted here for simplicity.
We next define the following locally conforming space :
Note that if , then for each . We equip with the following discrete norm
The degrees of freedom for can be defined below.
(VD1) For each edge , we have
(VD2) For each , we have
Finally, we define a finite dimensional subspace of by
With the above preparations, we can now derive our staggered DG method by following [44, 45]. Multiplying (2.1) by and performing integration by parts, we can obtain
where the staggered continuity property of is integrated into the derivation.
Similarly, multiplying (2.2) by and performing integration by parts yield
(2.6) 
Then we exploit the interface condition (2.5) in (2.6) and recast the above formulation as
As for the fracture model (2.4), we multiply by and replace the jump term by utilizing (2.5), which implies
Thereby we obtain the following discrete formulation for the model problem (2.1)(2.4): Find such that
(2.7) 
where the bilinear forms are defined by
Summing up the equations in (2.7) yields the following formulation: Find such that
(2.8) 
Integration by parts reveals the following adjoint property
(2.9) 
Remark 2.2.
In the derivation we employ the interface conditions (2.5) to replace all the terms corresponding to on the fracture by and , which is different from existing methods such as the hybrid highorder method and mixed finite element method [17, 31]. Our methodology is based on the fact that the degrees of freedom for bulk pressure (SD3) are defined with respect to the primal edges on the fracture. We also emphasize that the velocity can be made both locally and globally mass conservative by a suitable postprocessing (cf. [23]). Moreover, the use of conforming finite element to discretize the equations in the fracture is made just for simplicity, other discretization techniques can be exploited.
Lemma 2.1.
Under Assumption (A), we have the following infsup condition
(2.10) 
Proof.
The proof for this lemma follows similar idea as Theorem 3.2 of [25], we simplify the proof by direct applications of the degrees of freedom (VD1)(VD2). In addition, our proof here only relies on Assumption (A) thanks to the modified norm defined for and .
Let . It suffices to find such that
Recall that
(2.11) 
We define by using degrees of freedom (VD1)
and (VD2)
which together with (2.11) yields
On the other hand, scaling arguments imply
This completes the proof.
∎
Finally, we introduce the following interpolation operators, which play an important role in later analysis. We define the interpolation operator
byand the interpolation operator by
The definition of the interpolation operators implies that
Notice that if is continuous on , then is also continuous on . The interpolation operators and satisfy: (1) they are locally defined for each element ; (2) for and , we have and .
The following error estimates are clearly satisfied on the reference element by using (1) and (2) (see [30]).
where and are the corresponding variables of and on the reference element . In addition, under Assumption (A), the maximum angles in are uniformly bounded away from (although shape regularity is not guaranteed). Then we can proceed as Theorem 2.1 of [4] to obtain the following anisotropic error estimates for
(2.12) 
Here, generic constants possibly depend on in Assumption (A) but not on in Assumption (B). We next introduce the standard nodal interpolation operator , which satisfies for
where .
3 Error analysis
In this section, we present the unique solvability of the discrete system (2.7) and the convergence estimates for all the variables involved under Assumption (A). As error of is coupled with energy error of , it will yield suboptimal convergence if standard interpolation operator for is exploited. As such, we propose to employ the Ritz projection which enables us to achieve the optimal convergence estimates.
Theorem 3.1 (stability).
Under Assumption (A), the discrete system (2.7) admits a unique solution . Furthermore, there exists a positive constant independent of but possibly depending on and the problem data such that
(3.1) 
Proof.
Since (2.8) is a square linear system, existence follows from uniqueness, thus, it suffices to show uniqueness. Taking in (2.8) yields
On the other hand, an application of the discrete PoincaréFriedrichs inequality on anisotropic meshes (cf. [33]) leads to
In view of the infsup condition (2.10), the discrete adjoint property (2.9), and (2.7), we have
Moreover, and the Poincaré inequality imply that
Combining the above estimates with Young’s inequality, we can infer that
which gives the desired estimate (3.1). Here, depends on the permeability and . The uniqueness follows immediately by setting . ∎
Here, we introduce the Ritz projection , which is defined by
(3.2) 
It is wellposed by the Riesz representation theorem. Then taking in (3.2) yields
which implies
Next, we show the error estimate for
Comments
There are no comments yet.