# A mixed discontinuous Galerkin method with symmetric stress for Brinkman problem based on the velocity-pseudostress formulation

The Brinkman equations can be regarded as a combination of the Stokes and Darcy equations which model transitions between the fast flow in channels (governed by Stokes equations) and the slow flow in porous media (governed by Darcy's law). The numerical challenge for this model is the designing of a numerical scheme which is stable for both the Stokes-dominated (high permeability) and the Darcy-dominated (low permeability) equations. In this paper, we solve the Brinkman model in n dimensions (n = 2, 3) by using the mixed discontinuous Galerkin (MDG) method, which meets this challenge. This MDG method is based on the pseudostress-velocity formulation and uses a discontinuous piecewise polynomial pair P_k+1^S-P_k (k ≥ 0), where the stress field is symmetric. The main unknowns are the pseudostress and the velocity, whereas the pressure is easily recovered through a simple postprocessing. A key step in the analysis is to establish the parameter-robust inf-sup stability through specific parameter-dependent norms at both continuous and discrete levels. Therefore, the stability results presented here are uniform with respect to the permeability. Thanks to the parameter-robust stability analysis, we obtain optimal error estimates for the stress in broken H( div)-norm and velocity in L^2-norm. Furthermore, the optimal L^2 error estimate for pseudostress is derived under certain conditions. Finally, numerical experiments are provided to support the theoretical results and to show the robustness, accuracy, and flexibility of the MDG method.

## Authors

• 2 publications
• 9 publications
• 127 publications
12/07/2021

### Improved a priori error estimates for a discontinuous Galerkin pressure correction scheme for the Navier-Stokes equations

The pressure correction scheme is combined with interior penalty discont...
11/12/2019

### Tutorial on Hybridizable Discontinuous Galerkin (HDG) Formulation for Incompressible Flow Problems

A hybridizable discontinuous Galerkin (HDG) formulation of the linearize...
06/21/2021

### Energy stability analysis of turbulent incompressible flow based on the triple decomposition of the velocity gradient tensor

In the context of flow visualization a triple decomposition of the veloc...
02/12/2020

### New analysis of Galerkin-mixed FEMs for incompressible miscible flow in porous media

Analysis of Galerkin-mixed FEMs for incompressible miscible flow in poro...
12/11/2021

### A strongly mass conservative method for the coupled Brinkman-Darcy flow and transport

In this paper, a strongly mass conservative and stabilizer free scheme i...
08/21/2020

### A unified framework of continuous and discontinuous Galerkin methods for solving the incompressible Navier–Stokes equation

In this paper, we propose a unified numerical framework for the time-dep...
12/15/2020

### Uniformly well-posed hybridized discontinuous Galerkin/hybrid mixed discretizations for Biot's consolidation model

We consider the quasi-static Biot's consolidation model in a three-field...
##### 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 equations (cf. Brinkman (1949)),

 νκ––−1u−2νdiv(ε–(u))+∇p =f in Ω, (1.1a) divu =0 in Ω, (1.1b)

which model the flow of a fluid through a complex porous medium occupying domain

with a high-contrast permeability tensor

, can be seen as a mixture of Darcy and Stokes equations. Here, is the fluid viscosity, denotes the velocity field, is the strain rate, is the pressure field, and is the volume force. This model arises from applications in many fields, such as groundwater hydrology, biomedical engineering, petroleum industry, and environmental science (cf. Wehrspohn (2005); Ligaarden et al. (2010); Vafai (2010)). From (1.1), we can see that the Brinkman problem becomes Stokes-dominated when the permeability tensor is getting large, and it becomes Darcy-dominated when is quite small.

It is well known that the usual Darcy stable element pairs may diverge for Stokes flow and vice versa. Therefore, the numerical challenge for solving the Brinkman model is to construct a stable discretization method for both the Stokes and the Darcy equations. As shown in Mardal et al. (2002), the Darcy stable finite element pairs, for example, the Raviart-Thomas element (cf. Raviart and Thomas (1977)) leads to non-convergent results as the Brinkman model becomes Stokes-dominated; on the other hand, the usual Stokes stable finite element pairs, such as Mini-element (cf. Arnold et al. (1984a)), - element, nonconforming Crouzeix-Raviart (CR) finite element (cf. Crouzeix and Raviart (1973)), will diverge when the Brinkman equations turn into Darcy-dominated. Therefore, many researchers pay attention to developing stable and accurate numerical methods for solving Brinkman equations. One way to circumvent this difficulty is to modify the existing Stokes or Darcy elements to make them work well for the Brinkman model. In Burman and Hansbo (2005), inspired by discontinuous Galerkin (DG) method, a stabilized CR finite element method is constructed by adding a penalty term. In addition, a generalization of classical Mini-element is studied in Juntunen and Stenberg (2010); stabilized equal-order finite elements are proposed and analyzed in Braack and Schieweck (2011); (hybridized) interior penalty DG scheme with -conforming finite elements is investigated in Könnö and Stenberg (2011). Another approach is to develop new numerical schemes for solving Brinkman equations, for examples, pseudostress-based mixed finite element methods (cf. Barrios et al. (2012); Gatica et al. (2014)), weak Galerkin methods (cf. Mu et al. (2014); Zhai et al. (2016)), virtual element method (cf. Cáceres et al. (2017)), hybrid high-order method (cf. Botti et al. (2018)) and hybridizable discontinuous Galerkin method (cf. Gatica and Sequeira (2018)).

To study hydrodynamics, different formulations, like velocity-pressure, stress-velocity-pressure, pseudostress-velocity formulations, have been introduced and analyzed. The velocity-pressure formulation has been extensively studied in the computation of incompressible Newtonian flows (cf. Brezzi and Fortin (1991); Boffi et al. (2013)). However, the study of numerical methods for the stress-based and pseudostress-based formulations (cf. Farhloul and Fortin (1993); Behr et al. (1993); Cai et al. (2010); Gatica et al. (2010)) has become a very active research area because of the arising interest in non-Newtonian flows. The main advantage of the stress-based and pseudostress-based formulation is that it provides a unified framework for both the Newtonian and the non-Newtonian flows. In addition, physical quantity like the stress can be computed directly instead of by taking derivatives of the velocity, which avoids degrading of accuracy in the process of numerical differentiation. Precise computation of the stress is of paramount importance for the hydraulic fracturing problem as the crack propagation is determined by the stress field. While a formulation comprising the stress as a fundamental unknown is unavoidable for non-Newtonian flows in which the constitutive law is nonlinear, the drawback of the stress-velocity-pressure formulation is the increase in the number of unknowns. To avoid this disadvantage, we focus on the pseudostress-velocity formulation. Last but not least, we need to mention that the pressure field can be easily obtained by a simple postprocessing without affecting the accuracy of the approximation.

Due to the flexibility in constructing the local shape function spaces and the ability to capture non-smooth or oscillatory solutions effectively, DG methods have been applied to solve many problems in scientific computing and engineering, such as conservation laws (cf. Bey et al. (1995); Chen and Shu (2017)), Darcy flow (cf. Brezzi et al. (2005); Aizinger et al. (2018)), Navier-Stokes (or Stokes) equations (cf. Bassi and Rebay (1997); Kaya and Rivière (2005)), variational inequalities (cf. Wang et al. (2010); Brenner et al. (2012); Wang et al. (2014)) and much more. Besides, DG methods also enjoy the following advantages: (i) locally (and globally) conservative; (ii) easy to implement adaptivity; (iii) suitable for parallel computing. We refer to Cockburn et al. (2000); Brezzi et al. (2000); Arnold et al. (2002); Hong et al. (2019) for more discussion about DG methods.

In this paper, we construct a mixed discontinuous Galerkin (MDG) method with - element pair for solving the Brinkman equations based on the pseudostress-velocity formulation. The main results of this article include that: (i) The MDG scheme with symmetric stress field is uniformly stable and efficient for both Darcy-dominated and Stokes-dominated flows; (ii) Under specific parameter-dependent norms, the parameter-robust stability results of both continuous and discrete schemes are obtained; (iii) For , we get the optimal convergence order for the stress in broken -norm and velocity in -norm; (iv) When and the Stokes pair - is stable, we obtain the optimal error estimate for the pseudostress.

The rest of the paper is organized as follows. In Section 2, we introduce the pseudostress-velocity formulation for the Brinkman model and present some preliminary results. In Section 3, the MDG scheme is introduced and the well-posedness is obtained. We show the stability of the discrete scheme and prove optimal error estimates for both velocity and pressure in Section 4. In Section 5, numerical examples are provided to confirm the theoretical findings and to illustrate the performance of the mixed DG scheme. Finally, we give a short summary in Section 6.

## 2 Brinkman model in pseudostress-velocity formulation

In this section, we introduce the Brinkman model in the pseudostress-velocity formulation and provide the parameter-robust stability analysis of the continuous problem. First, we give the notation.

### 2.1 Notation

Given ( or ), we denote the space of real matrices of order by , and define as the space of real symmetric matrices. For matrices and , we write as usual

 τ––t=(τji),tr(τ––)=n∑i=1τii,τ––d=τ––−1ntr(τ––)I–,τ––:ζ–=n∑i,j=1τijζij, (2.1)

where

is the identity matrix.

For a subdomain and integer , we denote the scalar-valued Sobolev spaces by with the norm and seminorm . When , coincides with the Lebesgue spaces , which is equipped with the usual -inner product and -norm . The -inner product (or duality pairing) on is denoted by

. We denote the vector-valued spaces, tensor-valued function spaces and symmetric-tensor-valued spaces whose entries are in

by , and , respectively. In particular, , and . Then, we introduce the following space

 H––(div,D;M)={τ––∈L––2(D;M):divτ––∈L2(D)},

equipped with the norm for all . Here, differential operators are applied row by row, i.e., the -th row of is the divergence of the -th row vector of the matrix . Similarly, the -th row of the matrix in the definition of is the gradient (written as a row) of the -th component of the vector . We also define .

In the present context, Green’s formula takes the form

 (ε–(v),τ––)D=−(divτ––,v)D+⟨τ––nD,v⟩∂D, (2.2)

where is the exterior unit normal to . If is chosen as , we abbreviate it by using and , and similar rule follows for the spaces and norms mentioned above.

### 2.2 Brinkman model

Let be a bounded and simply connected polygonal domain with Lipschitz boundary . In this paper, we consider the permeability of the form , with the purpose of facilitating parameter-robust stability analysis. We could also take by a non-dimensionalization procedure (see Remark 2.1 below). With these simplifications, we find that for the unique solution (, ) of the Brinkman model (1.1), (, , ) solves the equations

 σ––=2ε–(u)−pI– in Ω, (2.3a) κ−1u−divσ––=f in Ω, (2.3b) divu=0 in Ω, (2.3c) u=g on Γ, (2.3d) ∫Ωpdx=0. (2.3e)

Additionally, due to the incompressibility condition, we assume that satisfies the compatibility condition , where stands for the unit outward normal on .

###### Remark 2.1.

If , by taking , and , we could eliminate the parameter , i.e.

 ˜σ––=2ε–(u)−˜pI– in Ω, κ−1u−div˜σ––=˜f in Ω, divu=0 in Ω, u=g on Γ, ∫Ω˜pdx=0.

As a result, we can get the same conclusions with the problem (2.3).

As described in Section 1, in order to keep the strengths and improve the weaknesses of the stress-velocity-pressure formulation, by the incompressible condition, the problem (2.3) can be rewritten equivalently as the pseudostress-velocity formulation (cf. Gatica and Sequeira (2018)):

 σ––d=2ε–(u) in Ω, (2.4a) κ−1u−divσ––=f in Ω, (2.4b) u=g on Γ, (2.4c) ∫Ωtr(σ––)dx=0, (2.4d)

where the pressure can be obtained by the postprocessing formula

 p=−1ntr(σ––)in Ω. (2.5)

There are two reasons for eliminating the pressure. An obvious one is to reduce one variable and, hence, many degrees of freedom in the discrete system. A more important reason is that we can use economic and accurate stable elements and develop fast solvers for the resulting discrete system so that computational cost will be greatly reduced.

Set and . Then, the variational formulation of (2.4) reads as follows: given and , find such that

 a(σ––,τ––)+b(τ––,u) =⟨τ––n,g⟩Γ ∀τ––∈Σ––, (2.6a) b(σ––,v)−s(u,v) =−(f,v) ∀v∈V. (2.6b)

Here, the bilinear forms are defined by

 a(σ––,τ––)=12(σ––d,τ––d),b(τ––,v)=(divτ––,v)ands(u,v)=(κ−1u,v).

Notice that, by Green’s formula, the equation (2.6a) contains both the equation in and the boundary condition , with in particular the incompressibility condition .

Furthermore, by the definition of , it is easy to check that

 ∥τ––∥20=∥τ––d∥20+1n∥tr(τ––)∥20, (2.7) ∥tr(τ––)∥0≤√n∥τ∥0. (2.8)

Throughout the paper, we use the abbreviation () for the inequality (), where the letter denotes a positive constant independent of the parameters , , and the mesh size , and may stand for different values at its different occurrences.

### 2.3 Well-posedness of the continuous problem

Due to the large variation of the permeability tensor, in order to show that our analysis is independent of the parameters , and are endowed with the norm

 ∥τ––∥2Σ–– =(τ––d,τ––d)+˜κ(divτ––,divτ––) ∀τ––∈Σ––, (2.9a) ∥v∥2V =˜κ−1(v,v) ∀v∈V, (2.9b)

where . We note that is indeed a norm due to the fact that for all (cf. Brezzi and Fortin (1991); Boffi et al. (2013)).

We introduce a new bilinear form on , i.e.

 A((σ––,u),(τ––,v))=a(σ––,τ––)+b(τ––,u)−b(σ––,v)+s(u,v). (2.10)

Then, the problem (2.6) can be transformed into the following problem: Find , such that

 A((σ––,u),(τ––,v))=F((τ––,v)), (2.11)

where . By Cauchy-Schwarz inequality, we have the boundedness of .

###### Lemma 2.2.

The bilinear form satisfies

 A((σ––,u),(τ––,v))≲(∥σ––∥Σ––+∥u∥V)(∥τ––∥Σ––+∥v∥V) ∀(σ––,u),(τ––,v)∈Σ––×V. (2.12)

Next, we show the inf-sup condition of at continuous level.

###### Lemma 2.3.

For any , we have

 inf(σ––,u)∈Σ––×Vsup(τ––,v)∈Σ––×VA((σ––,u),(τ––,v))(∥σ––∥Σ––+∥u∥V)(∥τ––∥Σ––+∥v∥V)≳1. (2.13)
###### Proof..

For any , there exists and a positive constant such that (cf. Arnold et al. (1984b); Brezzi and Fortin (1991))

 divσ––∗=u  in Ω,and∥σ––∗∥div≤C0∥u∥0.

Set and . Then, it is straightforward to show that

 ˜σ––∈Σ––,div˜σ––=u  in Ω,and∥˜σ––∥div≤C0∥u∥0. (2.14)

We take and , where the non-negative coefficients , and will be specified later. Then, according to Cauchy-Schwarz inequality and (2.14), one gets

 A((σ––,u),(τ––,v)) =12(σ––d,σ––d+α˜σ––d)+(div(σ––+α˜σ––),u)−(divσ––,δ1u−δ2divσ––)+(κ−1u,δ1u−δ2div–σ) =12∥σ––d∥2+α2(σ––d,˜σ––)+(1−δ1−κ−1δ2)(divσ––,u)+(α+κ−1δ1)∥u∥20+δ2∥divσ––∥20 ≥(12−αε4)∥σ––d∥2+(1−δ1−κ−1δ2)(div–σ,u)+(α+κ−1δ1−αC204ε)∥u∥20+δ2∥divσ––∥20.

From above inequality, let , , and . Then, it holds that

 12−αε4=14,1−δ1−κ−1δ2=0, δ2=˜κ2≳˜κ,α+κ−1δ1−αC204ε=1C20+κ−1(1−˜κ2κ)≳˜κ−1.

Here, we use the fact that due to the definition of . Then, we obtain

 A((σ––,u),(τ––,v))≳(∥σ––∥Σ––+∥u∥V)2.

Next, taking , and in and , by (2.14) and the fact that , one finds

 ∥τ––∥2Σ–– =∥σ––d+α˜σ––d∥20+˜κ(divσ––+αdiv˜σ––,divσ––+αdiv˜σ––) ≲∥σ––∥2Σ––+∥˜σ––d∥20+˜κ(u,u) ≲(∥σ––∥Σ––+∥u∥V)2, ∥v∥2V =˜κ−1∥δ1u−δ2divσ––∥20 ≲˜κ−1(u,u)+˜κ−1(˜κdivσ––,˜κdivσ––) ≲(∥u∥V+∥σ––∥Σ––)2.

Then, we finish this proof. ∎

From Lemma 2.2 and Lemma 2.3, we get the well-posedness of the problem (2.11).

###### Theorem 2.4.

Given and , the problem (2.11) has a unique solution .

## 3 MDG method

In this section, we formulate the MDG method for the Brinkman problem in the pseudostress-velocity formulation and show that it has a unique solution.

### 3.1 Derivation of the MDG scheme

Let be a family of quasi-regular decomposition of the domain into triangles (tetrahedrons), be the diameter of the element and . We denote the union of the boundaries of all the by , is the set of all the interior edges and is the set of boundary edges. Let and be the broken gradient and divergence operators whose restrictions on each element are equal to and , respectively. In addition, given an integer , we denote by the space of polynomials defined in of total degree at most . Recall the notation for vector-valued, tensor-valued and symmetric-tensor-valued function spaces, we have , , and . Construct the discontinuous finite element spaces and by

 Σ––h={τ––h∈L––2(Ω;S):τ––h∈P––Sk+1(K)∀K∈Th, ∫Ωtr(τ––h)dx=0}, (3.1a) Vh={vh∈L2(Ω):vh∈Pk(K)∀K∈Th}. (3.1b)

The norm of is defined by

 ∥τ––h∥2Σ––h=∥τ––dh∥20+˜κ∥divhτ––h∥20+|τ––h|2∗, (3.2)

where , is the length of edge , and denotes the -norm on edge .

For an interior edge shared by elements and , we define the unit normal vectors and on pointing exterior to and , respectively. Similarly, we define vector-valued functions and tensor-valued functions . Then define the averages and the jumps , on by

 {v}=12(v++v−),⟦v⟧=12(v+⊗n++v−⊗n−+n+⊗v++n−⊗v−), {τ––}=12(τ––++τ––−),[τ––]=12(τ––+n++τ––−n−),

where is a matrix with as its -th element. On boundary edge , we set

 {v}=v,⟦v⟧=12(v⊗n+n⊗v), {τ––}=τ––,[τ––]=τ––n.

For any tensor-valued function and vector-valued function , a straightforward computation shows that

 ∑T∈Th∫∂Kτ––nK⋅vds=∑e∈Eih∫e[τ––]⋅{v}ds+∑e∈Eh∫e{τ––}:⟦v⟧ds. (3.3)

Let us derive the MDG scheme for problem (2.4). Multiplying (2.4a) by a test function and (2.4b) by a test function , respectively, integrating on any element and applying the Green’s formula, we obtain

 12(σ––d,τ––dh)K+(divτ––h,u)K−⟨–τhnK,u⟩∂K=0 ∀τ––h∈Σ––h, (3.4a) (κ−1u,vh)K+(σ––,ε–h(vh))K−⟨σ––nK,vh⟩∂K=(f,vh)K ∀vh∈Vh. (3.4b)

Here, .

Then, we approximate and by and , respectively, and the trace of and on element edge by the numerical fluxes and . Summing on all , we get

 12(σ––dh,τ––dh)+(divhτ––h,uh)−⟨τ––hnK,ˆuh⟩∂Th=0 ∀τ––h∈Σ––h, (3.5a) (κ−1uh,vh)+(σ––h,ε–h(vh))−⟨ˆσ––hnK,vh⟩∂Th=(f,vh) ∀vh∈Vh. (3.5b)

By Green’s formula and (3.3), we have

 12(σ––dh,τ––dh)+(divhτ––h,uh)−∫Eih[τ––h]⋅{ˆuh}ds−∫Eh{τ––h}:⟦ˆuh⟧ds=0 ∀τ––h∈Σ––h, (3.6a) (κ−1uh,vh)−(divhσ––h,vh)+∫Eih[σ––h−ˆσ––h]⋅{vh}ds+∫Eh{σ––h−ˆσ––h}:⟦vh⟧ds=(f,vh) ∀vh∈Vh. (3.6b)

We define the numerical fluxes and by

 ˆσ––h={σ––h}andˆuh={uh}−ηehe[σ––h] on e∈Eih, (3.7) ˆσ––h=σ––handˆuh=g on e∈E∂h, (3.8)

where the penalty parameter . For simplicity, we choose in the analysis.

With such choices and symmetry properties of , the MDG method of the problem (2.4) is to find such that

 ah(σ––h,τ––h)+bh(τ––h,uh)=⟨τ––h