A Hybrid High-Order discretisation of the Brinkman problem robust in the Darcy and Stokes limits

Abstract

In this work, we develop and analyse a novel Hybrid High-Order discretisation of the Brinkman problem. The method hinges on hybrid discrete velocity unknowns at faces and elements and on discontinuous pressures. Based on the discrete unknowns, we reconstruct inside each element a Stokes velocity one degree higher than face unknowns, and a Darcy velocity in the Raviart–Thomas–Nédélec space. These reconstructed velocities are respectively used to formulate the discrete versions of the Stokes and Darcy terms in the momentum equation, along with suitably designed penalty contributions. The proposed construction is tailored to yield optimal error estimates that are robust throughout the entire spectrum of local (Stokes- or Darcy-dominated) regimes, as identified by a dimensionless number which can be interpreted as a friction coefficient. The singular limit corresponding to the Darcy equation is also fully supported by the method. Numerical examples corroborate the theoretical results. This paper also contains two contributions whose interest goes beyond the specific method and application treated in this work: an investigation of the dependence of the constant in the second Korn inequality on star-shaped domains and its application to the study of the approximation properties of the strain projector in general Sobolev seminorms.

Introduction

In this work, we develop and analyse a novel Hybrid High-Order (HHO) method for the Brinkman problem robust across the entire range of (Stokes- or Darcy-dominated) local regimes.

Let $\Omega \subset {\mathbb{R}}^d$, $d\in \{2,3\}$, denote a bounded connected open polygonal (if $d=2$) or polyhedral (if $d=3$) set that does not have cracks, i.e., it lies on one side of its boundary $\partial \Omega$. Let two functions $\mu :\Omega \to \mathbb{R}$ and $\nu :\Omega \to \mathbb{R}$ be given corresponding, respectively, to the fluid viscosity and to the ratio between the viscosity and the permeability of the medium. In what follows, we assume that there exist real numbers $\underline{\mu },\overline{\mu }$ and $\underline{\nu },\overline{\nu }$ such that, almost everywhere in $\Omega$,

$$0<\underline{\mu }\le \mu \le \overline{\mu },0\le \underline{\nu }\le \nu \le \overline{\nu }.$$

Let $\mathit{f}:\Omega \to {\mathbb{R}}^d$ and $g:\Omega \to \mathbb{R}$ denote volumetric source terms. The Brinkman problem reads: Find the velocity $\mathit{u}:\Omega \to {\mathbb{R}}^d$ and the pressure $p:\Omega \to \mathbb{R}$ such that

$$-\mathit{\nabla }\cdot \left(2\mu {\mathit{\nabla }}_{\mathrm{s}}\mathit{u}\right)+\nu \mathit{u}+\mathit{\nabla }p=\mathit{f}\text{in }\Omega ,$$$$\mathit{\nabla }\cdot \mathit{u}=g\text{in }\Omega ,$$$$\mathit{u}=\mathbf{0}\text{on }\partial \Omega ,$$$${\int }_{\Omega }p=0,$$

where ${\mathit{\nabla }}_{\mathrm{s}}$ denotes the symmetric part of the gradient. The PDE (2) locally behaves like a Stokes or a Darcy problem depending on the value of a dimensionless parameter, which can be interpreted as a local friction coefficient. Our goal is to handle both situations robustly, while keeping the usual convergence properties of HHO methods.

The literature on the discretisation of problem (2) is vast, and giving a detailed account lies out of the scope of the present work. As noticed in [1], the construction of a finite element which is uniformly well-behaved for both the Stokes and Darcy problems is not trivial. Some choices tailored to the Stokes problem fail to convergence in the Darcy limit (as is the case for the unstabilised Crouzeix–Raviart finite element [2]), or experience a loss of convergence and, possibly, a lack of convergence for the divergence of the velocity (as is the case for the Taylor–Hood element [3] or the minielement [4]). Concerning the Crouzeix–Raviart element, a possible fix was proposed in [5] based on jump penalisation terms inspired by Discontinuous Galerkin methods. In [6], the same authors study a discretisation based on piecewise linear velocities and piecewise constant pressures for which (generalised) inf–sup stability is obtained through pressure stabilisation. Stabilised equal-order finite elements are also proposed and analysed in [7]. A generalisation of the classical minielement is studied in [8], where uniform a priori and a posteriori error estimates are derived. The use of Darcy-tailored, $H(\mathrm{div};\Omega )$-conforming finite element methods is investigated in [9], where the continuity of the tangential component of the velocity across interfaces is enforced via symmetric interior penalty terms. Finite element methods have also been developed starting from weak formulations different from the one discussed in Section 2 below. Vorticity–velocity–pressure formulations are considered, e.g., in [[10], [11]]. Finally, new generation technologies have been recently proposed for the discretisation of problem (2). We cite, in particular, the isogeometric divergence-conforming B-splines of [12], the Weak Galerkin method of [13], the two-dimensional Virtual Element methods of [[14], [15]] (see also the related work [16]), and the multiscale hybrid-mixed method of [17].

In the HHO method studied here, for a given polynomial degree $k\ge 1$, the discrete unknowns for the velocity are vector-valued polynomials of total degree $\le k$ over the mesh faces and of degree $\le l≔max(k-1,1)$ inside the mesh elements. The discrete unknowns for the pressure are scalar-valued polynomials of degree $\le k$ inside each element. Based on the discrete velocity unknowns, we reconstruct, inside each mesh element $T$: (i) a Stokes velocity inspired by [18] which yields the strain projector of degree $(k+1)$ inside $T$ when composed with the local interpolator and (ii) a Darcy velocity in the local Raviart–Thomas–Nédélec space [[19], [20]] of degree $k$. The Stokes and Darcy velocity reconstructions are used to formulate the discrete counterparts of the first and second terms in (2a). Coercivity is ensured by stabilisation terms that penalise the difference between the discrete unknowns and the interpolate of the corresponding reconstructed velocity. Owing to this finely tailored construction, the resulting method behaves robustly across the entire range of local (Stokes- or Darcy-dominated) regimes.

We carry out an exhaustive analysis of the method. We first show in Theorem 11 that the method is inf–sup stable and, based on this result, that the discrete problem is well-posed. We next prove in Theorem 12 an estimate in $h^{k+1}$ (with $h$ denoting, as usual, the meshsize) for the energy-norm of the error defined as the difference between the discrete solution and the interpolate of the continuous solution. This estimate is robust in the sense that the multiplicative constant in the right-hand side: (i) is prevented from exploding in both the Stokes- and Darcy-limits by cutoff factors; (ii) has an explicit dependence on the local friction coefficient that shows how the relative importance of the Stokes- and Darcy-contributions varies according to the local regime; (iii) does not depend on the pressure, thereby ensuring robustness when $\mathit{f}$ has large irrotational part (see [21] and references therein for further insight into this point). The Darcy velocity reconstruction in the Raviart–Thomas–Nédélec space plays a key role in achieving the aforementioned robust features while retaining optimal convergence. We point out that, to the best of our knowledge, estimates for the Brinkman problem where the various local regimes are identified by a dimensionless number are new, and they contribute to shedding new light on aspects of this problem that had often been previously treated only in a more qualitative fashion. Finally, it is worth mentioning that the theoretical results extend to the Darcy problem (corresponding to $\mu =0$ and $\underline{\nu }>0$) thanks to a stabilisation term that strengthens the coercivity norm for the Darcy term in (2a); see Remark 14 and the numerical tests in Section 5.

Besides the results specific to the Brinkman problem, this paper also contains two important contributions of more general interest. The first contribution is a study of the dependence of the constant in the second Korn inequality for polytopal domains that are star-shaped with respect to every point of a ball. We show, in particular, that this type of inequality holds uniformly inside each mesh element when considering regular mesh sequences, a key point to prove stability and error estimates for discretisation methods. The second contribution of general interest, linked to the latter point, are optimal approximation results for the strain projector, stated in Theorem 24 and Corollary 26, which extend [18, Lemma 2] to more general Sobolev seminorms. The proof hinges on the framework of [22, Section 2.1] for the study of projectors on local polynomial spaces, based in turn on the classical theory of [23].

The rest of the paper is organised as follows. In Section 2 we recall a classical weak formulation of problem (2). In Section 3 we discuss the discrete setting: mesh, local and broken polynomial spaces, and $L^2$-orthogonal projectors thereon. In Section 4 we describe the construction underlying the HHO method, formulate the discrete problem, and state the main results (whose proofs are postponed to Section 6). Numerical results are collected in Section 5. The paper is completed by an Appendix made of two sections. Appendix A.1 is dedicated to proving a uniform Korn inequality for star-shaped polytopal sets. This inequality is used in Appendix A.2 to study the approximation properties of the strain projector on local polynomial spaces for such sets. The material is structured so that multiple levels of reading are possible. Readers mainly interested in the numerical recipe and results can focus on Sections 2 to 5. Those interested in the details of the convergence analysis can additionally consult Section 6 and, possibly, Appendix A.