The Brinkman equations (cf. Brinkman (1949)),
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.
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
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 inby , and , respectively. In particular, , and . Then, we introduce the following space
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
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
Additionally, due to the incompressibility condition, we assume that satisfies the compatibility condition , where stands for the unit outward normal on .
If , by taking , and , we could eliminate the parameter , i.e.
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)):
where the pressure can be obtained by the postprocessing formula
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
Here, the bilinear forms are defined by
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
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
We introduce a new bilinear form on , i.e.
Then, the problem (2.6) can be transformed into the following problem: Find , such that
where . By Cauchy-Schwarz inequality, we have the boundedness of .
The bilinear form satisfies
Next, we show the inf-sup condition of at continuous level.
For any , we have
Set and . Then, it is straightforward to show that
We take and , where the non-negative coefficients , and will be specified later. Then, according to Cauchy-Schwarz inequality and (2.14), one gets
From above inequality, let , , and . Then, it holds that
Here, we use the fact that due to the definition of . Then, we obtain
Next, taking , and in and , by (2.14) and the fact that , one finds
Then, we finish this proof. ∎
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
The norm of is defined by
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
where is a matrix with as its -th element. On boundary edge , we set
For any tensor-valued function and vector-valued function , a straightforward computation shows that
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
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
By Green’s formula and (3.3), we have
We define the numerical fluxes and by
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