The Brinkman problem can model fluid motion in porous media with fractures. A distinctive feature of Brinkman model is that it can behave like Stokes or Darcy problem by tuning the parameters related to the fluid viscosity and the medium permeability, respectively. But this at the same time brings additional technical difficulties for the design of numerical schemes that are robust in both Stokes and Darcy regimes. In general, the traditional Stokes stable elements will suffer from loss of convergence when the Brinkman problem becomes Darcy dominating, vice versa. To overcome this issue, numerous numerical schemes have been developed for the Birnkman problem [22, 24, 4, 5, 2, 17]. Darcy’s law is widely employed to model flow in porous media, becomes unreliable for Reynolds numbers greater than one. The Forchheimer model  accounts for faster flows by including a nonlinear inertial term and has been extensively studied [20, 26, 16, 25, 29, 28]. The Brinkman-Forchheimer model [27, 7, 23] combines the advantages of the two models and can be used for fast flows in highly porous media.
, and inspires study for various partial differential equations arising from practical applications[8, 19, 18, 12, 11, 14]. Recently, the concept of staggered DG method is integrated into polygonal methods and have been successfully designed for numerous mathematical models that have important physical applications [33, 35, 36, 31, 32, 21, 37, 34]. A uniformly stable staggered DG method has been established for the Brinkman problem by relaxing the tangential continuity of velocity in . The numerical experiments presented therein verify that the proposed scheme is robust with respect to viscosity and the accuracy of velocity remains almost the same for various values of viscosity. This feature is desirable in practical applications, therefore, the purpose of this paper is to extend the staggered DG method developed in  to the unsteady Darcy-Forchheimer-Brinkman problem.
Our formulation is based on velocity gradient-velocity-pressure, and the finite element spaces for these three variables enjoy staggered continuity properties. Specifically, the finite element space for velocity is continuous in the normal direction over the dual edges and the finite element space for pressure is continuous over the primal edges. Thus, our approach can be viewed as Darcy tailored method. We also emphasize that relaxing tangential continuity for velocity is the crux in the design of our uniformly stable scheme and the spaces exploited in [33, 34] will lead to loss of convergence when viscosity approaches zero. There are several desirable features of this approach, which can be summarized as follows: First, our method can be flexibly applied to fairly general meshes possibly including hanging nodes. Second, all the variables of physical interest can be calculated simultaneously. Third, it is uniformly robust with respect to the Brinkman parameter in both the Stokes and Darcy regimes, and the accuracy of velocity remains almost the same with various values of the Brinkman parameter. In this paper, we analyze the convergence for both the semi-discrete scheme and the fully discrete scheme, where backward Euler is employed for the time discretization. We can achieve optimal rates of convergence, and the convergence of velocity is independent of the Brinkman parameter. Finally, we perform several numerical experiments to verify the proposed theories and we can observe that our method is uniformly robust with respect to the parameters, as expected, the accuracy of velocity remains almost the same for various values of the parameters.
The rest of the paper is organized as follows. In the next section, we present the continuous weak formulation, in addition, the unique solvability and stability is proved. In section 3, we describe the discrete formulation for the unsteady Darcy-Forchheimer-Brinkman problem. Then unique solvability of the discrete formulation and the error estimates are presented in section 4. Several numerical experiments are carried out in section 5 to illustrate that the proposed method is uniformly robust with respect to the parameters. Finally, a conclusion is given.
2 The continuous formulation
We consider the following unsteady Darcy-Forchheimer-Brinkman model:
where is the Darcy coefficient, is the Brinkman coefficient, is the Forchheimer coefficient and is the external body force. In addition, we also assume that there exist real numbers , and such that and . The unknowns are the velocity and the pressure , which are functions of the spatial variable and temporal variable in ( is a finite time). Here for .
We introduce an additional unknown , then the above system of equations can be recast into the following first order system of equations:
Here for simplicity we assume for the subsequent analysis. We introduce some notations that will be used later. For a set , we denote the scalar product in by , namely , we use the same symbol for the inner product in and in . More precisely for . When coincides with , the subscript will be dropped unless otherwise mentioned. We denote by the scalar product inand , and denote the usual Sobolev space provided the norm and semi-norm , . If we usually write and , and . In addition, we need spaces of vector values functions such as and with the norms
Integration by parts yields the following weak formulation: Find such that
The bilinear forms given above satisfy the following inf-sup condition:
and (cf. )
To ease later analysis, we define for any
We describe some properties for , which will play an important role for later analysis. One can refer to  for more details. First, we have
For fixed , the mapping defined by (2.6) is monotone from into :
For fixed , the mapping defined by (2.6) is coercive in
The mapping is hemi-continuous in ; for fixed and in , the mapping
is continuous from into .
There exists a unique solution to (2.3). In addition, the following estimate holds
Since is monotone, coercive and hemi-continuous (cf. Lemmas 2.1-2.3), in addition, the inf-sup condition (2.4) and (2.5) hold, we can follow [30, 6] to show that there exists a solution to (2.3) and we omit the proof for simplicity. Next, we show that the solution is unique. Assume that the solution of (2.3) is not unique. Let with be two solutions corresponding to the same data. Then, taking (2.3) with and summing up the resulting equations, we can infer that
Integrating in time from to and using , we obtain
Therefore, we can infer that and , which implies that (2.3) has a unique solution.
It follows from the definition of that
thereby we can infer that
Integrating over time and using the fact that imply
Therefore, the proof is completed.
3 Description of staggered DG method
In this section, we introduce the discrete formulation for the unsteady Darcy-Forchheimer-Brinkman problem (2.1). To this end, we first introduce the construction of our staggered DG spaces, in line with this we then present the construction of staggered DG method. To begin, we construct three meshes: the primal mesh , the dual mesh , and the primal simplicial submeshes . For a polygonal domain , consider a general mesh (of ) that consists of nonempty connected close disjoint subsets of :
We let be the set of all primal edges in this partition and be the subset of all interior edges, that is, the set of edges in that do not lie on . We construct the primal submeshes as a triangular subgrid of the primal grid: for an element , elements of are obtained by connecting the interior point to all vertices of (see Figure 1). We use to denote the set of all the dual edges generated by this subdivision process. For each triangle , we let be the diameter of , be the length of edge , and . In addition, we define and . We rename the primal element by , which is assumed to be star-shaped with respect to a ball of radius , where is a positive constant. In addition, we assume that for every edge , it satisfies , where is a positive constant. More discussions about mesh regularity assumptions for general meshes can be referred to . The construction for general meshes is illustrated in Figure 1, where the black solid lines are edges in and the red dotted lines are edges in .
Finally, we construct the dual mesh. For each interior edge , we use to denote the dual mesh, which is 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 1.
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. In addition, we use to denote the corresponding unit tangent vector. We also introduce some notations that will be employed throughout this paper. 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. In the following, we use and to denote the element-wise gradient and divergence operators, respectively.
We now define jump terms which will be used throughout the paper. For each triangle in such that , we let be the outward unit normal vector on . The sign of with respect to on is then given by
For a double-valued scalar quantity , let .The jump across an edge can then be defined as:
where and are the two triangles sharing the common edge . For , we let . Similarly, for a vector quantity and a matrix quantity , we let and , then the jumps and across an edge are defined as:
In addition, for , we define and . In the sequel, we use to denote a positive constant which may have different values at different occurrences.
Now we are ready to introduce the finite dimensional spaces for our staggered DG method by following . We first define the following finite element space for velocity:
In this space, we define
Next, we define the following finite element space for velocity gradient:
Then, we define the following locally -conforming finite element space for pressure:
which is equipped by
Finally, we define the following space which is employed to enforce the weak continuity of the velocity gradient over the dual edge
where the bilinear forms are defined by
Performing integration by parts reveals the following adjoint properties
To facilitate the analysis, we define the subspace of by
Based on the definition of and the discrete formulation (3.2), we can conclude that . Therefore, we can reformulate our discrete formulation (3.2) and obtain the following equivalent formulation: Find such that
for all .
For later analysis, we state the following inf-sup condition
The inf-sup condition (3.6
) implies the existence of an interpolation operatorsuch that
To facilitate later analysis, we also need to define the following projection operator. Let be defined by
and let be defined by
Furthermore, the definitions of and imply directly that
4 Error analysis
In this section, we first perform error analysis and obtain rates of convergence for the semi-discrete scheme. Then we introduce the fully discrete scheme by using backward Euler for the time discretization, and a error analysis is established for the resulting fully discrete scheme.
4.1 Error analysis for semi-discrete scheme
In this subsection we establish the unique solvability and convergence estimates for the semi-discrete scheme. To this end, we first introduce the following lemma, which states the unique solvability and stability.
There exists a unique solution to (3.4), in addition, the following estimate holds
According to the definition of , we can easily see that
Thereby we can infer from (4.2) that
Integrating over time and using the fact that yield
Therefore, the proof is completed.