# A uniformly robust staggered DG method for the unsteady Darcy-Forchheimer-Brinkman problem

In this paper we propose and analyze a uniformly robust staggered DG method for the unsteady Darcy-Forchheimer-Brinkman problem. Our formulation is based on velocity gradient-velocity-pressure and the resulting scheme can be flexibly applied to fairly general polygonal meshes. We relax the tangential continuity for velocity, which is the key ingredient in achieving the uniform robustness. We present well-posedness and error analysis for both the semi-discrete scheme and the fully discrete scheme, and the theories indicate that the error estimates for velocity are independent of pressure. Several numerical experiments are presented to confirm the theoretical findings.

## Authors

• 12 publications
• 3 publications
• 27 publications
• ### A pressure robust staggered discontinuous Galerkin method for the Stokes equations

In this paper we propose a pressure robust staggered discontinuous Galer...
07/01/2020 ∙ by Lina Zhao, et al. ∙ 0

• ### A new staggered DG method for the Brinkman problem robust in the Darcy and Stokes limits

In this paper we propose a novel staggered discontinuous Galerkin method...
11/20/2019 ∙ by Lina Zhao, et al. ∙ 0

• ### Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows

Proper orthogonal decomposition (POD) stabilized methods for the Navier-...
05/30/2020 ∙ by Julia Novo, et al. ∙ 0

• ### Locking free staggered DG method for the Biot system of poroelasticity on general polygonal meshes

In this paper we propose and analyze a staggered discontinuous Galerkin ...
11/10/2020 ∙ by Lina Zhao, et al. ∙ 0

• ### Adaptive staggered DG method for Darcy flows in fractured porous media

Modeling flows in fractured porous media is important in applications. O...
06/08/2020 ∙ by Lina Zhao, et al. ∙ 0

• ### Guaranteed upper bounds for the velocity error of pressure-robust Stokes discretisations

This paper improves guaranteed error control for the Stokes problem with...
08/13/2020 ∙ by Philip L. Lederer, et al. ∙ 0

• ### Dynamics of fracturing saturated porous media and self-organization of rupture

Analytical solutions and a vast majority of numerical ones for fracture ...
02/26/2019 ∙ by C. Peruzzo, et al. ∙ 0

##### 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 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 [15] 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.

Staggered discontinuous Galerkin method as a new generation of numerical schemes is initially proposed by Chung and Engquist to solve wave propagation problems [9, 10]

, 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 [31]. 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 [31] 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:

 ut−ϵΔu+αu+β|u|u+∇p =f in Ω,t∈[0,T], (2.1) divu =0 in Ω,t∈[0,T], u =0 on ∂Ω,t∈[0,T], u(x,0) =u0 in Ω,

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:

 L−√ϵ∇u =0 in Ω,t∈[0,T], (2.2) ut−√ϵdivL+αu+β|u|u+∇p =f in Ω,t∈[0,T], divu =0 in Ω,t∈[0,T], u =0 on ∂Ω,t∈[0,T], u(x,0) =u0 in Ω.

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 in

(or duality pairing), for a scalar, vector, or tensor functions. Given an integer

and , 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

 ∥ψ∥L2(0,T;Hm(Ω))=(∫T0∥ψ(t)∥2Hm(Ω)dt)1/2,∥ψ∥C(0,T;Hm(Ω))=max0≤t≤T∥ψ(t)∥Hm(Ω).

Integration by parts yields the following weak formulation: Find such that

 (L,G)−√ϵ(∇u,G) =0 ∀G∈L3/2(Ω)2×2, (2.3) (ut,v)+√ϵ(L,∇v)+(αu,v)+(β|u|u,v)−(p,divv) =(f,v) ∀v∈W1,30(Ω)2, (divu,q) =0 ∀q∈L3/2(Ω),

where

 W1,30(Ω)2={v∈W1,3(Ω)2,v=0on∂Ω}.

The bilinear forms given above satisfy the following inf-sup condition:

 infv∈W1,30(Ω)2(∇v,G)∥G∥L3/2(Ω)∥v∥W1,3(Ω)≥C (2.4)

and (cf. [1])

 infq∈L2(Ω)supv∈W1,30(Ω)2(∇⋅v,q)Ω∥q∥L3/2(Ω)∥v∥W1,3(Ω)≥C. (2.5)

To ease later analysis, we define for any

 N(v)=αv+β|v|v. (2.6)

We describe some properties for , which will play an important role for later analysis. One can refer to [16] for more details. First, we have

 |N(v)−N(w)|≤αmax|v−w|+βmax|v−w|(|v|+|w|),∀v,w∈L3(Ω)2. (2.7)
###### Lemma 2.1.

For fixed , the mapping defined by (2.6) is monotone from into :

 ∀u,v∈L3(Ω)2,∫Ω(N(u+uℓ)−N(v+uℓ))⋅(u−v)dx≥αmin∥u−v∥2L2(Ω). (2.8)
###### Lemma 2.2.

For fixed , the mapping defined by (2.6) is coercive in

 lim∥u∥L3(Ω)→∞(1∥u∥L3(Ω)∫ΩN(u+ul)⋅udx)=∞.
###### Lemma 2.3.

The mapping is hemi-continuous in ; for fixed and in , the mapping

 γ→∫ΩN(ul+u+γv)⋅vdx

is continuous from into .

###### Theorem 2.1.

There exists a unique solution to (2.3). In addition, the following estimate holds

 ∫t0(∥L∥2L2(Ω)+αmin2∥u∥2L2(Ω))ds+12∥u(t)∥2L2(Ω)≤∫t012α% min∥f∥2L2(Ω)ds. (2.9)
###### Proof.

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

 ∥L1−L2∥2L2(Ω)+12∂t∥u1−u2∥2L2(Ω)+(N(u1)−N(u2),u1−u2)=0,

which yields

 ∥L1−L2∥2L2(Ω)+12∂t∥u1−u2∥2L2(Ω)+∥u1−u2∥2L2(Ω)≤0.

Integrating in time from to and using , we obtain

 12∥u1(t)−u2(t)∥2L2(Ω)+∫t0(∥L1−L2∥2L2(Ω)+∥u1−u2∥2L2(Ω))≤0.

Therefore, we can infer that and , which implies that (2.3) has a unique solution.

Next, we will show the stability estimate (2.9). Taking , and in (2.3), we can get

 ∥L∥2L2(Ω)+12ddt∥u∥2L2(Ω)+(N(u),u)=(f,u)≤12αmin∥f∥2L2(Ω)+αmin2∥u∥2L2(Ω).

It follows from the definition of that

 (N(u),u)≥αmin∥u∥2L2(Ω),

thereby we can infer that

 ∥L∥2L2(Ω)+12ddt∥u∥2L2(Ω)+αmin∥u∥2L2(Ω)≤12αmin∥f∥2L2(Ω)+αmin2∥u∥2L2(Ω),

which yields

 ∥L∥2L2(Ω)+12ddt∥u∥2L2(Ω)+αmin2∥u∥2L2(Ω)≤12αmin∥f∥2L2(Ω).

Integrating over time and using the fact that imply

 ∫t0(∥L∥2L2(Ω)+α% min2∥u∥2L2(Ω))ds+12∥u(t)∥2L2(Ω)≤∫t012αmin∥f∥2L2(Ω)ds.

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 :

 ¯Ω=⋃E∈TuE.

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 [3]. 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

 δi=ni⋅n={1 if ni=n on e,−1 if ni=−n on e.

For a double-valued scalar quantity , let .The jump across an edge can then be defined as:

 [ϕ]=δ1ϕ1+δ2ϕ2,

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:

 [ϕ⋅n]=δ1(ϕ1⋅n)+δ2(ϕ2⋅n),[Φn]=δ1(Φ1n)+δ2(Φ2n). (3.1)

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 [31]. We first define the following finite element space for velocity:

 Uh={v:v|τ∈Pk(τ)2;τ∈Th;v⋅n is continuous over e∈Fp}.

In this space, we define

 ∥v∥3Z2 =∥∇hv∥3L3(Ω)+∑e∈Fuh−2e∥[v]∥3L3(e)+∑e∈Fph−2e∥[(v⋅t)t]∥3L3(e).

Next, we define the following finite element space for velocity gradient:

 Wh={G:G|τ∈Pk(τ)2×2;τ∈Th;Gn is continuous over e∈F0u}.

Then, we define the following locally -conforming finite element space for pressure:

 Ph={q:q|τ∈Pk(τ);τ∈Th;q is % continuous over e∈F0u;∫Ωqdx=0},

which is equipped by

 ∥q∥3/23/2,h =∑τ∈Th∫τ|∇q|3/2dx+∑e∈Fph−1/2e∫e|[q]|3/2ds.

Finally, we define the following space which is employed to enforce the weak continuity of the velocity gradient over the dual edge

 ˆUh={ˆv:ˆv|e∈Pk(e)2;ˆv⋅n∣e=0∀e∈Fp}.

Following [31], we can formulate our staggered DG formulation for the unsteady Darcy-Forchheimer-Brinkman problem (2.1): Find such that

 (Lh,G)=√ϵB∗h(uh,G)+√ϵT∗h(ˆuh,G)∀G∈Wh,(∂tuh,v)+√ϵBh(Lh,v)+(N(uh),v)+b∗h(ph,v)=(f,v)∀v∈Uh,−bh(uh,q)=0∀q∈Ph,Th(Lh,ˆv)=0∀ˆv∈ˆUh, (3.2)

where the bilinear forms are defined by

 B∗h(v,G) =−∫Ωv⋅divhGdx+∑e∈Fp∫e(v⋅n)n⋅[Gn]ds, Bh(G,v) =∫ΩG⋅∇hvdx−∑e∈Fu∫e[v]⋅(Gn)ds−∑e∈Fp∫e[(v⋅t)t⋅(Gn)]ds, T∗h(ˆv,G) =∑e∈Fp∫eˆv⋅[Gn]ds, Th(G,ˆv) =∑e∈Fp∫e[Gn]⋅ˆvds, b∗h(q,v) =−∫Ωqdivhvdx+∑e∈Fu∫eq[v⋅n]ds, bh(v,q) =∫Ωv⋅∇hqdx−∑e∈Fp∫ev⋅n[q]ds.

Performing integration by parts reveals the following adjoint properties

 Bh(G,v)=B∗h(v,G)∀(G,v)∈Wh×Uh,bh(v,q)=b∗h(q,v)∀(v,q)∈Uh×Ph,Th(G,ˆv)=T∗h(ˆv,G)∀(G,ˆv)∈Wh×ˆUh. (3.3)

To facilitate the analysis, we define the subspace of by

 ˆWh:={G∈Wh:∫e[Gn]⋅^vds=0∀^v∈ˆUh,∀e∈Fp}.

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

 (Lh,G)=√ϵB∗h(uh,G),(∂tuh,v)+√ϵBh(Lh,v)+b∗h(ph,v)+(N(uh),v)=(f,v),−bh(uh,q)=0 (3.4)

for all .

For later analysis, we state the following inf-sup condition

 infq∈Phsupv∈Uhbh(v,q)∥v∥L3(Ω)|||q|||3/2,h≥C (3.5)

and

 infv∈UhsupG∈ˆWhBh(G,v)∥G∥L3/2(Ω)∥v∥Z2≥C. (3.6)

We remark that the proof of the above inf-sup conditions can follow the techniques employed in [31, 32] and we omit it for simplicity.

The inf-sup condition (3.6

) implies the existence of an interpolation operator

such that

 Bh(L−ΠhL,v)=0∀v∈Uh. (3.7)

By the standard theory for polynomial preserving operators (cf. [13, 10]), we obtain

 ∥L−ΠhL∥L2(Ω) ≤Chk+1|L|Hk+1(Ω). (3.8)

To facilitate later analysis, we also need to define the following projection operator. Let be defined by

 (Ihq−q,ϕ)e =0∀ϕ∈Pk(e),∀e∈Fu, (Ihq−q,ϕ)τ =0∀ϕ∈Pk−1(τ),∀τ∈Th

and let be defined by

 ((Jhv−v)⋅n,φ)e =0∀φ∈Pk(e), ∀e∈Fp, (Jhv−v,ϕ)τ =0∀ϕ∈Pk−1(τ)2, ∀τ∈Th.

It is easy to see that and are well defined polynomial preserving operators. In addition, the following approximation properties hold for (cf. [13, 10])

 ∥v−Jhv∥L2(Ω) ≤Chk+1|v|Hk+1(Ω), (3.9) ∥v−Jhv∥L4(Ω) ≤Chk+1|v|Wk+1,4(Ω). (3.10)

Furthermore, the definitions of and imply directly that

 bh(v−Jhv,q) =0∀q∈Ph, (3.11) b∗h(q−Ihq,v) =0∀v∈Uh. (3.12)

## 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.

###### Lemma 4.1.

There exists a unique solution to (3.4), in addition, the following estimate holds

 ∫t0(∥Lh∥2L2(Ω)+α% min∥uh∥2L2(Ω))ds+12∥uh(t)∥2L2(Ω)≤C∫t0∥f∥2L2(Ω)ds. (4.1)
###### Proof.

We can proceed similarly to Theorem 2.1 to infer that there exists a unique solution to (3.4).

Now we show the stability estimate (4.1). Taking , and in (3.4), and summing up the resulting equations yields

 ∥Lh∥2L2(Ω)+(∂tuh,uh)+(N(uh),uh)≤∥f∥L2(Ω)∥uh∥L2(Ω). (4.2)

According to the definition of , we can easily see that

 (N(uh),uh)≥αmin∥uh∥2L2(Ω).

Thereby we can infer from (4.2) that

 ∥Lh∥2L2(Ω)+12ddt∥uh∥2L2(Ω)+αmin∥uh∥2L2(Ω)≤12αmin∥f∥2L2(Ω)+αmin2∥uh∥2L2(Ω).

Integrating over time and using the fact that yield

 ∫t0(∥Lh∥2L2(Ω)+αmin∥uh∥2L2(Ω))+12∥uh(t)∥2L2(Ω)≤C∫t0∥f∥2L2(Ω).

Therefore, the proof is completed.

###### Lemma 4.2.

Let be the weak solution of (2.3) and be the numerical solution of (3.4), then the following identity holds

 (4.3)
###### Proof.

Replacing by in (3.4) yields the following error equations

 (L−Lh,G)=√ϵB∗h(u−uh,G),(∂t(u−uh),v)+√ϵBh(L−Lh,v)+b∗h(p−ph,v)+(N(u)−N(uh),v)=0,bh(u−uh,q)=0 (4.4)

for all .

Taking , and in (4.4) and adding the resulting equations, then we can infer from (3.7), (3.11) and (3.12) that

 (∂t(u−uh),Jhu−uh)+(L−Lh,ΠhL−Lh)+(N(