Staggered DG method for coupling of the Stokes and Darcy-Forchheimer problems

06/16/2019 ∙ by Lina Zhao, et al. ∙ Yonsei University The Chinese University of Hong Kong 0

In this paper we develop a staggered discontinuous Galerkin method for the Stokes and Darcy-Forchheimer problems coupled with the Beavers-Joseph-Saffman conditions. The method is defined by imposing staggered continuity for all the variables involved and the interface conditions are enforced by switching the roles of the variables met on the interface, which eliminate the hassle of introducing additional variables. This method can be flexibly applied to rough grids such as the highly distorted grids and the polygonal grids. In addition, the method allows nonmatching grids on the interface thanks to the special inclusion of the interface conditions, which is highly appreciated from a practical point of view. A new discrete trace inequality and a generalized Poincaré-Friedrichs inequality are proved, which enables us to prove the optimal convergence estimates under reasonable regularity assumptions. Finally, several numerical experiments are given to illustrate the performances of the proposed method, and the numerical results indicate that the proposed method is accurate and efficient, in addition, it is a good candidate for practical applications.



There are no comments yet.


page 22

page 24

page 26

page 27

page 28

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The fluid flow between porous media and free-flow 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 Beavers-Joseph-Saffman 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 Darcy-Forchheimer law accounts for this nonlinear behavior. Mixed methods for the steady Darcy-Forchheimer flow was studied in more general setting in [31] and unsteady case was analyzed in [27, 32]. A nonconforming primal mixed method for the Darcy-Forchheimer was introduced in [22] and a block-centered 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 Darcy-Forchheimer problems. In addition, very few results are available for coupling of the Stokes and Darcy-Forchheimer problems even on triangular meshes, and it is not an easy task to design numerical schemes for this coupled model. In [39] a stabilized Crouzeix-Raviart 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 Darcy-Forchheimer 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 Darcy-Forchheimer 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 Darcy-Forchheimer 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 Darcy-Forchheimer 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 Darcy-Forchheimer model, is considered. When it is imposed on the porous medium with homogeneous Dircihlet boundary condition on the equations read:



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:


where denotes the viscosity of the fluid.

On the interface, we prescribe the following interface conditions


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 Beaver-Joseph-Saffman 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 semi-norm , . 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.

Figure 1: Coupled domain with interface .

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)


Here, is the identity matrix.

Therefore, the weak formulation for the coupled Stokes and Darcy-Forchheimer problems reads: find and such that


2.3 Description of staggered DG method

Figure 2: Schematic of the primal mesh , the dual mesh and the primal simplexes.

This section gives a detailed derivation of staggered DG method for (2.1)-(2.8c). Then the inf-sup condition for the bilinear forms related to the Darcy-Forchheimer 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 non-overlapping 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 quasi-uniform 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.

We also make the following mesh regularity assumptions (cf. [3, 10, 41]).

Assumption 2.1.

We assume the existence of a constant such that

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

  2. Each element in is star-shaped 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 Darcy-Forchheimer region, we adjust the commonly employed norms for staggered DG space, instead, we introduce the following mesh-dependent 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


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


Multiplying (2.1) by and (2.2) by , it then follows from integration by parts


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


In addition, we define the following bilinear forms corresponding to the aforementioned derivations

Then, the staggered DG formulation for the coupled Stokes and Darcy-Forchheimer problems reads: find and such that


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