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

TL;DR

This paper introduces a stable mixed discontinuous Galerkin method for solving the Brinkman equations, effectively handling both Stokes and Darcy flow regimes with proven stability and optimal error estimates.

Contribution

The paper develops a parameter-robust mixed discontinuous Galerkin method with symmetric stress for the Brinkman problem, ensuring stability and accuracy across flow regimes.

Findings

The method is stable for both high and low permeability regimes.

Optimal error estimates are achieved for stress and velocity.

Numerical experiments confirm robustness and accuracy.

Abstract

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 $\underline{\bm{\mathcal{P}}}_{k+1}^{\mathbb{S}}$-$\bm{\mathcal{P}}_k$ $(k \geq 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 $\underline{\bm{H}}(\bm{{\rm div}})$-norm and velocity in $\bm{L}^2$-norm. Furthermore, the optimal $\underline{\bm{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.

Equations

$$\nu {\underline{\mathbf{\kappa }}}^{-1}\mathbf{u}-2\nu \mathbf{div}(\underline{\mathbf{\epsilon }}(\mathbf{u}))+\nabla p$$

$$\mathrm{div}\mathbf{u}$$

$${\underline{\mathbf{\tau }}}^t=({\tau }_{ji}),\mathrm{tr}(\underline{\mathbf{\tau }})=\sum _{i=1}^n{\tau }_{ii},{\underline{\mathbf{\tau }}}^d=\underline{\mathbf{\tau }}-\frac{1}{n}\mathrm{tr}(\underline{\mathbf{\tau }})\underline{\mathbf{I}},\underline{\mathbf{\tau }}:\underline{\mathbf{\zeta }}=\sum _{i,j=1}^n{\tau }_{ij}{\zeta }_{ij},$$

$$\underline{\mathbf{H}}(\mathrm{div},D;\mathbb{M})=\{\underline{\mathbf{\tau }}\in {\underline{\mathbf{L}}}^2(D;\mathbb{M}):\mathbf{div}\underline{\mathbf{\tau }}\in {\mathbf{L}}^2(D)\},$$

$$(\underline{\mathbf{\epsilon }}(\mathbf{v}),\underline{\mathbf{\tau }}{)}_D=-(\mathbf{div}\underline{\mathbf{\tau }},\mathbf{v}{)}_D+⟨\underline{\mathbf{\tau }}{\mathbf{n}}_D,\mathbf{v}{⟩}_{\partial D},$$

$$\underline{\mathbf{\sigma }}=2\underline{\mathbf{\epsilon }}(\mathbf{u})-p\underline{\mathbf{I}}$$

$${\kappa }^{-1}\mathbf{u}-\mathbf{div}\underline{\mathbf{\sigma }}=\mathbf{f}$$

$$\mathrm{div}\mathbf{u}=0$$

$$\mathbf{u}=\mathbf{g}$$

$${\int }_{\Omega }pd\mathbf{x}=0.$$

$$\underline{\tilde{\mathbf{\sigma }}}=2\underline{\mathbf{\epsilon }}(\mathbf{u})-\tilde{p}\underline{\mathbf{I}}$$

$${\kappa }^{-1}\mathbf{u}-\mathbf{div}\underline{\tilde{\mathbf{\sigma }}}=\tilde{\mathbf{f}}$$

$$\mathrm{div}\mathbf{u}=0$$

$$\mathbf{u}=\mathbf{g}$$

$${\int }_{\Omega }\tilde{p}d\mathbf{x}=0.$$

$${\underline{\mathbf{\sigma }}}^d=2\underline{\mathbf{\epsilon }}(\mathbf{u})$$

$${\kappa }^{-1}\mathbf{u}-\mathbf{div}\underline{\mathbf{\sigma }}=\mathbf{f}$$

$$\mathbf{u}=\mathbf{g}$$

$${\int }_{\Omega }\mathrm{tr}(\underline{\mathbf{\sigma }})d\mathbf{x}=0,$$

$$p=-\frac{1}{n}\mathrm{tr}(\underline{\mathbf{\sigma }})\text{in }\Omega .$$

$$a(\underline{\mathbf{\sigma }},\underline{\mathbf{\tau }})+b(\underline{\mathbf{\tau }},\mathbf{u})$$

$$b(\underline{\mathbf{\sigma }},\mathbf{v})-s(\mathbf{u},\mathbf{v})$$

$$a(\underline{\mathbf{\sigma }},\underline{\mathbf{\tau }})=\frac{1}{2}({\underline{\mathbf{\sigma }}}^d,{\underline{\mathbf{\tau }}}^d),b(\underline{\mathbf{\tau }},\mathbf{v})=(\mathbf{div}\underline{\mathbf{\tau }},\mathbf{v})\mathrm{and}s(\mathbf{u},\mathbf{v})=({\kappa }^{-1}\mathbf{u},\mathbf{v}).$$

$$\Vert \underline{\mathbf{\tau }}{\Vert }_0^2=\Vert {\underline{\mathbf{\tau }}}^d{\Vert }_0^2+\frac{1}{n}\Vert \mathrm{tr}(\underline{\mathbf{\tau }}){\Vert }_0^2,$$

$$\Vert \mathrm{tr}(\underline{\mathbf{\tau }}){\Vert }_0\le \sqrt{n}\Vert \mathbf{\tau }{\Vert }_0.$$

$$\Vert \underline{\mathbf{\tau }}{\Vert }_{\underline{\mathbf{\Sigma }}}^2$$

$$\Vert \mathbf{v}{\Vert }_{\mathbf{V}}^2$$

$$A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))=a(\underline{\mathbf{\sigma }},\underline{\mathbf{\tau }})+b(\underline{\mathbf{\tau }},\mathbf{u})-b(\underline{\mathbf{\sigma }},\mathbf{v})+s(\mathbf{u},\mathbf{v}).$$

$$A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))=F((\underline{\mathbf{\tau }},\mathbf{v})),$$

$$A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))\lesssim (\Vert \underline{\mathbf{\sigma }}{\Vert }_{\underline{\mathbf{\Sigma }}}+\Vert \mathbf{u}{\Vert }_{\mathbf{V}})(\Vert \underline{\mathbf{\tau }}{\Vert }_{\underline{\mathbf{\Sigma }}}+\Vert \mathbf{v}{\Vert }_{\mathbf{V}})\forall (\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v})\in \underline{\mathbf{\Sigma }}\times \mathbf{V}.$$

$$\underset{(\underline{\mathbf{\sigma }},\mathbf{u})\in \underline{\mathbf{\Sigma }}\times \mathbf{V}}{\inf }\underset{(\underline{\mathbf{\tau }},\mathbf{v})\in \underline{\mathbf{\Sigma }}\times \mathbf{V}}{\sup }\frac{A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))}{(\Vert \underline{\mathbf{\sigma }}{\Vert }_{\underline{\mathbf{\Sigma }}}+\Vert \mathbf{u}{\Vert }_{\mathbf{V}})(\Vert \underline{\mathbf{\tau }}{\Vert }_{\underline{\mathbf{\Sigma }}}+\Vert \mathbf{v}{\Vert }_{\mathbf{V}})}\gtrsim 1.$$

$$\mathbf{div}{\underline{\mathbf{\sigma }}}^{\ast }=\mathbf{u}\mathrm{in}\Omega ,\mathrm{and}\Vert {\underline{\mathbf{\sigma }}}^{\ast }{\Vert }_{\mathbf{div}}\le C_0\Vert \mathbf{u}{\Vert }_0.$$

$$\widetilde{\underline{\mathbf{\sigma }}}\in \underline{\mathbf{\Sigma }},\mathbf{div}\widetilde{\underline{\mathbf{\sigma }}}=\mathbf{u}\mathrm{in}\Omega ,\mathrm{and}\Vert \widetilde{\underline{\mathbf{\sigma }}}{\Vert }_{\mathbf{div}}\le C_0\Vert \mathbf{u}{\Vert }_0.$$

$$A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))$$

$$=\frac{1}{2}\Vert {\underline{\mathbf{\sigma }}}^d{\Vert }^2+\frac{\alpha }{2}({\underline{\mathbf{\sigma }}}^d,\widetilde{\underline{\mathbf{\sigma }}})+(1-{\delta }_1-{\kappa }^{-1}{\delta }_2)(\mathbf{div}\underline{\mathbf{\sigma }},\mathbf{u})+(\alpha +{\kappa }^{-1}{\delta }_1)\Vert \mathbf{u}{\Vert }_0^2+{\delta }_2\Vert \mathbf{div}\underline{\mathbf{\sigma }}{\Vert }_0^2$$

$$\ge (\frac{1}{2}-\frac{\alpha \epsilon }{4})\Vert {\underline{\mathbf{\sigma }}}^d{\Vert }^2+(1-{\delta }_1-{\kappa }^{-1}{\delta }_2)(\mathbf{div}\underline{\mathbf{\sigma }},\mathbf{u})+(\alpha +{\kappa }^{-1}{\delta }_1-\frac{\alpha C_0^2}{4\epsilon })\Vert \mathbf{u}{\Vert }_0^2+{\delta }_2\Vert \mathbf{div}\underline{\mathbf{\sigma }}{\Vert }_0^2.$$

$$\frac{1}{2}-\frac{\alpha \epsilon }{4}=\frac{1}{4},1-{\delta }_1-{\kappa }^{-1}{\delta }_2=0,$$

$${\delta }_2=\frac{\tilde{\kappa }}{2}\gtrsim \tilde{\kappa },\alpha +{\kappa }^{-1}{\delta }_1-\frac{\alpha C_0^2}{4\epsilon }=\frac{1}{C_0^2}+{\kappa }^{-1}(1-\frac{\tilde{\kappa }}{2\kappa })\gtrsim {\tilde{\kappa }}^{-1}.$$

$$A((\underline{\mathbf{\sigma }},\mathbf{u}),(\underline{\mathbf{\tau }},\mathbf{v}))\gtrsim (\Vert \underline{\mathbf{\sigma }}{\Vert }_{\underline{\mathbf{\Sigma }}}+\Vert \mathbf{u}{\Vert }_{\mathbf{V}}{)}^2.$$

Full text

**A mixed discontinuous Galerkin method with symmetric stress for

Brinkman problem based on the velocity-pseudostress formulation**

Yanxia Qian111School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, P. R. China. Email: yxqian0520@163.com,  Shuonan Wu222School of Mathematical Sciences, Peking University, Beijing 100871, P. R. China. The work of this author is partially supported by the startup grant from Peking University. Email: snwu@math.pku.edu.cn,  Fei Wang333School of Mathematics and Statistics & State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, P. R. China. The work of this author was partially supported by the National Natural Science Foundation of China (Grant No. 11771350). Email: feiwang.xjtu@xjtu.edu.cn (Corresponding author)

Abstract. 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 $\underline{\bm{\mathcal{P}}}_{k+1}^{\mathbb{S}}$-$\bm{\mathcal{P}}_{k}$ $(k\geq 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 $\underline{\bm{H}}(\bm{{\rm div}})$-norm and velocity in $\bm{L}^{2}$-norm. Furthermore, the optimal $\underline{\bm{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.

Keywords. Brinkman model, mixed discontinuous Galerkin method, pseudostress, parameter-robust stability

1 Introduction

The Brinkman equations (cf. [18]),

[TABLE]

which model the flow of a fluid through a complex porous medium occupying domain $\Omega$ with a high-contrast permeability tensor $\underline{\bm{\kappa}}$, can be seen as a mixture of Darcy and Stokes equations. Here, $\nu>0$ is the fluid viscosity, $\bm{u}$ denotes the velocity field, $\underline{\bm{\varepsilon}}(\bm{u})=(\underline{\bm{\nabla}}\bm{u}+(\underline{\bm{\nabla}}\bm{u})^{t})/2$ is the strain rate, $p$ is the pressure field, and $\bm{f}$ is the volume force. This model arises from applications in many fields, such as groundwater hydrology, biomedical engineering, petroleum industry, and environmental science (cf. [51, 41, 47]). From (1.1), we can see that the Brinkman problem becomes Stokes-dominated when the permeability tensor $\underline{\bm{\kappa}}$ is getting large, and it becomes Darcy-dominated when $\underline{\bm{\kappa}}$ 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 [42], the Darcy stable finite element pairs, for example, the Raviart-Thomas element (cf. [44]) 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. [3]), $\bm{\mathcal{P}_{2}}$-$\mathcal{P}_{0}$ element, nonconforming Crouzeix-Raviart (CR) finite element (cf. [24]), 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 [19], 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 [37]; stabilized equal-order finite elements are proposed and analyzed in [11]; (hybridized) interior penalty DG scheme with $\bm{H}(\rm div)$-conforming finite elements is investigated in [39]. Another approach is to develop new numerical schemes for solving Brinkman equations, for examples, pseudostress-based mixed finite element methods (cf. [5, 26]), weak Galerkin methods (cf. [43, 53]), virtual element method (cf. [20]), hybrid high-order method (cf. [10]) and hybridizable discontinuous Galerkin method (cf. [29]).

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. [15, 9]). However, the study of numerical methods for the stress-based and pseudostress-based formulations (cf. [25, 7, 21, 27]) 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. [8, 22]), Darcy flow (cf. [16, 1]), Navier-Stokes (or Stokes) equations (cf. [6, 38]), variational inequalities (cf. [49, 12, 48]) and much more. Besides, DG methods also enjoy the following advantages: (i) locally (and globally) conservative; (ii) easy to implement $hp$ adaptivity; (iii) suitable for parallel computing. We refer to [23, 17, 2, 32] for more discussion about DG methods.

In this paper, we construct a mixed discontinuous Galerkin (MDG) method with $\underline{\bm{\mathcal{P}}}_{k+1}^{\mathbb{S}}$-$\bm{\mathcal{P}}_{k}$ 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 $k\geq 0$, we get the optimal convergence order for the stress in broken $\underline{\bm{H}}(\bm{{\rm div}})$-norm and velocity in $\bm{L}^{2}$-norm; (iv) When $k\geq n$ and the Stokes pair $\bm{\mathcal{P}}_{k+2}$-$\mathcal{P}_{k+1}$ is stable, we obtain the optimal $\underline{\bm{L}}^{2}$ 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 $n$ ($n=2$ or $3$), we denote the space of real matrices of order $n\times n$ by $\mathbb{M}$, and define $\mathbb{S}\subset\mathbb{M}$ as the space of real symmetric matrices. For matrices $\underline{\bm{\tau}}=(\tau_{ij})\in\mathbb{M}$ and $\underline{\bm{\zeta}}=(\zeta_{ij})\in\mathbb{M}$, we write as usual

[TABLE]

where $\underline{\bm{I}}$ is the identity matrix.

For a subdomain $D\subset\mathbb{R}^{n}$ and integer $m\geq 0$, we denote the scalar-valued Sobolev spaces by $H^{m}(D)=W^{m,2}(D)$ with the norm $\|\cdot\|_{m,D}$ and seminorm $|\cdot|_{m,D}$. When $m=0$, $H^{0}(D)$ coincides with the Lebesgue spaces $L^{2}(D)$, which is equipped with the usual $L^{2}$-inner product $(\cdot,\cdot)_{D}$ and $L^{2}$-norm $||\cdot||_{0,D}$. The $L^{2}$-inner product (or duality pairing) on $\partial D$ is denoted by $\langle\cdot,\cdot\rangle_{\partial D}$. We denote the vector-valued spaces, tensor-valued function spaces and symmetric-tensor-valued spaces whose entries are in $H^{m}(D)$ by $\bm{H}^{m}(D)$, $\underline{\bm{H}}^{m}(D)$ and $\underline{\bm{H}}^{m}(D;\mathbb{S})$, respectively. In particular, $\bm{H}^{0}(D)=\bm{L}^{2}(D)$, $\underline{\bm{H}}^{0}(D)=\underline{\bm{L}}^{2}(D)$ and $\underline{\bm{H}}^{0}(D;\mathbb{S})=\underline{\bm{L}}^{2}(D;\mathbb{S})$. Then, we introduce the following space

[TABLE]

equipped with the norm $\|\underline{\bm{\tau}}\|_{\bm{{\rm div}},D}=(\|\underline{\bm{\tau}}\|^{2}_{0,D}+\|\mathbf{div}\underline{\bm{\tau}}\|^{2}_{0,D})^{1/2}$ for all $\underline{\bm{\tau}}\in\,\underline{\bm{H}}(\bm{{\rm div}};\mathbb{M})$. Here, differential operators are applied row by row, i.e., the $i$-th row of $\bm{{\rm div}}\underline{\bm{\sigma}}$ is the divergence of the $i$-th row vector of the matrix $\underline{\bm{\sigma}}$. Similarly, the $i$-th row of the matrix $\underline{\bm{\nabla}}\bm{u}$ in the definition of $\underline{\bm{\varepsilon}}(\bm{u})$ is the gradient (written as a row) of the $i$-th component of the vector $\bm{u}$. We also define $\underline{\bm{H}}(\bm{{\rm div}},D;\mathbb{S})=\{\underline{\bm{\tau}}\in\underline{\bm{L}}^{2}(D;\mathbb{S}):\mathbf{div}\underline{\bm{\tau}}\in\bm{L}^{2}(D)\}$.

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

[TABLE]

where $\bm{n}_{D}$ is the exterior unit normal to $\partial D$. If $D$ is chosen as $\Omega$, we abbreviate it by using $(\cdot,\cdot)$ and $\langle\cdot,\cdot\rangle$, and similar rule follows for the spaces and norms mentioned above.

2.2 Brinkman model

Let $\Omega\subset\mathbb{R}^{n}$ be a bounded and simply connected polygonal domain with Lipschitz boundary $\Gamma$. In this paper, we consider the permeability of the form $\underline{\bm{\kappa}}=\kappa\underline{\bm{I}}$, with the purpose of facilitating parameter-robust stability analysis. We could also take $\nu=1$ by a non-dimensionalization procedure (see Remark 2.1 below). With these simplifications, we find that for the unique solution ($\bm{u}$, $p$) of the Brinkman model (1.1), ($\underline{\bm{\sigma}}$, $\bm{u}$, $p$) solves the equations

[TABLE]

Additionally, due to the incompressibility condition, we assume that $\bm{g}$ satisfies the compatibility condition $\int_{\Gamma}\bm{g}\cdot\bm{n}ds=0$, where $\bm{n}$ stands for the unit outward normal on $\Gamma$.

Remark 2.1**.**

If $\nu\neq 1$, by taking $\underline{\widetilde{\bm{\sigma}}}=\underline{\bm{\sigma}}/\nu$, $\widetilde{\bm{f}}=\bm{u}/{\nu}$ and $\widetilde{p}=p/{\nu}$, we could eliminate the parameter $\nu$, i.e.

[TABLE]

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. [29]):

[TABLE]

where the pressure $p$ can be obtained by the postprocessing formula

[TABLE]

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 $\underline{\bm{\Sigma}}=\{\underline{\bm{\tau}}\in\underline{\bm{H}}(\bm{{\rm div}};\mathbb{S}):\int_{\Omega}{\rm tr}(\underline{\bm{\tau}})d\bm{x}=0\}$ and $\bm{V}=\bm{L}^{2}(\Omega)$. Then, the variational formulation of (2.4) reads as follows: given $\bm{f}\in\bm{L}^{2}(\Omega)$ and $\bm{g}\in\bm{H}^{1/2}(\Gamma)$, find $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{V}$ such that

[TABLE]

Here, the bilinear forms are defined by

[TABLE]

Notice that, by Green’s formula, the equation (2.6a) contains both the equation $2\underline{\bm{\varepsilon}}(\bm{u})=\underline{\bm{\sigma}}^{d}$ in $\Omega$ and the boundary condition $\bm{u}=\bm{g}$, with in particular the incompressibility condition ${\rm div}\,\bm{u}=\frac{1}{2}{\rm tr}(\underline{\bm{\sigma}}^{d})=0$.

Furthermore, by the definition of $\underline{\bm{\tau}}^{d}$, it is easy to check that

[TABLE]

Throughout the paper, we use the abbreviation $x\lesssim y$ ($x\gtrsim y$) for the inequality $x\leq Cy$ ($x\geq Cy$), where the letter $C$ denotes a positive constant independent of the parameters $\kappa$, $\nu$, and the mesh size $h$, 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 $\kappa$, $\underline{\bm{\Sigma}}$ and $\bm{V}$ are endowed with the norm

[TABLE]

where $\widetilde{\kappa}=\min\{\kappa,1\}$. We note that $\|\cdot\|_{\underline{\bm{\Sigma}}}$ is indeed a norm due to the fact that $\|\underline{\bm{\tau}}\|_{0}^{2}\lesssim\|\underline{\bm{\tau}}^{d}\|_{0}^{2}+\|\mathbf{div}\underline{\bm{\tau}}\|_{0}^{2}$ for all $\underline{\bm{\tau}}\in\underline{\bm{\Sigma}}$ (cf. [15, 9]).

We introduce a new bilinear form $A((\underline{\bm{\sigma}},\bm{u}),(\underline{\bm{\tau}},\bm{v}))$ on $(\underline{\bm{\Sigma}}\times\bm{V},\underline{\bm{\Sigma}}\times\bm{V})$, i.e.

[TABLE]

Then, the problem (2.6) can be transformed into the following problem: Find $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{V}$, such that

[TABLE]

where $F((\underline{\bm{\tau}},\bm{v}))=(\bm{f},\bm{v})+\langle\underline{\bm{\tau}}\bm{n},\bm{g}\rangle_{\Gamma}$. By Cauchy-Schwarz inequality, we have the boundedness of $A(\cdot,\cdot)$.

Lemma 2.2**.**

The bilinear form $A(\cdot,\cdot)$ satisfies

[TABLE]

Next, we show the inf-sup condition of $A(\cdot,\cdot)$ at continuous level.

Lemma 2.3**.**

For any $(\underline{\bm{\sigma}},\bm{u}),(\underline{\bm{\tau}},\bm{v})\in\underline{\bm{\Sigma}}\times\bm{V}$, we have

[TABLE]

Proof..

For any $\bm{u}\in\bm{V}$, there exists ${\underline{\bm{\sigma}}}^{*}\in\underline{\bm{H}}(\bm{{\rm div}};\mathbb{S})$ and a positive constant $C_{0}>0$ such that (cf. [4, 15])

[TABLE]

Set $\gamma=\frac{1}{|\Omega|}\int_{\Omega}\mathrm{tr}({\underline{\bm{\sigma}}}^{*})d\bm{x}$ and $\widetilde{\underline{\bm{\sigma}}}={\underline{\bm{\sigma}}}^{*}-\frac{\gamma}{n}\underline{\bm{I}}$. Then, it is straightforward to show that

[TABLE]

We take $\underline{\bm{\tau}}=\underline{\bm{\sigma}}+\alpha\widetilde{\underline{\bm{\sigma}}}$ and $\bm{v}=\delta_{1}\bm{u}-\delta_{2}\mathbf{div}\underline{\bm{\sigma}}$, where the non-negative coefficients $\alpha$, $\delta_{1}$ and $\delta_{2}$ will be specified later. Then, according to Cauchy-Schwarz inequality and (2.14), one gets

[TABLE]

From above inequality, let $\delta_{1}=1-\frac{\widetilde{\kappa}}{2\kappa}$, $\delta_{2}=\frac{\widetilde{\kappa}}{2}$, $\varepsilon=\frac{C_{0}^{2}}{2}$ and $\alpha=\frac{2}{C_{0}^{2}}$. Then, it holds that

[TABLE]

Here, we use the fact that $\frac{1}{2}\leq\delta_{1}=1-\frac{\widetilde{\kappa}}{2\kappa}<1$ due to the definition of $\widetilde{\kappa}$. Then, we obtain

[TABLE]

Next, taking $\alpha$, $\delta_{1}$ and $\delta_{2}$ in $\underline{\bm{\tau}}$ and $\bm{v}$, by (2.14) and the fact that $\widetilde{\kappa}\leq 1\leq\widetilde{\kappa}^{-1}$, one finds

[TABLE]

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 $\bm{f}\in\bm{L}^{2}(\Omega)$ and $\bm{g}\in\bm{H}^{1/2}(\Gamma)$, the problem (2.11) has a unique solution $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{V}$.

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 $\{\mathcal{T}_{h}\}_{h}$ be a family of quasi-regular decomposition of the domain $\overline{\Omega}$ into triangles (tetrahedrons), $h_{K}$ be the diameter of the element $K\in\mathcal{T}_{h}$ and $h=\max\{h_{K}:K\in\mathcal{T}_{h}\}$. We denote the union of the boundaries of all the $K\in\mathcal{T}_{h}$ by $\mathcal{E}_{h}$, $\mathcal{E}^{i}_{h}$ is the set of all the interior edges and $\mathcal{E}_{h}^{\partial}=\mathcal{E}_{h}/\mathcal{E}^{i}_{h}$ is the set of boundary edges. Let $\underline{\bm{\nabla}}_{h}$ and $\mathbf{div}_{h}$ be the broken gradient and divergence operators whose restrictions on each element $K\in\mathcal{T}_{h}$ are equal to $\underline{\bm{\nabla}}$ and $\mathbf{div}$, respectively. In addition, given an integer $k\geq 0$, we denote by $\mathcal{P}_{k}(D)$ the space of polynomials defined in $D$ of total degree at most $k$. Recall the notation for vector-valued, tensor-valued and symmetric-tensor-valued function spaces, we have $\bm{\mathcal{P}}_{k}(D)=[\mathcal{P}_{k}(D)]^{n}$, $\underline{\bm{\mathcal{P}}}^{\mathbb{M}}_{k}(D)=[\mathcal{P}_{k}(D)]^{n\times n}$, and $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{k}(D)=\{\underline{\bm{\tau}}\in[\mathcal{P}_{k}(D)]^{n\times n}:\underline{\bm{\tau}}^{t}=\underline{\bm{\tau}}\}$. Construct the discontinuous finite element spaces $\underline{\bm{\Sigma}}_{h}$ and $\bm{V}_{h}$ by

[TABLE]

The norm of $\underline{\bm{\Sigma}}_{h}$ is defined by

[TABLE]

where $|\underline{\bm{\tau}}_{h}|_{*}^{2}=\sum_{e\in\mathcal{E}^{i}_{h}}h_{e}^{-1}\|[\underline{\bm{\tau}}_{h}]\|_{e}^{2}$, $h_{e}$ is the length of edge $e$, and $\|\cdot\|_{e}$ denotes the $L^{2}$-norm on edge $e$.

For an interior edge $e\in\mathcal{E}^{i}_{h}$ shared by elements $K^{+}$ and $K^{-}$, we define the unit normal vectors $\bm{n}^{+}$ and $\bm{n}^{-}$ on $e$ pointing exterior to $K^{+}$ and $K^{-}$, respectively. Similarly, we define vector-valued functions $\bm{v}^{\pm}=\bm{v}|_{\partial K^{\pm}}$ and tensor-valued functions $\underline{\bm{\tau}}^{\pm}=\underline{\bm{\tau}}|_{\partial K^{\pm}}$. Then define the averages $\{\cdot\}$ and the jumps $\llbracket\cdot\rrbracket$, $[\cdot]$ on $e\in\mathcal{E}^{i}_{h}$ by

[TABLE]

where $\bm{v}\otimes\bm{w}$ is a matrix with $v_{i}w_{j}$ as its $(i,j)$-th element. On boundary edge $e\in\mathcal{E}_{h}^{\partial}$, we set

[TABLE]

For any tensor-valued function $\underline{\bm{\tau}}$ and vector-valued function $\bm{v}$, a straightforward computation shows that

[TABLE]

Let us derive the MDG scheme for problem (2.4). Multiplying (2.4a) by a test function $\underline{\bm{\tau}}_{h}$ and (2.4b) by a test function $\bm{v}_{h}$, respectively, integrating on any element $K\in\mathcal{T}_{h}$ and applying the Green’s formula, we obtain

[TABLE]

Here, $\underline{\bm{\varepsilon}}_{h}(\bm{v}_{h})=(\underline{\bm{\nabla}}_{h}\bm{v}_{h}+(\underline{\bm{\nabla}}_{h}\bm{v}_{h})^{t})/2$.

Then, we approximate $\underline{\bm{\sigma}}$ and $\bm{u}$ by $\underline{\bm{\sigma}}_{h}\in\underline{\bm{\Sigma}}_{h}$ and $\bm{u}_{h}\in\bm{V}_{h}$, respectively, and the trace of $\underline{\bm{\sigma}}$ and $\bm{u}$ on element edge by the numerical fluxes $\widehat{\underline{\bm{\sigma}}}_{h}$ and $\widehat{\bm{u}}_{h}$. Summing on all $K\in\mathcal{T}_{h}$, we get

[TABLE]

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

[TABLE]

We define the numerical fluxes $\widehat{\underline{\bm{\sigma}}}_{h}$ and $\widehat{\bm{u}}_{h}$ by

[TABLE]

where the penalty parameter $\eta_{e}=\mathcal{O}(1)$. For simplicity, we choose $\eta_{e}=1$ in the analysis.

With such choices and symmetry properties of $\underline{\bm{\tau}}_{h}$, the MDG method of the problem (2.4) is to find $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})\in\underline{\bm{\Sigma}}_{h}\times\bm{V}_{h}$ such that

[TABLE]

where

[TABLE]

3.2 Well-posedness of the MDG method

In this subsection, we show the well-posedness of the MDG scheme (3.9). First, we give some inequalities by lemmas.

The first lemma is a discrete analogy of [9, Proposition 9.1.1], which indicates the well-posedness of the discrete norm (3.2).

Lemma 3.1** (Lemma 3.3 in [28]).**

For every $\underline{\bm{\tau}}_{h}\in\underline{\bm{\Sigma}}_{h}$, it holds

[TABLE]

Lemma 3.2** (cf. [2]).**

There exists positive constants $C_{1}$ and $C_{2}$ such that

[TABLE]

Then, we define

[TABLE]

and

[TABLE]

Equivalently, the problem (3.9) can be rewritten as the following problem: Find $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})\in\underline{\bm{\Sigma}}_{h}\times\bm{V}_{h}$, such that

[TABLE]

In what follows, we prove the well-posedness of problem (3.13).

Lemma 3.3**.**

The bilinear form $A_{h}(\cdot,\cdot)$ satisfies

[TABLE]

Proof..

By the definition of norm (3.2) and Cauchy-Schwarz inequality, we have

[TABLE]

where we use the trace and inverse inequalities and the fact that $\widetilde{\kappa}\leq 1$. ∎

In order to prove the discrete inf-sup condition, we introduce the space (cf. [52])

[TABLE]

For any $\bm{u}_{h}\in\bm{V}_{h}$, there exists a constant $C_{3}>0$ and $\underline{\bm{\sigma}}_{h}^{*}\in\underline{\bm{\Sigma}}_{h}^{\rm NC}$ (cf. [52]) such that

[TABLE]

We are now in the position to show the inf-sup condition of $A_{h}(\cdot,\cdot)$.

Lemma 3.4**.**

For any $(\underline{\bm{\sigma}}_{h},\bm{u}_{h}),(\underline{\bm{\tau}}_{h},\bm{v}_{h})\in\underline{\bm{\Sigma}}_{h}\times\bm{V}_{h}$, it holds

[TABLE]

Proof..

Set $\gamma=\frac{1}{|\Omega|}\int_{\Omega}\mathrm{tr}(\underline{\bm{\sigma}}_{h}^{*})d\bm{x}$ and $\widetilde{\underline{\bm{\sigma}}}_{h}=\underline{\bm{\sigma}}_{h}^{*}-\frac{\gamma}{n}\underline{\bm{I}}$. Then, it is straightforward to show that

[TABLE]

We take $\underline{\bm{\tau}}_{h}=\underline{\bm{\sigma}}_{h}+\alpha\widetilde{\underline{\bm{\sigma}}}_{h}$ and $\bm{v}_{h}=\delta_{1}\bm{u}_{h}-\delta_{2}\mathbf{div}_{h}\underline{\bm{\sigma}}_{h}$, where the non-negative undetermined coefficients $\alpha$, $\delta_{1}$ and $\delta_{2}$ will be specified in the following analysis. Thanks to (3.17) and $\int_{\mathcal{E}_{h}^{i}}[\widetilde{\underline{\bm{\sigma}}}_{h}]\cdot\{\bm{u}_{h}\}ds=0$ obtained from the property of $\underline{\bm{\Sigma}}_{h}^{\rm NC}$, one gets

[TABLE]

According to Cauchy-Schwarz inequality, (3.17), the trace and inverse inequalities, we obtain

[TABLE]

Combining with above inequalities, we have

[TABLE]

Now, we need to choose appropriate parameters $\delta_{1}$, $\delta_{2}$, $\varepsilon$ and $\alpha$, such that

[TABLE]

We could take $\delta_{2}=\frac{\widetilde{\kappa}}{\max\{2,3C_{1}(1+C_{2}),2C_{3}\}}$, $\delta_{1}=1-\frac{\delta_{2}}{\kappa}$, $\varepsilon=1$ and $\alpha=\frac{1}{2C_{3}}$, by which the first three requirements above meet easily. Furthermore, we have

[TABLE]

Using $\frac{1}{2}\leq\delta_{1}<1$ and $\widetilde{\kappa}^{-1}=\max\{1,\kappa^{-1}\}$, we show the last inequality into two cases:

[TABLE]

From the above, we obtain

[TABLE]

Next, taking $\alpha$, $\delta_{1}$ and $\delta_{2}$ in $\underline{\bm{\tau}}_{h}$ and $\bm{v}_{h}$, due to (3.17), the fact that $\delta_{1}<1$ and $\widetilde{\kappa}\leq 1\leq\widetilde{\kappa}^{-1}$ , it holds that

[TABLE]

Then, we finish this proof. ∎

From Lemma 3.3 and Lemma 3.4, the well-posedness of the problem (3.13) can be obtained.

Theorem 3.5**.**

The mixed DG scheme (3.13) has a unique solution $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})\in\underline{\bm{\Sigma}}_{h}\times\bm{V}_{h}$.

Remark 3.6**.**

For the bilinear form $A_{h}(\bm{u}_{h},\bm{v}_{h})$ of (3.13), when the penalty term $\int_{\mathcal{E}_{h}^{i}}\frac{1}{h_{e}}[\underline{\bm{\sigma}}_{h}]\cdot[\underline{\bm{\tau}}_{h}]ds$ is replaced by $\sum_{e\in\mathcal{E}_{h}^{i}}\int_{\Omega}r_{e}([\underline{\bm{\sigma}}_{h}])\cdot r_{e}([\underline{\bm{\tau}}_{h}])d\bm{x}$, we can obtain another mixed DG scheme, which is the dual form of the method of Brezzi et al. [6]. Here, $r_{e}:(L^{2}({\cal E}_{h}))^{d}\rightarrow\bm{V}_{h}$ is the lifting operator (cf. [6, 2, 50]). The well-posedness of the corresponding scheme can be proved similarly.

4 Error estimates

In this section, we aim to derive the error estimates for the MDG scheme (3.9). First, we show the consistency of the MDG scheme, which naturally leads to an error estimate by the inf-sup condition.

4.1 Error estimate in energy norm for the pseudostress and

velocity

Lemma 4.1**.**

Let the solution $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{H}^{1}(\Omega)$, then

[TABLE]

Proof..

Since $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{H}^{1}(\Omega)$, we have $[\underline{\bm{\sigma}}]=0$ and $\llbracket\bm{u}\rrbracket=0$ on each edge $e\in\mathcal{E}_{h}^{i}$, then

[TABLE]

and

[TABLE]

Then, from (3.9), we have

[TABLE]

By (3.12), we complete the proof. ∎

Theorem 4.2**.**

Let $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{H}^{1}(\Omega)$ be the solution of (2.11), and $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})\in\underline{\bm{\Sigma}}_{h}\times\bm{V}_{h}$ be the solution of (3.9). Then, we have

[TABLE]

Proof..

From Lemma 2.3, Lemma 3.4 and Lemma 4.1, we obtain that for any $\underline{\bm{\tau}}_{h}\in\underline{\bm{\Sigma}}_{h}$ and $\bm{v}_{h}\in\bm{V}_{h}$,

[TABLE]

Then, the triangle inequality indicates the estimate (4.4). ∎

Recall that $\widetilde{\kappa}^{-1}\geq 1$ and the definition of $\|\cdot\|_{\bm{V}}$ in (2.9b), the above theorem shows the parameter-robust error estimate of $\bm{u}_{h}$ in $\bm{L}^{2}$ norm by the standard interpolation theory (cf. [46]).

Theorem 4.3**.**

Assume that the solution of (2.11) satisfies $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{H}}^{k+2}(\Omega)\times\bm{H}^{k+1}(\Omega)$. Then, the solution of the mixed DG problem (3.9) satisfies for any $\kappa>0$,

[TABLE]

Further, if $\kappa\gtrsim 1$ (high permeability case), we have

[TABLE]

Here, the hidden constants in (4.6) and (4.7) are both independent of $\kappa$.

Remark 4.4**.**

By Lemma 3.1, we have that if $\kappa\gtrsim 1$, then $\|\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h}\|_{0}\lesssim h^{k+1}(|\underline{\bm{\sigma}}|_{k+2}+|\bm{u}|_{k+1})$.

4.2 Parameter-robust error estimate of pseudostress

In this subsection, we show a parameter-robust error estimate of pseudostress for arbitrary permeability, which fills the gap of (4.7) in Theorem 4.3. The result hinges on the parameter-robust estimate of velocity given in (4.6).

Let $\bm{P}_{h}:\bm{V}\rightarrow\bm{V}_{h}$ denote the $\bm{L}^{2}$-orthogonal projection defined by

[TABLE]

For $\underline{\bm{\sigma}}\in\underline{\bm{H}}^{k+2}(\Omega)$, let $\underline{\bm{\sigma}}_{I}^{\rm SZ}\in\underline{\bm{H}}^{1}(\Omega)$ be the Scott-Zhang interpolation (cf. [46]) that satisfies

[TABLE]

Next, we modify the Scott-Zhang interpolation by ${\underline{\bm{\sigma}}}_{I}=\underline{\bm{\sigma}}_{I}^{\rm SZ}-\frac{1}{n|\Omega|}\int_{\Omega}\mathrm{tr}(\underline{\bm{\sigma}}_{I}^{\rm SZ})\underline{\bm{I}}d\bm{x}\in\underline{\bm{H}}^{1}(\Omega)\cap\underline{\bm{\Sigma}}_{h}$. Using the property of Scott-Zhang interpolation in (4.9) and the fact that $\int_{\Omega}\mathrm{tr}(\underline{\bm{\sigma}})d\bm{x}=0$, we have

[TABLE]

Hence, the modified Scott-Zhang interpolation has the same approximation as the standard one, i.e.,

[TABLE]

We also denote $\underline{\bm{e}}_{\underline{\bm{\sigma}}}=\underline{\bm{\sigma}}_{I}-\underline{\bm{\sigma}}_{h}$, ${\bm{e}}_{\bm{u}}=\bm{P}_{h}\bm{u}-\bm{u}_{h}$.

Theorem 4.5**.**

Assume that the solution of (2.11) satisfies $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{H}}^{k+2}(\Omega)\times\bm{H}^{k+1}(\Omega)$. Then, the solution of the mixed DG problem (3.9) satisfies for any $\kappa>0$,

[TABLE]

where the hidden constant is independent of $\kappa$.

Proof..

Taking $\underline{\bm{\tau}}_{h}=\underline{\bm{e}}_{\underline{\bm{\sigma}}}$ and $\bm{v}_{h}=\bm{e}_{\bm{u}}$ in the error equation (4.3), we have

[TABLE]

Subtracting the above equations, we get

[TABLE]

Using the property of Scott-Zhang interpolation (4.10), the estimate (4.6), Cauchy-Schwarz inequality, and trace inequality, we have

[TABLE]

Taking above inequalities into the right hand side of (4.12), we have

[TABLE]

which leads to the desired estimate (4.11) by triangle inequality. ∎

Remark 4.6**.**

As a byproduct in the proof of above theorem, it can be seen that, under the condition of Theorem 4.5,

[TABLE]

which implies that $\bm{u}_{h}\to\bm{P}_{h}\bm{u}$ as $\kappa\to 0$.

We then have the error estimate of pressure in the following theorem.

Theorem 4.7**.**

Assume that the solution of (2.11) satisfies $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{H}}^{k+2}(\Omega)\times\bm{H}^{k+1}(\Omega)$. Then, the solution of the mixed DG problem (3.9) satisfies

[TABLE]

Proof..

Taking $\bm{v}_{h}=\mathbf{div}_{h}{\underline{\bm{e}}}_{\underline{\bm{\sigma}}}$ in (4.3b), we obtain

[TABLE]

Using Theorem 4.5 and (4.13), one gets

[TABLE]

In light of (4.10), we apply Lemma 3.1 to have

[TABLE]

which leads to the desired estimate. ∎

4.3 Improved error estimates in $L^{2}$ norm for the

pseudostress and pressure

In this subsestion, following a similar argument in [50], we show that the $L^{2}$ error estimates for pesudostress and pressure are both optimal when the Stokes finite element pair $\bm{\mathcal{P}}^{c}_{k+2}$-$\mathcal{P}_{k+1}$ ($k\geq n$) is stable. Here, $\mathcal{P}^{c}$ represents the conforming polynomial.

Now, we recall the classical BDM projection $\pi_{h}^{c}$ (cf. [14] for two-dimension case and [13] for three-dimensional case). The projection $\pi_{h}^{c}:\bm{H}({\rm div})\to\bm{W}_{h}=\{\bm{v}\in\bm{H}({\rm div}):\bm{v}|_{K}\in\bm{\mathcal{P}}_{k+1}(K)\ \forall K\in\mathcal{T}_{h}\}$ is defined by

[TABLE]

Here, $\bm{W}_{h,*}(K)=\{\bm{z}\in\bm{\mathcal{P}}_{k+1}(K):\bm{z}\cdot\bm{n}=0\ {\rm on}\ e\in\partial K\ {\rm and}\ (\bm{z},\nabla q_{h})_{K}=0\ \forall\,q_{h}\in\mathcal{P}_{k}(K)\}$.

Then, based on the projection (4.15), on each element $K\in\mathcal{T}_{h}$, we define a function $\widetilde{\underline{\bm{\sigma}}}_{h}$ as the only element of $\underline{\bm{\mathcal{P}}}^{\mathbb{M}}_{k+1}(K)$ by $\widehat{\underline{\bm{\sigma}}}_{h}$ and $\underline{\bm{\sigma}}_{h}$ in (3.5).

[TABLE]

where $\underline{\bm{\Sigma}}_{h,*}^{c}(K)=\{\underline{\bm{\theta}}\in\underline{\bm{\mathcal{P}}}^{\mathbb{M}}_{k+1}(K):\underline{\bm{\theta}}\bm{n}=\mathbf{0}\ {\rm on}\ e\in\partial K\ {\rm and}\ (\underline{\bm{\theta}},\underline{\bm{\nabla}}\bm{v}_{h})_{K}=0\ \forall\,\bm{v}_{h}\in\bm{\mathcal{P}}_{k}(K)\}$.

The system (4.16) can be regarded as the row-wise BDM projection. According to the definition of $\pi_{h}^{c}$ and the fact that the normal component of the numerical trace for the flux is single-valued, we have the following lemma.

Lemma 4.8**.**

The function $\widetilde{\underline{\bm{\sigma}}}_{h}$ in (4.16) is well-defined,

[TABLE]

Proof..

According to definition of BDM projection, the fact that the normal component of the numerical trace for the flux is single-valued and $\int_{\Omega}{\rm tr}(\underline{\bm{\sigma}}_{h})d\bm{x}=0$, the well-posedness and (4.17) are directly available. Setting $\underline{\bm{\sigma}}=\widetilde{\underline{\bm{\sigma}}}_{h}-\underline{\bm{\sigma}}_{h}$, from (4.16), we know

[TABLE]

By the standard scaling argument, we have (4.18). ∎

Then, by the similar argument in [30, 50], we symmetrize $\widetilde{\underline{\bm{\sigma}}}_{h}$ to establish the $\underline{\bm{L}}^{2}$ error estimate for pseudostress variable. With the help of the stable Stokes pair $\bm{\mathcal{P}}^{c}_{k+2}$-$\mathcal{P}_{k+1}$ ($k\geq n$), one finds the following result. We refer the reader to [50] for detailed discussion.

Lemma 4.9** (cf. [50]).**

Assume that the Stokes pair $\bm{\mathcal{P}}^{c}_{k+2}$-$\mathcal{P}_{k+1}$ ($k\geq n$) is stable on the decomposition $\mathcal{T}_{h}$. For $\widetilde{\underline{\bm{\sigma}}}_{h}$ given in (4.16), there exists $\widetilde{\underline{\bm{\tau}}}_{h}\in\underline{\bm{\Sigma}}_{h}^{c}$ such that $\underline{\bm{\sigma}}_{h,*}=\widetilde{\underline{\bm{\sigma}}}_{h}+\widetilde{\underline{\bm{\tau}}}_{h}\in\underline{\bm{H}}(\mathbf{div};\mathbb{S})$,

[TABLE]

Next, we begin to show the optimal $\underline{\bm{L}}^{2}$ error estimate. In [35], the conforming mixed element $\underline{\bm{\mathcal{P}}}^{c}_{k+1}$-$\bm{\mathcal{P}}_{k}$ ($k\geq n$) is constructed on simplicial grids. Moreover, when $k\geq n$, there exists a projection $\Pi_{h}^{c}$ satisfying ([33])

[TABLE]

Theorem 4.10**.**

Let the solutions $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{\Sigma}}\times\bm{H}^{1}(\Omega)$ and $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})$ be the solutions of MDG problems (3.9). Under the condition of Lemma 4.9, we have

[TABLE]

Further, the following superconvergence holds

[TABLE]

Proof..

According to (3.5b), (4.16) and Lemma 4.8, for any $\bm{v}_{h}\in\bm{V}_{h}$, we have

[TABLE]

Applying Lemma 4.9, there exist $\widetilde{\underline{\bm{\tau}}}_{h}\in\underline{\bm{\Sigma}}_{h}^{c}$ such that the symmetrized variable $\underline{\bm{\sigma}}_{h,*}=\widetilde{\underline{\bm{\sigma}}}_{h}+\widetilde{\underline{\bm{\tau}}}_{h}$ is piecewise $\underline{\bm{\mathcal{P}}}^{\mathbb{M}}_{k+1}(K)$ and $\underline{\bm{\sigma}}_{h,*}\in\underline{\bm{H}}(\mathbf{div};\mathbb{S})$. Then,

[TABLE]

By (4.2), (4.8), (4.20a) and (4.23), it holds that

[TABLE]

Taking $\bm{v}_{h}=\bm{e}_{\bm{u}}$ in (4.24) and $\underline{\bm{\tau}}_{h}=\Pi_{h}^{c}\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h,*}$ in (4.3a), due to the $\underline{H}(\mathbf{div})$ conformity of $\Pi_{h}^{c}\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h,*}$, we obtain

[TABLE]

which implies that

[TABLE]

Then, it holds

[TABLE]

Combining the above inequality with Lemma 4.8, Lemma 4.9 and (4.11), we get

[TABLE]

which finishes the proof. ∎

Remark 4.11**.**

For two dimensional case, when $k\geq 2$, the conforming mixed element pair $\underline{\bm{\mathcal{P}}}^{c}_{k+1}$-$\bm{\mathcal{P}}_{k}$ is constructed in [34] and the Scott-Vogelius elements $\bm{\mathcal{P}}^{c}_{k+2}$-$\mathcal{P}_{k+1}$ are stable (cf. [45, 31]). Therefore, the optimal $\underline{\bm{L}}^{2}$ estimate for pseudostress holds when $k\geq 2$ in two dimensional case.

Theorem 4.12**.**

Assume that the solution of (2.11) satisfies $(\underline{\bm{\sigma}},\bm{u})\in\underline{\bm{H}}^{k+2}(\Omega)\times\bm{H}^{k+1}(\Omega)$. Under the condition of Theorem 4.10, the solution of the mixed DG problem (3.9) satisfies

[TABLE]

Further, we have the optimal error estimate for pressure $p$ in $L^{2}$-norm

[TABLE]

Proof..

Taking $\bm{v}_{h}=\bm{{\rm div}}_{h}(\Pi_{h}^{c}\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h,*})$ in (4.24), by the superconvergence result in (4.22), we obtain

[TABLE]

Here, we use the fact that $\kappa^{-1}\leq\widetilde{\kappa}^{-1}$. Again, by Lemma 4.9, Lemma 4.8 and Theorem 4.5, we have

[TABLE]

whence, by (4.20b), (4.21) and (4.28),

[TABLE]

Note here that $\Pi_{h}^{c}\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h,*}\in\underline{\bm{\Sigma}}_{h}\cap\underline{\bm{H}}(\mathbf{div},\Omega;\mathbb{S})$, using Lemma 3.1, (4.28), one gets

[TABLE]

which yields

[TABLE]

Note that (2.5), we can define the numerical solution of pressure by the postprocessed approximation $p_{h}=-\frac{1}{n}{\rm tr}(\underline{\bm{\sigma}}_{h})$. The optimal $L^{2}$ estimate for pressure (4.27) then follows from the fact that $\|p-p_{h}\|_{0}\lesssim\|\mathrm{tr}(\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h})\|_{0}$. ∎

5 Numerical examples

In this section, we present some numerical results to illustrate the reliability, accuracy, and flexibility of the MDG method (3.13). The numerical results presented below are obtained by using Fenics software (cf. [40]). For simplicity, we consider the triangular meshes in the two-dimensional case and the discontinuous Galerkin pair $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{k+1}$-$\bm{\mathcal{P}}_{k}$ for all the numerical examples.

Example 1 is employed to illustrate the performance of the MDG scheme (3.13) for different permeability with the polynomial degrees $k=0,1,2$. For the variable permeability, Example 2 is used to test the accuracy of the MDG scheme (3.13) with different viscosity. Example 3 and Example 4 are utilized to show the behavior of MDG scheme (3.13) for the Brinkman problem in a region with different contrast permeability.

Example 1**.**

Consider the steady Brinkman problem (2.3) in a square domain $(0,1)\times(0,1)$ with a homogeneous boundary condition that $\bm{u}=\bm{0}$ on $\Gamma$. The right hand side function $\bm{f}$ and the exact stress function $\underline{\bm{\sigma}}$ are selected such that the exact solution is given by

[TABLE]

This example aims at testing the accuracy and reliability of the MDG method for fixed viscosity and different permeability. Set $1/h=4$, $8$, $16$, $32$ and $\nu=1$. We compute the numerical solutions $(\underline{\bm{\sigma}}_{h},\bm{u}_{h})$ on uniform meshes with $\kappa^{-1}=10^{-3},10^{0}$, and $10^{3}$. The numerical results of $\|\bm{u}-\bm{u}_{h}\|_{0}$, $\|\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h}\|_{\underline{\bm{\Sigma}}_{h}}$, $\|\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h}\|_{0}$ and $\|p-p_{h}\|_{0}$ for finite element pairs $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{k+1}$-$\bm{\mathcal{P}}_{k}$ $(k=0,1{\rm\ and\ }2)$ are given in Table 1–Table 3, respectively. The numerical results confirm the optimal convergence orders, which are consistent with the theoretical results developed in Section 4. We can see that the MDG scheme is very stable with respected to different permeability.

Example 2**.**

In this example, choose $\Omega=(0,1)\times(0,1)$ for the steady Brinkman problem (2.3). The right hand side function $\bm{f}$, the exact stress function $\underline{\bm{\sigma}}$ and boundary condition $\bm{g}$ are selected such that the exact solution is given by

[TABLE]

This example aims at testing the accuracy and reliability of the MDG method for different viscosity and variable permeability. The MDG finite element pair $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{2}$-$\bm{\mathcal{P}}_{1}$ is employed in the numerical discretization on uniform meshes. The parameter $\kappa^{-1}=1000(\sin(\pi x)+1.1)$ and set $1/h=4$, $8$, $16$, $32$. Then, for $\nu=10^{-2}$, $10^{-1}$ and $10^{0}$, we present the numerical results of $\|\bm{u}-\bm{u}_{h}\|_{0}$, $\|\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h}\|_{\underline{\bm{\Sigma}}_{h}}$, $\|\underline{\bm{\sigma}}-\underline{\bm{\sigma}}_{h}\|_{0}$ and $\|p-p_{h}\|_{0}$ in Table 4.

Examples 3 and 4 do not have analytical solutions, so we do not list the convergence order as shown in the first two examples. We mention that the similar test of the profile of $\kappa^{-1}$ can be found in other literature [36, 43]. In the following two examples, a mesh $100\times 100$ is used and the data setting is designed as follows: $\Omega=(0,1)\times(0,1)$, $\bm{f}=\bm{0}$, $\nu=10^{-2}$ and $\bm{g}=(1,0)^{t}$. The Brinkman problems (3.13) are solved in a region with different contrast permeability.

Example 3**.**

For this test case, the profile of $\kappa^{-1}$ is plotted in Figure 1 with $\kappa^{-1}=10,10^{3},10^{5}$ in yellow region and $\kappa^{-1}=1$ in the purple region. ([43]).

For $\kappa^{-1}=10,10^{3},10^{5}$ in the yellow region and $\kappa^{-1}=1$ in purple region of 1, the first and the second components of the velocity obtained by MDG method with $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{2}$-$\bm{\mathcal{P}}_{1}$ element are presented in Figure 2 and Figure 3, respectively. The stress intensity and pressure profiles are showed in Figure 4 and Figure 5.

Example 4**.**

In this example, the profile of $\kappa^{-1}$ is plotted in Figure 1 with $\kappa^{-1}=10,10^{3},10^{5}$ in the purple region and $\kappa^{-1}=1$ in yellow region. ([36, 43]).

For $\kappa^{-1}=10,10^{3},10^{5}$ in the purple region and $\kappa^{-1}=1$ in yellow region of 1, the first and the second components of the velocity obtained by MDG method with $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{2}$-$\bm{\mathcal{P}}_{1}$ element are presented in Figure 6 and Figure 7, respectively. The stress intensity and pressure profiles are showed in Figure 8 and Figure 9.

From Figure 2 and Figure 3, we can see that the velocity of fluid in the purple region is faster and in the yellow region become slower as the contrast permeability increases. While, from Figure 6 and Figure 7, we can see that the velocity of fluid in the yellow region become faster and in the purple region become slower with the contrast permeability increasing. Figure 8 and Figure 9 show that the high contrast permeability gives rise to the large velocity difference, and the intensity and pressure increase rapidly.

Both Example 3 and Example 4 indicate that the contrast permeability is higher, the change of the velocity, pressure, and stress is greater. And they show the robustness, accuracy, and flexibility of the MDG method for the Brinkman problem.

6 Summary

In this paper, the mixed discontinuous Galerkin method with $\underline{\bm{\mathcal{P}}}^{\mathbb{S}}_{k+1}$-$\bm{\mathcal{P}}_{k}$ element pair is constructed and studied for solving the Brinkman equations based on the pseudostress-velocity formulation. The well-posedness of the MDG scheme is proved by the generalized Brezzi theory, and a priori error analysis is established. For any $k\geq 0$, we prove the optimal convergence order for the stress in broken $\underline{\bm{H}}(\mathbf{div})$ norm and velocity in $\bm{L}^{2}$ norm. Furthermore, the $\underline{\bm{L}}^{2}$ error estimate for the pseudostress is also investigated under certain conditions. Numerical examples confirm the theoretical results. In summary, the proposed MDG method for Brinkman equation has following main advantages: (i) it is uniformly stable and efficient from the Darcy limit to the Stokes limit; (ii) it provides accurate approximation to both the symmetric stress and the velocity; (iii) it is locally conservative for the physical quantities.

Acknowledgments. We thank Professor Suchuan Dong (Purdue University) for the helpful discussions.