1 Introduction
The fluid flow between porous media and freeflow zones have extensive applications in hydrology, environment science and biofluid dynamics, which motivates a lot of researchers to derive suitable mathematical and numerical models for the fluid movement. The system can be viewed as a coupled problem with two physical systems interacting across an interface. The simplest mathematical formulation for the coupled problem is the coupled Stokes and Darcy flows problems, which imposes the Stokes equations in the Stokes region and the Darcy’s law in the porous media region, coupled with proper interface conditions. The most suitable and popular interface conditions are called BeaversJosephSaffman conditions, which are proposed by Saffman [37]. The coupled model has received great attention and a large number of works have been devoted to the coupled Stokes and Darcy flows problem, see, for example [18, 29, 33, 34, 9, 21, 23, 30]. However, Darcy’s law only provides linear relationship between the gradient of pressure and velocity, which usually fails for complex physical problems. Forchheimer [20] conducted flow experiments in sandpacks and recognized that for moderate Reynolds numbers (Re approximately), Darcy’s law is not adequate. He found that the pressure gradient and Darcy velocity should satisfy nonlinear relationship and the DarcyForchheimer law accounts for this nonlinear behavior. Mixed methods for the steady DarcyForchheimer flow was studied in more general setting in [31] and unsteady case was analyzed in [27, 32]. A nonconforming primal mixed method for the DarcyForchheimer was introduced in [22] and a blockcentered finite difference method was considered in [36, 35].
Developing robust numerical schemes on general meshes has drawn great attention in recent years, and numerous numerical schemes have been developed on general meshes. Among all the methods, we mention in particular the virtual element (VEM) method, the hybrid high order (HHO) method and the weak Galerkin (WG) method [3, 17, 38, 8]. However, to the best of our knowledge, none of the existing numerical schemes applicable on general meshes have been exploited to simulate the coupling of the Stokes and DarcyForchheimer problems. In addition, very few results are available for coupling of the Stokes and DarcyForchheimer problems even on triangular meshes, and it is not an easy task to design numerical schemes for this coupled model. In [39] a stabilized CrouzeixRaviart element method is developed and optimal order error estimates are achieved. In [1] a fully mixed method is proposed, whereas the authors can only obtain a suboptimal convergence rate. Staggered discontinuous Gaerkin (DG) method is initially developed by Chung and Engquist [11, 12]
, since then it has been successfully applied to a wide range of partial differential equations arising from practical applications, see, e.g.,
[13, 14, 26, 25, 15, 40, 42]. Recently, staggered DG method has been successfully designed for Darcy law and the Stokes equations, respectively, on general quadrilateral and polygonal meshes [41, 43], and the numerical results proposed therein indicate that it is accurate and robust for rough grids. Staggered DG method has many distinctive features such as local and global conservations, optimal convergence estimates and robustness to mesh distortion. Furthermore, the method favors adaptive mesh refinement by allowing hanging nodes. Thus, staggered DG method is a good candidate for problems under real applications, and it will be beneficial if one can apply this method to physical problems. Therefore, the goal of this paper is to initiate a study on staggered DG method for the coupled Stokes and DarcyForchheimer problems that can be flexibly applied to general quadrilateral and polygonal meshes, in addition, arbitrary polynomial orders are allowed.The key idea of staggered DG method is to divide the initial partition (general quadrilateral or polygonal meshes) into the union of triangles by connecting the interior point to all the vertices of the initial partition. Then the corresponding basis functions for the associated variables with staggered continuity are defined on the resulting triangulations. The interface conditions are imposed by switching the roles of the variables met on the interface. By doing so, no additional variables are introduced on the interface, which is different from the method proposed in [1]. The resulting method allows nonmatching grids on the interface due to the special inclusion of the interface conditions. In addition to the aforementioned features, another distinctive property is that the staggered continuity property naturally gives an interelement flux term in contrast to other DG methods, where the numerical fluxes or penalty parameters have to be carefully chosen. Note that developing stable numerical schemes for the Stokes and DarcyForchheimer problems is still in its infancy, and the present method is well suited to general quadrilateral and polygonal meshes, which will definitely inspire more works in this direction.
The primary difficulty for the convergence estimates results from the terms appearing on the interface. To overcome this issue, a novel discrete trace inequality and a generalized PoincaréFriedrichs inequality are proposed. Our analysis can be carried out by exploiting these two inequalities and the monotone property of the nonlinear operator under certain regularity assumptions without any restrictions on the source terms. The main novelties of this paper are as follows: First, this is the first result for coupling of the Stokes and DarcyForchheimer problems which can be flexibly applied to general quadrilateral and polygonal meshes and well adapted to incorporate hanging nodes in the construction of the method. Second, we can prove the optimal convergence estimates by proving some new discrete trace inequality and generalized PoincaréFriedrichs inequality.
The rest of the paper is organized as follows. In the next section we derive staggered DG method for the coupled Stokes and DarcyForchheimer problems and present some preliminary lemmas. Then in Section 3 we present the convergence estimates, where the discrete trace inequality and the generalized PoincaréFriedrichs inequality are the crucial ingredients. Several numerical experiments are carried out in Section 4 to validate the theoretical results. Finally, the conclusion is given in Section 5.
2 Description of staggered DG method
2.1 Preliminaries
Let and be two bounded and simply connected polygonal domains in such that and . Then, let , and denote by
the unit normal vector pointing from
to and by the unit normal vector pointing from to , on the interface we have . In addition, represents the unit tangential vector along the interface . Figure 1 gives a schematic representation of the geometry.When kinematic effects surpass viscous effects in a porous medium, the Darcy velocity and the pressure gradient does not satisfy a linear relation. Instead, a nonlinear approximation, known as the DarcyForchheimer model, is considered. When it is imposed on the porous medium with homogeneous Dircihlet boundary condition on the equations read:
(2.1)  
(2.2)  
(2.3) 
where
is the permeability tensor, assumed to be uniformly positive definite and bounded,
is the density of the fluid, is its viscosity and is a dynamic viscosity, all assumed to be positive constants. In addition, and are source terms. We remark that in this context we exploit homogeneous Dirichlet boundary condition, in fact, we can also consider homogeneous Neumann boundary condition, i.e., and the arguments used in this paper are still true.The fluid motion in is described by the Stokes equations:
(2.4)  
(2.5)  
(2.6)  
(2.7) 
where denotes the viscosity of the fluid.
On the interface, we prescribe the following interface conditions
(2.8a)  
(2.8b)  
(2.8c) 
Condition (2.8a) represents continuity of the fluid velocity’s normal components, (2.8b) represents the balance of forces acting across the interface, and (2.8c) is the BeaverJosephSaffman condition [2]. The constant is given and is usually obtained from experimental data.
Before closing this section, we introduce some notations that will be used later. Let . By , we denote the inner product in . We use the same symbol for the scalar product in and in . We denote by the inner product in , and its vector and tensor versions. Given an integer and , and denote the usual Sobolev space provided the norm and seminorm , . If we usually write and , and . In the sequel, we use to denote a generic positive constant which may have different values at different occurrences.
2.2 Weak formulation
This section presents the derivation of the weak formulation for the model problems (2.1)(2.8c). To begin, set
Then multiplying (2.4) by , (2.1) by , integration by parts and adding these two equations imply
Simple algebraic calculation yields
which gives by employing (2.8b) and (2.8c)
(2.9) 
Here, is the identity matrix.
Therefore, the weak formulation for the coupled Stokes and DarcyForchheimer problems reads: find and such that
(2.10) 
2.3 Description of staggered DG method
This section gives a detailed derivation of staggered DG method for (2.1)(2.8c). Then the infsup condition for the bilinear forms related to the DarcyForchheimer region is proved motivated by the methodology employed in [12]. First, we introduce the staggered meshes and the spaces exploited in the definition of staggered DG method. Following [41, 43], we first let () be the initial partition of the domain into nonoverlapping quadrilateral or polygonal meshes. We require that be aligned with . We also let be the set of all edges in the initial partion and be the subset of all interior edges of . The union of all the nodal points in the initial partition is denoted as and . For each quadrilateral or polygonal mesh in the initial partition , we select an interior point and create new edges by connecting to the vertices of the primal element. This process will divide into the union of triangles, where the triangle is denoted as , and we rename the union of these triangles by . The union of all is denoted as and . We remark that is the quadrilateral or polygonal mesh in the initial partition. Moreover, we will use to denote the set of all the new edges generated by this subdivision process and use to denote the resulting quasiuniform triangulation, on which our basis functions are defined. In addition, we define and . For each triangle , we let be the diameter of and , and we define . Also, we let to denote the length of edge . This construction is illustrated in Figure 2, where black solid lines are edges in and red dotted lines are edges in .
Moreover, the union of all for is denoted as . For each interior edge , we use to denote the union of the two triangles in sharing the edge , and for each boundary 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.
Assumption 2.1.
We assume the existence of a constant such that

For every element and every edge , it satisfies , where denotes the diameter of .

Each element in is starshaped with respect to a ball of radius .
We remark that the above assumptions ensure that the triangulation is shape regular.
For a scalar or vector function belonging to broken Sobolev space, its jump on is defined as
where and , are the two triangles in having the edge . For the boundary edges, we simply define .
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. Next, we will introduce the staggered DG spaces. First, we define the following locally conforming staggered DG space for
Notice that, if , then for each edge
. The degrees of freedom for this space can be described as:
(UD1) For each edge , we have
(UD2) For each , we have
The discrete norm for the space is given by
We define the following norm for
Due to the nonlinear term in the DarcyForchheimer region, we adjust the commonly employed norms for staggered DG space, instead, we introduce the following meshdependent norms for , which flavors our analysis
We next define the following staggered DG space
Note that if , then for each . The following degrees of freedom can be defined correspondingly.
(VD1) For each edge , we have
(VD2) For each , we have
We define the following discrete norm and discrete seminorm for
In addition, is equipped by
Scaling arguments imply there exists a positive constant such that
(2.11) 
Finally, locally conforming finite element space for pressure is defined as:
with norm
Norm equivalence yields
Now we are ready to derive the discrete formulation for (2.1)(2.8). Multiplying (2.4) by , (2.5) by , (2.6) by and integration by parts yield
(2.12) 
Multiplying (2.1) by and (2.2) by , it then follows from integration by parts
(2.13) 
Adding the first equation of (2.12) and (2.13), and employing (2.9), we have
Adding the second equations of (2.12) and (2.13), we can get
For the convenience of the presentation, we introduce the nonlinear operator defined by
(2.14) 
In addition, we define the following bilinear forms corresponding to the aforementioned derivations
Then, the staggered DG formulation for the coupled Stokes and DarcyForchheimer problems reads: find and such that
(2.15) 
Next, we define the following projection operators, which are defined by employing the special property of staggered DG method. Let be defined by
In addition, is defined by
Straightforwardly, we have
Comments
There are no comments yet.