Modelling adjacent free flow and porous media flow is important for a range of applications, e.g., transport of drugs via blood flow in vessels in biomedical engineering, and transport of pollutants via surface/groundwater flow in environmental engineering. The problem can be stated as a system of partial differential equations, with free flow governed by the Stokes equations and porous media flow governed by Darcy’s equations. The interactions at the boundary between the free flow and porous media flow regions were specified byBeavers and Joseph (1967); Saffman (1971), and were mathematically justified in Mikelic and Jäger (2000). We refer to Discacciati and Quarteroni (2009) for an overview of the model.
Well-posedness of the weak formulation of the Stokes–Darcy problem can be found in Discacciati et al. (2002) for the primal form, and in Layton et al. (2002) for the primal–mixed form. Many different finite element and mixed finite element methods have been proposed to discretize the Stokes–Darcy problem for both formulations, e.g. Discacciati et al. (2002); Layton et al. (2002); Burman and Hansbo (2005); Cao et al. (2010); Gatica et al. (2009); Márquez et al. (2015). Other devised finite element methods include discontinuous Galerkin (DG) methods Çeşmelioğlu and Rivière (2009); Girault and Rivière (2009); Girault et al. (2014); Rivière (2005); Rivière and Yotov (2005), and hybridizable discontinuous Galerkin (HDG) methods Egger and Waluga (2013); Fu and Lehrenfeld (2018); Gatica and Sequeira (2017); Igreja and Loula (2018).
We develop a numerical scheme for which the velocity field is divergence-conforming on the whole domain and for which mass is conserved pointwise. Finite element methods that satisfy these properties were proposed in Girault et al. (2014); Kanschat and Rivière (2010) where they are referred to as ‘strongly conservative’. For the Stokes region Girault et al. (2014); Kanschat and Rivière (2010)
used a divergence-conforming DG space for the velocity and a standard DG space for the pressure. In the Darcy region they used a mixed finite element method. It is well-known, however, that DG methods can be expensive due to the large number of degrees-of-freedom on a given mesh compared to other methods. To reduce the number of globally coupled degrees-of-freedom,Fu and Lehrenfeld (2018) proposed an HDG method for the Stokes region using a divergence-conforming finite element space for the velocity. Their method results in less globally coupled degrees-of-freedom compared to standard HDG methods as they only enforce continuity of the tangential direction of the facet velocity. Additionally, to reduce the problem size even further, they applied the ‘projected jumps’ method in which the polynomial degree of the tangential facet velocity is reduced by one compared to the cell velocity approximation (see also Lehrenfeld and Schöberl (2016)).
In this paper we propose an embedded–hybridized discontinuous Galerkin (EDG–HDG) finite element method of the primal–mixed formulation of the Stokes–Darcy problem on the whole domain. The EDG–HDG method uses a continuous trace velocity approximation and a discontinuous trace pressure approximation. Due to the continuous trace velocity approximation, the number of globally coupled degrees-of-freedom of the EDG–HDG method is fewer than for a traditional HDG method. However, the main motivation for an EDG–HDG discretization is not that the problem size is smaller, but that ‘continuous’ discretizations are generally better suited to fast iterative solvers. This was demonstrated for the Stokes problem in Rhebergen and Wells (2018), where CPU time and iteration count to convergence was reduced significantly compared to a hybridized method using only discontinuous facet approximations. We will show that the EDG–HDG method proposed in this paper is pointwise mass-conserving and that the resulting velocity field is divergence-conforming. We present furthermore an analysis of the proposed EDG–HDG method for the Stokes–Darcy problem, proving well-posedness, and optimal a priori
The remainder of this paper is organized as follows. In section 2 we briefly introduce the coupled Stokes–Darcy problem. The EDG–HDG method to this problem is presented in section 3. Consistency, continuity and well-posedness are shown in section 4 while the main results of this paper, an a priori error analysis, is presented in section 5. Numerical simulations support our theoretical results in section 6, and conclusions are drawn in section 7.
2 The Stokes–Darcy system
Let be a bounded polygonal domain with , boundary and boundary outward unit normal . We assume that is divided into two non-overlapping regions, and , such that and and are a union of polygonal subdomains. We denote the polygonal interface between and by . Furthermore, we define the external boundary of by , and the external boundary of by . See fig. 1.
Given the kinematic viscosity , forcing term , permeability constant and the source/sink term , the Stokes–Darcy system for the velocity field and pressure is given by
is the strain rate tensor and
is the characteristic function that has the value 1 inand 0 in . We will also frequently denote the velocity and pressure in by and , respectively, for .
denote the unit normal vector on the interface between the two domains,, pointing outwards from . On the interface we prescribe
3 The embedded–hybridized discontinuous Galerkin method
For , let be a triangulation of into non-overlapping cells . For brevity, we consider the case of matching meshes at the interface . Furthermore, let . The diameter of a cell is denoted by and denotes the maximum diameter over all cells. The outward unit normal vector on the boundary of a cell, , is denoted by . An interior facet is shared by adjacent cells, and , while a boundary facet is part of that lies on . The set and union of all facets are denoted by and , respectively. By we denote the set of all facets that lie on . For , we denote by the set of all facets that lie in . Finally, for , we denote the union of all facets in by .
We consider the following discontinuous Galerkin finite element function spaces on ,
where denotes the space of polynomials of degree on domain . On and , we consider the finite element spaces:
Note that and do not necessarily coincide on the interface . Note also that functions in are continuous on , while functions in are discontinuous on , for .
For notational purposes, we introduce the spaces , and for . Function pairs in , and , for , will be denoted by , and . Furthermore, we set .
The bilinear form is defined as
and where is a penalty parameter. The bilinear forms and are defined as:
3.3 Properties of the numerical scheme
Setting and in eq. 5 demonstrates cell-wise momentum balance eq. 1a subject to weak satisfaction of the boundary condition provided by , and a cell-wise statement of Darcy’s law eq. 1b subject to weak satisfaction of the boundary condition provided by and a Neumann boundary condition on . Setting and in eq. 5 shows that the formulation imposes normal continuity weakly across facets of the ‘numerical’ Stokes stress tensor:
where is the identity tensor. Setting and in eq. 5 and noting that , the numerical scheme imposes pointwise mass conservation, i.e.,
where is the standard -projection operator onto . Finally, setting , and with in eq. 5 and noting that on each , we find that is -conforming, i.e.,
where is the usual jump operator and the unit normal vector on .
4 Consistency, continuity and well-posedness for Stokes–Darcy
To prove consistency, continuity and stability we require extended function spaces and appropriate norms. We introduce
and . We let be the trace space of restricted to and be the trace space of restricted to . We introduce the trace operator to restrict functions in to , and the trace operators to restrict functions in to . We remark that on the interface . Where no ambiguity arises we omit the subscript when using the trace operator. For notational purposes we also introduce , and
For we denote by and the restriction of, respectively, and to .
We use various norms throughout, which are defined now. On we define
where denotes the standard -norm on domain , and
On we introduce
Note that the norms and are equivalent on , see (Wells, 2011, eq. (5.5)).
Finally, for with and , we define
We will make use of various standard estimates. In particular, use will be made of the trace inequalities for , (Di Pietro and Ern, 2012, Lemma 1.46, Remark 1.47)
and the following straightforward extensions of the continuous trace inequality (Brenner and Scott, 2010, Theorem 1.6.6)
where are independent of .
We consider each form in the definition of separately. Since on cell boundaries in , and integration by parts,
Using smoothness of , single-valuedness of , and eq. 2c, we note that
Combining with eq. 19,
Applying integration by parts, noting that on cell boundaries in , with , and on , results in
4.3 Coercivity and continuity of and
In this section we show that is coercive on for sufficiently large penalty parameter . We furthermore prove continuity of and .
Lemma 2 (Coercivity)
There exists a constant , independent of , and a constant such that for and for all ,
The result follows by definition of . ∎
Lemma 3 (Continuity)
There exists a generic constant , independent of , such that for all ,
4.4 The inf-sup condition
To present the inf-sup condition it will be convenient to split the velocity-pressure coupling term in eq. 9a as follows
The main result of this section is the following theorem.
Theorem 1 (inf-sup condition)
There exists a constant , independent of , such that for any ,
To prove this result we require the definition of two interpolation operators. For the velocity we require the following BDM interpolation operator(Hansbo and Larson, 2002, Lemma 7).
Lemma 4 (BDM interpolation operator)
If the mesh consists of triangles () or tetrahedra () there is an interpolation operator , where is the space of functions in restricted to cells in , such that for all ,
for all .
for all , where is an edge () or face () of .
, where is the usual jump operator.
with and .
We also require an interpolation operator , for example the Scott–Zhang interpolant (Brenner and Scott, 2010, Theorem 4.8.12), with the property that for , ,
where is a generic constant independent of .
Additionally, we require the following two auxiliary results.
There exists a constant , independent of , such that for any ,