A hybridizable discontinuous Galerkin method for Stokes/Darcy coupling on dissimilar meshes

Abstract

1. Introduction

During the last decade, the development of new non-body-fitted numerical methods for partial differential equations (PDEs) has become of interest in the community, especially with a focus on high-order schemes. One of the most popular is the cut finite element method (CutFEM). Roughly speaking, the CutFEM method considers a background grid where the domain is immersed and a Nitsche’s approach is employed to impose the transmission conditions in the elements cut by the interface. A review can be found in Burman et al. (2015) and recent works have also proposed conservative CutFEM schemes (Larson & Zahedi, 2023; Frachon et al., 2024). CutFEM requires special quadrature rules to compute the integrals over the interface, in contrast with the recently developed |$\phi $|-FEM method (Duprez & Lozinski, 2020; Duprez et al., 2023a,b). The main idea there is to introduce an auxiliary variable that depends on the level-set function in such a way that the homogeneous Dirichlet boundary condition is automatically satisfied.

A different approach to handle transmission/boundary conditions with unfitted methods is based on a Taylor expansion of the function near the interface/boundary. In this direction, in the literature we can find two methods: the Shifted Boundary Method (SBM) (Main & Scovazzi, 2018; Atallah et al., 2022) and the Transfer Path Method (TPM) (Cockburn & Solano, 2012; Qiu et al., 2016; Oyarzúa et al., 2019). The former considers a primal formulation, and the residual of the Taylor expansion vanishes at discrete level. The TPM, on the other hand, is based on a mixed formulation where the residual of the Taylor expansion does not vanish but it involves the mixed variable that is then approximated by the numerical scheme. Our work focuses on the latter with the aim of demonstrating that the TPM can be a useful technique for handling situations where two meshes of different sizes are apart from each other. In addition, as a byproduct, our analysis also covers the case where the interface is fitted by the two meshes, allowing the presence of hanging nodes.

In several applications, the domain of interest |$\varOmega \subset \mathbb{R}^{n}$|⁠, |$n\in \{2,3\}$|⁠, is divided into subdomains where different governing equations are posed. It is not uncommon to mesh each subdomain separately using different meshsizes. For instance, in the case of solid–fluid interactions, the fluid equations are coupled to the elasticity equations via appropriate transmission conditions across an interface, and it is often desirable to have a finer mesh in the region occupied by the fluid compared to the meshsize used for discretizing the solid. When the domain of a PDE is discretized by the union of different computational subdomains, it is possible to identify two configurations. In the first one, the interface is not fitted by the triangulations, generating dissimilar meshes with gaps and overlaps appearing between the grids associated to each subdomains, as the one depicted in Fig. 1 (left). Therefore, the discrete interfaces of neighboring grids need to be properly connected. In the second configuration the interface is fitted by the grids, but it presents hanging nodes as portrayed in Fig. 1 (right). This causes a nonconformity at the intersection of the subdomains in which adjacent elements do not necessarily share a complete face or edge. This is why we prefer to consider a discontinuous Galerkin method (DG) to discretize the PDE. In particular, we focus on the hybridizable DG (HDG) method.

Fig. 1.

Left: example of dissimilar meshes separated by a uniform gap of size |$\delta $|⁠. Right: a piecewise polygonal interface separating two regions discretized by different meshes.

Open in new tabDownload slide

The HDG method introduced in Cockburn et al. (2009) discretizes a mixed formulation of the Darcy equation, where the only globally coupled degrees of freedom are those of the numerical traces on the boundaries between elements, while the remaining unknowns are obtained by solving local problems in each element. Specifically, at the continuous level, intra-element variables can be expressed in terms of inter-element unknowns by solving local problems on each element. These problems, referred to as local solvers, can be discretized using a DG method, leading to the family of methods known as HDG methods.

Furthermore, to the best of our knowledge, there are only two works that analyze the HDG method for non-conforming triangulations (Chen & Cockburn, 2012, 2014). In the first approach the authors perform an analysis for the convection-diffusion equations in nonconforming meshes. In particular, using polynomial approximations of degree |$k$| in all elements, they obtained suboptimal order of convergence |$h^{k+1/2}$| for the diffusive flux and optimal convergence |$h^{k+1}$| for the projection of the error in the scalar variable. The second approach is similar to the first one, but uses the so-called semimatching nonconforming meshes. Then, both optimal convergence for the diffusive flux and superconvergence of the projection of the error in the scalar variable is obtained.

In this work, with the aim of developing a high-order method to handle geometries with complex interfaces, we present an HDG method for coupled problems on dissimilar and nonconforming meshes. More precisely, we focus on the coupling of fluid and porous media flows across a discrete interface that does not necessarily match the true interface. To this end, we rely on the TPM (Cockburn et al., 2010b, 2012, 2014; Qiu et al., 2016) originally designed for handling boundary value problems in curved domains, but recently employed for coupling dissimilar meshes in the context of a single PDE in the entire domain (Manríquez et al., 2022; Solano et al., 2022). It is worth noting that, in the last decade, several HDG methods addressing the Stokes–Darcy coupling problem have emerged in the literature. For instance, Cesmelioglu et al. (2023), Cesmelioglu et al. (2020) present a strongly conservative discretization of the velocity–pressure formulation of the Stokes equations, coupled with Darcy-transport and Darcy equations, respectively. Likewise, the approach developed in Fu & Lehrenfeld (2018) guarantees exact mass conservation, where the Stokes flow is discretized using a divergence-conforming HDG method, while the Darcy flow is addressed by a mixed finite element method. Another contribution presented in Gatica & Sequeira (2017) employs a fully mixed formulation, selecting stress, vorticity, velocity and velocity trace as main variables in the Stokes domain, and velocity, pressure and pressure trace in the Darcy domain. HDG schemes for the coupled Navier–Stokes and Darcy equations have also been developed. For instance, in Cesmelioglu et al. (2024) a strongly conservative HDG for the time dependent case was proposed. All of these methods successfully demonstrate optimal convergence rates. Our work proposes and analyzes a new method for Stokes/Darcy coupling, following a similar scheme developed in Gatica & Sequeira (2017). More precisely, let |$\varOmega _{\mathrm{s}}$| and |$\varOmega _{\mathrm{d}}$| be bounded and simply connected polyhedral domains in |$\mathbb{R}^{n}$|⁠, |$n\in \{2,3\}$|⁠, with outward unit normal vectors |$\boldsymbol{\mathrm{n}}_{\mathrm{s}}$| and |$\boldsymbol{\mathrm{n}}_{\mathrm{d}}$|⁠, respectively, such that |${\mathcal{I}}:= \overline{\varOmega }_{\mathrm{s}} \cap \overline{\varOmega }_{\mathrm{d}}$| is the interface that separates them, and let |$\varGamma _{\mathrm{s}}:= \partial \varOmega _{\mathrm{s}}\setminus{\mathcal{I}}_{} $|⁠, |$\varGamma _{\mathrm{d}}:= \partial \varOmega _{\mathrm{d}}\setminus{\mathcal{I}}$|⁠. On |${\mathcal{I}}$|⁠, we also consider unit tangent vectors, |$\boldsymbol{m}_{1}$| when |$n=2$| and {|$\boldsymbol{m}_{1}, \boldsymbol{m}_{2}$|} when |$n=3$|⁠. The model consists of two separate groups of equations and a set of coupling terms. In the fluid region |$\varOmega _{\mathrm{s}}$|⁠, the governing equations are those of the Stokes problem, which can be written as follows:

where |$\nu>0$| is the fluid dynamic viscosity, |$\boldsymbol{\mathrm{f}}_{\mathrm{s}} \in \mathbf{L}^{2}(\varOmega _{\mathrm{s}})$| is the volumetric force acting on the fluid, |$\mathbf{u}_{\mathrm{s}}$| is the fluid velocity, |$\mathrm{L}_{\mathrm{s}}$| is the velocity gradient tensor, |$P_{\mathrm{s}}$| is the pressure, and |$\mathrm{I}$| is the |$n \times n$| identity matrix. In turn, in the porous medium region |$\varOmega _{\mathrm{d}}$| we consider the following Darcy model:

where |$\kappa $| is a tensor valued function, which describes the permeability of |$\varOmega _{\mathrm{d}}$|⁠, satisfies |$\kappa ^{\texttt{t}}=\kappa ,$| and has |$L^{\infty }(\varOmega _{\mathrm{d}})$| components, |${\mathrm{f}}_{\mathrm{d}} \in{L}^{2}(\varOmega _{\mathrm{d}})$| is a given source term such that |$\int _{\varOmega _{\mathrm{d}}}{\mathrm{f}}_{\mathrm{d}} =0$|⁠, and |$\mathbf{u}_{\mathrm{d}}$| and |${p}_{\mathrm{d}}$| denote the velocity and pressure, respectively. Also, we assume that there exist positive constants |$\underline{\kappa }$| and |$\overline{\kappa }$| such that |$\underline{\kappa }\leq \|\kappa \|_{\infty ,\varOmega _{\mathrm{d}}} \leq \overline{\kappa }$|⁠. Finally, the transmission conditions on |${\mathcal{I}}$| are given by

where |$\{\omega _{1},...,\omega _{n-1}\}$| is a set of positive frictional constants that can be determined experimentally. The first equation in (1.1c) is based on mass conservation, whereas the second one establishes the balance of normal forces and the Beavers–Joseph–Saffman law. Well-posedness of weak formulations associated to (1.1) can be found, for instance, in Gatica et al. (2009, 2011) and references therein.

We end this section by mentioning that our work continues the development of the TPM applied to dissimilar meshes started in Manríquez et al. (2022); Solano et al. (2022) and therefore we employ similar techniques. However, in this manuscript other issues arise due to the fact that two different equations are being coupled. Here, the primal variable of the Stokes system is coupled to the mixed variable of the Darcy equations and vice-versa. This adds more difficulties in designing and analyzing the scheme compared to the cases in Manríquez et al. (2022); Solano et al. (2022) where a single equation holds in the entire domain. We will come back to this point in Remark 2 after introducing the discrete setting.

The manuscript is organized as follows: in Section 2, we introduce the notation related to the discretization and transferring segments, as well as some preliminaries and definitions related to the computational domain and the approximation spaces. Next, in Section 3, the HDG method is introduced along with the proposed transmission conditions. In Section 4, we show the stability of the method and present the error estimates in Section 5. Several numerical experiments validating the performance of the method and confirming the rates of convergence are reported in Section 6. Finally, we conclude the study providing final remarks in Section 7.

2. Preliminaries

We begin by introducing some preliminary notations related to the geometric discretization, the approximation spaces and the HDG scheme. In turn, we introduce the main tools to address the discretization of the interface.

The computational domain. Let |$\varOmega ^{\mathrm{s}}_{h_{\mathrm{s}}}$| and |$\varOmega ^{\mathrm{d}}_{h_{\mathrm{d}}}$| be triangulations of the domains |$\varOmega _{\mathrm{s}}$| and |$\varOmega _{\mathrm{d}}$|⁠, with meshsizes |$h_{\mathrm{s}},\, h_{\mathrm{d}}>0$| and boundaries |$\mathcal{S}_{\mathrm{s},h_{\mathrm{s}}}, \mathcal{S}_{\mathrm{d},h_{\mathrm{d}}}$|⁠, respectively. Without loss of generality, we suppose |$h_{\mathrm{d}}\, \geq \,h_{\mathrm{s}}$| and drop the sub-index |$\star \in \{\mathrm{s},\mathrm{d}\}$| when there is no confusion; for example, we just write |$\varOmega ^{\mathrm{s}}_{h}$|⁠, and |$ \varOmega ^{\mathrm{d}}_{h}$| henceforth. We also denote the set of all faces of the triangulation |$\varOmega ^{\star }_{h}$| by |$\mathcal{E}^{\star }_{h}.$| Furthermore, since |$\overline{\varOmega }^{\,\mathrm{s}}_{h}\cap \overline{\varOmega }^{\,\mathrm{d}}_{h}$| is not necessarily equal to |${\mathcal{I}}$|⁠, then, for |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠, we define the discrete interfaces |${\mathcal{I}}^{\star }_{h}:= \mathcal{S}_{ \star , h}\setminus \varGamma _{\star ,h}$|⁠, where |$\varGamma _{\star , h}$| denotes the discretization of |$\varGamma _\star $|⁠. Bearing in mind the above, we consider outward normal vectors for the new interfaces |$\mathcal{I}^{\mathrm{s}}_{h}$| and |$\mathcal{I}^{\mathrm{d}}_{h}$|⁠, which will be denoted by |$\boldsymbol{\mathrm{n}}_{\mathrm{s}, h}$| and |$\boldsymbol{\mathrm{n}}_{\mathrm{d},h}$|⁠, respectively. The family of triangulations |$\left \{\varOmega _{h}^{\star }\right \}_{h>0}$| is assumed to be shape-regular, i.e., there exists a constant |$\kappa _{\star }>0$| such that for all elements |$K\in \varOmega _{h}^{\star }$| and all |$h>0$|⁠, |$h_{K}/\rho _{K}\leq \kappa _{\star }$|⁠, where |$h_{K}$| is the diameter of |$K$| and |$\rho _{K}$| is the diameter of the largest ball contained in |$K$|⁠. For every element |$K$|⁠, we denote by |$\boldsymbol{\mathrm{n}}_{K}$| the outward unit normal vector to |$K$|⁠, writing |$\boldsymbol{\mathrm{n}}$| instead of |$\boldsymbol{\mathrm{n}}_{K}$| when there is no confusion. In this work, we consider the two configurations depicted in Fig. 1. In the first one, a uniform gap of size |$\delta $| separates the two triangulations. In the second setting, the interface is piecewise flat and both meshes are fitted to it, but with different meshsizes.

Spaces and norms. We use the standard notation for Sobolev spaces and their associated norms and seminorms, where vector-valued functions and their corresponding spaces are denoted in bold face font, and roman font in the tensor-valued case. In addition, let |$D$| be an open bounded region of |$\mathbb{R}^{n}$| or |$\mathbb{R}^{n-1}$|⁠. We denote by |$(\cdot , \cdot )_{D}$| and |$\langle \cdot , \cdot \rangle _{\partial D}$| the |${L}^{2}(D)$| and |$ {L}^{2}(\partial D)$| inner products, respectively, with induced norms |$\| \cdot \|_{D}$| and |$\| \cdot \|_{\partial D}$|⁠. Given an integer |$k\, \geq \,0$|⁠, we use the usual notation to denote the space of polynomials of degree at most |$k$| as |$ {P}_{k}(D)$|⁠, and set |$\mathbf{P}_{k}(D)\,:= \,[ {P}_{k}(D)]^{n}$| and |$\mathrm{P}_{k}(D)\,:= \,[ {P}_{k}(D)]^{n \times n}$|⁠.

We introduce now the finite-dimensional spaces

for intra-element variables, and

for trace variables. We denote by |$ {N}_{h}^{\mathrm{d}}$| and |${\mathbf N}_{h}^{\mathrm{s}}$| the restrictions of |$ {M}_{h}^{\mathrm{d}}$| and |${\mathbf M}_{h}^{\mathrm{s}}$| to the discrete interfaces |$\mathcal{I}^{\mathrm{d}}_{h}$| and |$\mathcal{I}^{\mathrm{s}}_{h}$|⁠, respectively.

The mesh-dependent inner products are defined as

and their corresponding norms denoted by

To avoid proliferation of unimportant constants, we use the terminology |$a \lesssim b$| whenever |$a \leq C b$| and |$C$| is a positive constant independent of |$h$| and the gap between both discrete interfaces.

Transfer Paths. For |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠, we introduce a mapping |$\boldsymbol{\psi }_{\star }: {\mathcal{I}} \xrightarrow{} {\mathcal{I}}^{\star }_{h}$|⁠, such that for each point |$\boldsymbol{x} \in \mathcal{I}$|⁠, we associate a point |$\boldsymbol{x}_{\star }\,=\, \boldsymbol{\psi }_{\star }(\boldsymbol{x}) \in{\mathcal{I}}^{\star }_{h}$|⁠. We also define a mapping |$\boldsymbol{\psi }:{\mathcal{I}}^{\mathrm{d}}_{h} \xrightarrow{} {\mathcal{I}}^{\mathrm{s}}_{h}$| as |$\boldsymbol{\psi }=\boldsymbol{\psi }_{\mathrm{s}}\circ \boldsymbol{\psi }_{\mathrm{d}}^{-1}$|⁠, which means that for each |$\boldsymbol{x}_{\mathrm{d}} \in{\mathcal{I}}^{\mathrm{d}}_{h}$|⁠, we associate a point |$\boldsymbol{x}_{\mathrm{s}}=\boldsymbol{\psi }(\boldsymbol{x}_{\mathrm{d}}) \in{\mathcal{I}}^{\mathrm{s}}_{h}$|⁠. We denote by |$\sigma _{\star }(\boldsymbol{x}_{\star })$| the segment starting at |$\boldsymbol{x}_{\star }$| and ending at |$\boldsymbol{x}$|⁠, with unit tangent vector |$\boldsymbol{t}_{\star }$| and length |$|\sigma _{\star }(\boldsymbol{x}_{\star })|$|⁠. Then, keeping in mind the configuration of the interfaces, i.e., piece-wise polynomial if the meshes coincide or flat interfaces for the case with gap, as depicted in Fig. 1, it follows immediately that for each |$e\in{\mathcal{I}}^{\star }_{h}$|⁠, |$\boldsymbol{t}_{\star }=\boldsymbol{\mathrm{n}}_{\star ,h}= \boldsymbol{\mathrm{n}}_{\star }$| with |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠, and |$\boldsymbol{t}_{\mathrm{d}}=-\boldsymbol{t}_{\mathrm{s}}$|⁠. This means that the direction of the segment |$\sigma _{\star }(\boldsymbol{x}_{\star })$| must be parallel to the normals computed at its ends. Therefore, from now on, when no confusion arise, we will write |$\boldsymbol{\mathrm{n}}_{\star }$|⁠, to refer to the vector associated with |$\sigma _{\star }(\boldsymbol{x}_{\star })$|⁠, with |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠. Then, |$\sigma (\boldsymbol{x}_{\mathrm{d}})$| is the segment that starts at |$\boldsymbol{x}_{\mathrm{d}}$| and ends at |$\boldsymbol{x}_{\mathrm{s}}$|⁠, with unit tangent vector |$\boldsymbol{\mathrm{n}}_{\mathrm{d}}$| and length |$|\sigma (\boldsymbol{x}_{\mathrm{d}})|$|⁠. The segment |$\sigma (\boldsymbol{x}_{\mathrm{d}})$| is referred as the Transfer Path associated with |$\boldsymbol{x}_{\mathrm{d}}$| and is assumed to satisfy two conditions: it does not intersect the interior of another transfer path and its length |$|\sigma (\boldsymbol{x}_{\mathrm{d}})|$| is of order at most |$\max \{h_{\mathrm{s}}, h_{\mathrm{d}}\} = h_{\mathrm{d}}$|⁠.

The treatment of the pressure. Since the computational domain |$\varOmega ^{\mathrm{s}}_{h}$| does not necessarily coincide with the physical domain |$\varOmega _{\mathrm{s}}$|⁠, we adopt a similar approach to that developed in Solano & Vargas (2019) by introducing a decomposition |$ P_{\mathrm{s}}=\alpha _{\mathrm{s}}+{p}_{\mathrm{s}}$| that imposes a zero-mean condition on the pressure, with |$\alpha _{\mathrm{s}}:= \frac{1}{|\varOmega ^{\mathrm{s}}_{h}|}\int _{\varOmega ^{\mathrm{s}}_{h}}P_{\mathrm{s}}$| and |${p}_{\mathrm{s}} \in{L}^{2}_{0}(\varOmega ^{\mathrm{s}}_{h})$| (⁠|$ {L}^{2}(\varOmega ^{\mathrm{s}}_{h})$|-function with zero mean in |$\varOmega ^{\mathrm{s}}_{h}$|⁠). In turn, since |$P_{\mathrm{s}}$| will be eliminated from the system, we need to rewrite |$\alpha _{\mathrm{s}}$| in terms of |${p}_{\mathrm{s}}$|⁠. By using the fifth equation of (1.1a), we deduce that

Then, |$P_{\mathrm{s}}$| can be recovered after the approximation of |${p}_{\mathrm{s}}$| is computed.

Extrapolation operator. The region enclosed by |$\varOmega ^{\mathrm{s}}_{h}$| and |$\varOmega ^{\mathrm{d}}_{h}$| (shaded area in Fig. 1) is denoted by |$\varOmega _{h}^{\texttt{ext}}$|⁠. We notice that |$\varOmega _{h}^{\texttt{ext}}$| is not meshed and, as a consequence, we do not have an HDG approximation in there. That is why the HDG approximation of the velocity gradient |$\mathrm{L}_{\mathrm{s}}$|⁠, the pressure field |$p_{\mathrm{s}}$|⁠, and the flux |$\mathbf{u}_{\mathrm{d}}$|⁠, will be locally extrapolated from the computational domain |$\varOmega ^{\mathrm{d}}_{h} \cup \varOmega ^{\mathrm{s}}_{h}$| to |$\varOmega _{h}^{\texttt{ext}}$|⁠. More precisely, let |$q$| be a tensor, vector or scalar-valued polynomial function defined on an element |$K$| in |$\varOmega ^{\mathrm{d}}_{h} \cup \varOmega ^{\mathrm{s}}_{h}$| such that |$\overline{K}\cap \overline{\varOmega }_{h}^{\,\texttt ext}\neq \emptyset $|⁠. We define its extrapolation to |$\varOmega _{h}^{\texttt{ext}}$| as

Note that the extrapolation function |$\boldsymbol{E}_{q|_{K}}(\boldsymbol{y})$| is a function whose support includes |$\varOmega _{h}^{\texttt{ext}}$|⁠, and each element |$K$| has its own extrapolation function.

The HDG projections. Let |$(\mathrm{L}_{\mathrm{s}},\mathbf{u}_{\mathrm{s}}, {p}_{\mathrm{s}}) \in \mathrm{H}^{1}(\varOmega ^{\mathrm{s}}_{h})\times \mathbf{H}^{1}(\varOmega ^{\mathrm{s}}_{h})\times{H}^{1}(\varOmega ^{\mathrm{s}}_{h})$|⁠. For the Stokes problem, we recall from Cockburn et al. (2011) its HDG projection |$\varPi _{\mathrm{s}}(\mathrm{L}_{\mathrm{s}}, \mathbf{u}_{\mathrm{s}}, {p}_{\mathrm{s}})=(\varPi ^{\mathrm{s}}_{\mathrm{G}}\mathrm{L}_{\mathrm{s}}, \varPi ^{\mathrm{s}}_{{\mathbf V}} \mathbf{u}_{\mathrm{s}}, \varPi ^{\mathrm{s}}_{Q}p_{\mathrm{s}})$| as the element of |$\mathrm{G}_{h}^{\mathrm{s}} \times{\mathbf V}_{h}^{\mathrm{s}}{} \times Q_{h}^{\mathrm{s}}{}$| defined as follows: on an arbitrary element |$K$| of |$\varOmega ^{\mathrm{s}}_{h}$|⁠, the values of the projected function on |$K$| are determined by requiring that

where |$\tau>0$| is the stabilization parameter of the HDG method. For the case |$k=0$|⁠, the projections are well-defined considering only the second equation of (2.3b) and (2.3c). Furthermore, if |$(\mathrm{L}_{\mathrm{s}},\mathbf{u}_{\mathrm{s}},{p}_{\mathrm{s}}) \in \mathrm{H}^{l_{\sigma }+1}(\varOmega ^{\mathrm{s}}_{h})\times \mathbf{H}^{l_{\mathbf{u}_{\mathrm{s}}}+1}(\varOmega ^{\mathrm{s}}_{h}) \times{H}^{l_{\sigma }+1}(\varOmega ^{\mathrm{s}}_{h})$|⁠, for |$l_{\mathbf{u}_{\mathrm{s}}}$|⁠, |$l_{\sigma } \in [0,k]$|⁠, the above projection satisfies (cf. (Cockburn et al., 2011, Theorem 2.1)) the following properties:

for all |$K \in \varOmega _{h}^{\mathrm{s}}$|⁠, where |$\mathrm{I}$| is the identity tensor. Similarly, given |$(\mathbf{u}_{\mathrm{d}}, {p}_{\mathrm{d}}) \in \mathbf{H}^{1}{(\varOmega ^{\mathrm{d}}_{h})}\times{H}^{1}{(\varOmega ^{\mathrm{d}}_{h})}$| for the Darcy problem, we recall from Cockburn et al. (2010a) its HDG projection |$\varPi _{\mathrm{d}}( \mathbf{u}_{\mathrm{d}}, {p}_{\mathrm{d}})=(\varPi ^{\mathrm{d}}_{{\mathbf V}} \mathbf{u}_{\mathrm{d}}, \varPi ^{\mathrm{d}}_{Q}p_{\mathrm{d}}) $| as the element of |$ {\mathbf V}_{h}^{\mathrm{d}}{} \times Q_{h}^{\mathrm{d}}{}$| defined as the unique element-wise solution of

for every element |$K \in \varOmega ^{\mathrm{d}}_{h}$|⁠. Similarly to the projections in the Stokes domain, for the case |$k=0$|⁠, the projections in the Darcy domain are well-defined considering only (2.5b). Given constants |$l_{\mathbf{u}_{\mathrm{d}}}$|⁠, |$l_{{p}_{\mathrm{d}}} \in [0,k]$|⁠, if |$(\mathbf{u}_{\mathrm{d}}, {p}_{\mathrm{d}}) \in \mathbf{H}^{l_{\mathbf{u}_{\mathrm{d}}}+1}{(\varOmega ^{\mathrm{d}}_{h})}~\times ~ {H}^{l_{{p}_{\mathrm{d}}}+1}{(\varOmega ^{\mathrm{d}}_{h})}$|⁠, there hold (cf. Cockburn et al. (2010a))

for all |$K \in \varOmega _{h}^{\mathrm{d}}.$| We end this section by mentioning that |$ \boldsymbol{\mathcal{P}}^\star : \mathbf{L}^{2}(e) \xrightarrow{} {\mathbf P}_{k}(e)$| and |$ \mathcal{P}^\star : {L}^{2}(e) \xrightarrow{} P_{k}(e)$| are the respective |$\mathbf{L}^{2}$| and |$ {L}^{2}$| orthogonal projections for all facet |$e$| of |$\mathcal{E}_{h}^\star $|⁠. In abuse of notation, the global projections will be also denoted by |$\boldsymbol{\mathcal{P}}^\star $| and |$\mathcal{P}^\star $|⁠.

3. The HDG method

The HDG formulation of the coupled system (1.1a)-(1.1b) reduces to: Find |$(\mathrm{L}_{\mathrm{s}, h}, \mathbf{u}_{\mathrm{s},h}, p_{\mathrm{s},h}, \widehat{\mathbf{u}}_{\mathrm{s},h}, \mathbf{u}_{\mathrm{d},h}, p_{\mathrm{d},h}, \widehat{p}_{\mathrm{d},h}) \in \mathrm{G}_{h}^{\mathrm{s}} \times{\mathbf V}_{h}^{\mathrm{s}} \times Q_{h}^{\mathrm{s}} \times{\mathbf M}_{h}^{\mathrm{s}} \times{\mathbf V}_{h}^{\mathrm{d}} \times Q_{h}^{\mathrm{d}} \times{M}_{h}^{\mathrm{d}}$| such that

for all |$(\mathrm{G}_{\mathrm{s},h}, \mathbf{v}_{\mathrm{s},h}, q_{\mathrm{s},h}, \boldsymbol{\mu }_{\mathrm{s},h}, \mathbf{v}_{\mathrm{d},h}, q_{\mathrm{d},h}, {\mu }_{\mathrm{d},h}) \in \mathrm{G}_{h}^{\mathrm{s}}{} \times{\mathbf V}_{h}^{\mathrm{s}}{} \times Q_{h}^{\mathrm{s}} \times{\mathbf M}_{h}^{\mathrm{s}}{}\times{\mathbf V}_{h}^{\mathrm{d}}{} \times Q_{h}^{\mathrm{d}}{} \times{M}_{h}^{\mathrm{d}}{}$|⁠, where

and we recall that |$\tau $| is a positive stabilization function defined in |$\partial \varOmega ^{\mathrm{s}}_{h} \cup \partial \varOmega ^{\mathrm{d}}_{h}$|⁠, assumed to be uniformly bounded. For simplicity of the exposition, we assume |$\tau $| is constant everywhere. The above equations must be complemented with suitable transmission conditions across the interfaces |${\mathcal{I}}^{\mathrm{s}}_{h}$| and |${\mathcal{I}}^{\mathrm{d}}_{h}$|⁠, which we proceed to derive now and this constitutes the novelty of our work. Indeed, we propose the following conditions:

where |$\widetilde{\mathbf{u}}_{\mathrm{s},h}$|⁠, |$\widetilde{p}_{\mathrm{d},h}$|⁠, |$\widetilde{{\boldsymbol \sigma }}_{\mathrm{s},h} \boldsymbol{\mathrm{n}}_{\mathrm{s}} $| and |$ \widetilde{\mathbf{u}}_{\mathrm{d},h} $| stand for the approximations of |$\mathbf{u}_{\mathrm{s}}|_{{\mathcal{I}}^{\mathrm{s}}_{h}}$|⁠, |${p}_{\mathrm{d}}|_{{\mathcal{I}}^{\mathrm{d}}_{h}}$|⁠, |$(\nu \mathrm{L}_{\mathrm{s}}{}-{p}_{\mathrm{s}}\mathrm{I}) \boldsymbol{\mathrm{n}}_{\mathrm{s}}|_{{\mathcal{I}}^{\mathrm{s}}_{h}}$| and |$\mathbf{u}_{\mathrm{d}} |_{{\mathcal{I}}^{\mathrm{d}}_{h}} $|⁠, respectively, based on suitable extensions (constructed below) of |$\widehat{\mathbf{u}}_{\mathrm{s},h}$|⁠, |$\widehat{p}_{\mathrm{d},h}$|⁠, |$\widehat{{\boldsymbol \sigma }}_{\mathrm{s},h}{} \boldsymbol{\mathrm{n}}_{\mathrm{s}}$| and |$\widehat{\mathbf{u}}_{\mathrm{d},h}$| outside their corresponding computational domains. More precisely, employing the transferring technique of Cockburn & Solano (2012) (see also Cockburn et al. (2014); Manríquez et al. (2022); Solano et al. (2022)), the tilde variables are constructed as follows: let |$\boldsymbol{x}_{\star } \in{\mathcal{I}}_{h}^{\star }$| and its corresponding point |$\boldsymbol{x} \in{\mathcal{I}}$|⁠. Integrating the first equation of (1.1a) along the transferring path |${\sigma _{\mathrm{s}}(\boldsymbol{x}_{\mathrm{s}})}$| and using the first equation of (1.1c), we obtain

where |$\mathbf{u}_{\mathrm{s}}\circ \boldsymbol{\psi }_{\mathrm{d}}^{-1} (\boldsymbol{x}_{\mathrm{d}})= \mathbf{u}_{\mathrm{s}}(\boldsymbol{x}_{\mathrm{s}}) +|\sigma _{\mathrm{s}}(\boldsymbol{x}_{\mathrm{s}})|\int _{0}^{{1}}\mathrm{L}_{\mathrm{s}}(\boldsymbol{y}_{\mathrm{s}}(t))\boldsymbol{\mathrm{n}}_{\mathrm{s}} dt$|⁠, being |$\boldsymbol{y}_{\mathrm{s}}(t)=\boldsymbol{x}_{\mathrm{s}}+(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{s}})t$| with |$t\in [0,1]$| the parametrization of |${\sigma _{\mathrm{s}}(\boldsymbol{x}_{\mathrm{s}})}$|⁠. Similarly, integrating the second equation of (1.1b) along the connecting segment |${\sigma _{\mathrm{d}}(\boldsymbol{x}_{\mathrm{d}})}$| and using the second equation of (1.1c), it follows that

where |${p}_{\mathrm{d}}\circ \boldsymbol{\psi }_{\mathrm{s}}^{-1} (\boldsymbol{x}_{\mathrm{s}})={p}_{\mathrm{d}}(\boldsymbol{x}_{\mathrm{d}}) -|\sigma _{\mathrm{d}}(\boldsymbol{x}_{\mathrm{d}})|\int _{0}^{{1}}\kappa ^{-1} \mathbf{u}_{\mathrm{d}}(\boldsymbol{y}_{\mathrm{d}}(t)) \cdot \boldsymbol{\mathrm{n}}_{\mathrm{d}} dt$|⁠, being |$\boldsymbol{y}_{\mathrm{d}}(t)=\boldsymbol{x}_{\mathrm{d}}+(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{d}})t$| with |$t\in [0,1]$| the parametrization of |${\sigma _{\mathrm{d}}(\boldsymbol{x}_{\mathrm{d}})}$|⁠.

Hence, motivated by these expressions and based on the form of the HDG numerical fluxes (3.2a) and (3.2b), we define

where |$\alpha _{\mathrm{s},h}\,:= \,\dfrac{1}{|\varOmega _{\mathrm{s}}|}\int _{\varOmega _{\mathrm{s}}\setminus \varOmega ^{\mathrm{s}}_{h}}\boldsymbol{E}_{p_{\mathrm{s},h}}$|⁠, and |$\boldsymbol{E}$| denotes the local extrapolation defined in (2.2).

In the particular case of matching interfaces, namely |${\mathcal{I}}= \overline{\varOmega }^{\mathrm{s}}_{h} \cap \overline{\varOmega }^{\mathrm{d}}_{h}= {\mathcal{I}}^{\mathrm{s}}_{h}={\mathcal{I}}^{\mathrm{d}}_{h}$|⁠, the transmission conditions (3.3) become

and the resulting HDG formulation is very similar to the one presented in Gatica & Sequeira (2017).

4. Stability analysis

In this section we show a stability estimate associated with (3.1) and (3.2). For that, we recall some important estimates and assumptions required to carry out our analysis.

Further notation and auxiliary estimates. Let |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠. Given a face |$e \in{\mathcal{I}}^{\star }_{h}$| belonging to the element |${K}_{e} \in \varOmega _{h}^{\star }$|⁠, we define the extrapolation patch as |$ {K}_{e}^{\texttt{ext}}:= \big \{\boldsymbol{x}+\boldsymbol{\mathrm{n}}_{\star }t: 0\leq t \leq |\sigma _{\star }(\boldsymbol{x})|,\,\boldsymbol{x}\in e\big \} $|⁠, and denote by |$h_{e}^{\perp }$| (resp. |$\delta _{e}$|⁠) the largest distance of a point inside |$K_{e}$| (resp. |$K_{e}^{\texttt{ext}}$|⁠) to the plane determined by the face |$e$|⁠. In other words, |$ h_{e}^{\perp }=\max _{\boldsymbol{x}\in K_{e}}|\mathrm{dist}(\boldsymbol{x},e)|, \quad \delta _{e}=\max _{\boldsymbol{x} \in e}|\sigma _{\star }(\boldsymbol{x})|$|⁠, where |$\mathrm{dist}(\boldsymbol{x},e)$| denotes the distance from |$\boldsymbol{x}$| to the face |$e$|⁠. We note that |$\delta _{e}$| is a measure of the local size of the gap and |$\delta := \max _{e}\delta _{e}$| is an upper bound of the size of the gap. We define the ratio |$r_{e}:= \delta _{e}/h_{e}^{\perp }$| and, for |$e \in{\mathcal{I}}^{\mathrm{s}}_{h} \cup{\mathcal{I}}^{\mathrm{d}}_{h}$|⁠, |$\boldsymbol{\mathcal{N}}^{k}\,:= \,\big \{\boldsymbol{\mathrm{q}} \in \mathbf{P}_{k}(K_{e}^{\texttt{ext}}),\quad \boldsymbol{\mathrm{q}} \cdot \boldsymbol{ \mathrm{n}}_{e}\neq 0 \text{ on each} e \subset \partial{K}_{e}^{\texttt{ext}}\big \}$|⁠, where we denoted by |$\boldsymbol{\mathrm{n}}_{e}$| the interior normal vector to |$K_{e}^{\texttt{ext}}$| along the face |$e$|⁠, that is, the exterior normal vector to |$K_{e}$| pointing in the direction of |$K_{e}^{\texttt{ext}}$|⁠. We can then introduce the constants

As proved in Cockburn et al. (2014, Lemma A.2), these constants are independent of the meshsize, but depend on the polynomial degree |$k$|⁠. The superscripts in |$ C_{e}^{\texttt{ext}}$| and |$ C_{e}^{\texttt{inv}}$| refer to an extrapolation constant and an inverse inequality constant.

On the other hand, proceeding as in Cockburn et al. (2014), we introduce the following auxiliary functions: let |$e \in{\mathcal{I}}_{h}^{\star }$| that belongs to |$K_{e}$| and |$K_{e}^{\texttt{ext}}.$| For a function |$\boldsymbol{\mathrm{q}}$|⁠, we define

for |$\star \in \{\mathrm{s},\mathrm{d} \}$|⁠, where |$\boldsymbol{x}_{\mathrm{d}} \in e$| and |$\boldsymbol{x}_{\mathrm{s}} \in{\mathcal{I}}^{\mathrm{s}}_{h}$| are connected by the segment |$\sigma (\boldsymbol{x}_{\mathrm{d}})$|⁠. They satisfy (cf. (Cockburn et al., 2014, Lemma 5.2)),

Another important tool in the analysis of this method, which is based on the Taylor series expansion of a function defined on |${\mathcal{I}}$| around a point |${\mathcal{I}}_{h}^{\star }$|⁠, is the following lemma ((Solano et al., 2022, Lemma 2.1) or (Manríquez et al., 2022, Lemma 2)).

In addition, we recall the discrete trace inequality (cf. (Di Pietro & Ern, 2012, Lemma 1.46)), which establishes that if |$\phi $| is a scalar, vector or tensor-valued polynomial in |$K_{e}$|⁠, there exists a positive constant |$C_{e}^{\mathrm{tr}}$|⁠, independent of the meshsize but depends on the polynomial degree, such that

We stress that the identities and inequalities established throughout this section hold true for the tensor, vector or scalar-valued cases as required.

Assumptions. Here we outline a set of assumptions under which our analysis holds. More precisely, we assume that

• (A.1) |$\varOmega ^{\mathrm{s}}_{h} \cap \varOmega ^{\mathrm{d}}_{h}= \emptyset $|⁠,

• (A.2) the mappings |$\boldsymbol{\psi }_{\star }: {\mathcal{I}} \xrightarrow{} {\mathcal{I}}_{\star }$| for |$\star \in \{\mathrm{s},\mathrm{d}\}$|⁠, and |$\boldsymbol{\psi }: {\mathcal{I}}^{\mathrm{d}}_{h} \xrightarrow{} {\mathcal{I}}^{\mathrm{s}}_{h}$| are bijections,

• (A.3) |$4\,{\max \big \{1,\nu ^{-1},2\sum ^{n-1}_{\ell =1}\omega ^{-2}_{\ell }\big \}}\tau ^{-1/2}C_{\delta _{\star },h} \max _{e \in{\mathcal{I}}_{h}^{\star }} \big(\delta _{e}^{-10/14} \big)+{8}\max _{e \in{\mathcal{I}}_{h}^{\star }} \big(\widehat{C}_{e}^{\star } \delta _{e}^{12/7}h_{e}^{-3}(C_{e}^{\texttt{ext}})^{2}\tau ^{-1} \big)^{}\leq \, \frac{1}{4},$| where |$C_{\delta _{\star },h}$| depends on |$h$| and |$\delta $| (cf. Lemma 2), and |$\widehat{C}_{e}^{\star }$| is a positive constant appearing in the proof of Lemma 4,

• (A.4) |$C_{\delta ,h}\,:= \,\widetilde{C}_{1}C^{\mathcal{S}}_{\delta ,h}+\left ((C^{p}_{\delta ,h})^{2}+\widetilde{C}_{2}\right )\widetilde{C}_{is}$| is small enough, where |$C^{\mathcal{S}}_{\delta ,h}\,:= \,h^{2}+ \big(C^{u_{\mathrm{d}}, L_{\mathrm{s}}}_{\delta ,h}\big)^{2}+ \big(C^{\widehat{u}_{\mathrm{s}}, \widehat{p}_{\mathrm{d}}}_{\delta ,h}\big)^{2}$|⁠,

  $$ \begin{align*} \widetilde{C}_{1}&:= \max\big\{\nu, 1\big\}\displaystyle\sum_{\star \in \{\mathrm{s}, \mathrm{d}\}}\Bigg(4 \max_{e \in{\mathcal{I}}_{h}^{\star}}\left(\left(C^{\mathrm{tr}}_{e}\right)^{2}\delta_{e}^{2/7}h_{e}^{-1}\tau\right) +\displaystyle\max_{e\in{\mathcal{I}}_{h}^{\star}}\Big(\frac{(C_{e}^{\mathrm{tr}})^{-2}{h_{e}^{-\beta}}}{2}\Big)\Bigg),\\ \widetilde{C}_{2}&:= \frac{8}{\nu}\max_{e \in{\mathcal{I}}^{\mathrm{s}}_{h}} \left(\widehat{C}_{e}^{\mathrm{s}} \delta_{e}^{12/7}h_{e}^{-3}(C_{e}^{\texttt{ext}})^{2}\tau^{-1}\right)^{}+4\frac{C_{\alpha}^{2}}{\nu}\max_{e \in{\mathcal{I}}^{\mathrm{s}}_{h}}\left(\delta_{e}^{-2/7}\right),\\ C^{u_{\mathrm{d}}, L_{\mathrm{s}}}_{\delta,h}&\,:= \,\max\big\{\nu^{-1}, 1\big\}\displaystyle\sum_{\star \in \{\mathrm{s}, \mathrm{d}\}} \Bigg( \max_{e \in{\mathcal{I}}_{h}^{\star}} \left(\delta_{e}h_{e}^{-3/2}C^{\texttt{ext}}_{e} \right)+\max_{e \in{\mathcal{I}}_{h}^{\star}} \left(\delta_{e}^{1/2} C_{\delta_{\star,h}}^{1/2}\right)\Bigg),\\ C^{\widehat{u}_{\mathrm{s}}, \widehat{p}_{\mathrm{d}}}_{\delta,h}&\,:= \,\max\big\{\nu^{-1}, 1\big\}\displaystyle\sum_{\star \in \{\mathrm{s}, \mathrm{d}\}} \Bigg(\max_{e \in{\mathcal{I}}_{h}^{\star}}\left(\delta_{e}^{6/7}\tau^{-1/2}\right)+\max_{e \in{\mathcal{I}}_{h}^{\star}}\left(\delta_{e}^{5/14}\tau^{-1/2}\right)\Bigg)+{\delta},\\ C^{p}_{\delta,h}&\,:= \,\nu^{-1}\max_{e \in{\mathcal{I}}^{\mathrm{s}}_{h}} \left(\delta_{e}h_{e}^{-3/2}C^{\texttt{ext}}_{e} \right)+\nu^{-1}C_{\alpha}, \end{align*} $$

  |$\widetilde{C}_{is}$| is related to an inf-sup condition of Lemma 6, |$\beta $| is a non-negative parameter whose range will be determined later, and |$C_{\alpha }$| is positive constant, which will appear in the proof of Lemma 4.

Let us briefly comment on these assumptions. Note that Assumptions (A.1) and (A.2) hold true, for instance, in the illustrations of Fig. 1. In particular, the Assumption (A.1) ensures that there is no overlap between subdomains, and is intended to simplify the analysis, whereas the Assumption (A.2) is the key to ‘tie’ the interfaces |${\mathcal{I}}_{h}^{\mathrm{s}}$| and |${\mathcal{I}}_{h}^{\mathrm{d}}$| that cause the gap. The remaining assumptions are technical and define the relationship between the gap size |$\delta $| and the mesh size |$h$|⁠, which is essential for ensuring the convergence and optimality of the method. For example, (A.3) is always satisfied for |$h$| small enough if |$\delta \,\lesssim \, h^{7/4}$|⁠. To analyze the feasibility of other assumptions, let us write |$\delta =C_{\mathrm{g}}h^{1+\gamma }$| with |$ C_{\mathrm{g}}\geq 0$| and |$\gamma>0$| constants independent of the meshsize. Assumption (A.4) is satisfied for all |$\gamma> 3/4$| and |$\beta \in [0,2] \cap [0, 2\gamma -1)$|⁠, if |$h$| is small enough, as we will explain in Corollaries 2 and 3. These are the strongest assumptions, since they indicate that our analysis holds if the gap size is at most of order |$h^{7/4}$|⁠, however we will present numerical evidence suggesting that the method is still optimal when the gap is of order |$h$|⁠. Finally, we highlight that the remaining constants are defined in the subsequent results presented below. In turn, in order to begin with the analysis, we establish the following result.

In order to use this analysis to establish both, well-posedness and error bounds, we consider the problem (3.1), but (3.1a) and (3.1g) are replaced by

where |$\mathrm{J}_{\mathrm{s}} \in \mathrm{L}^{2}(\varOmega ^{\mathrm{s}}_{h})$|⁠, |$\mathbf{J}_{\mathrm{d}} \in \mathbf{L}^{2}(\varOmega ^{\mathrm{d}}_{h})$| are given functions orthogonal to polynomials of degree |$k-1$| and (3.3) is replaced by

for all |$({\mu }_{\mathrm{d},h}, \boldsymbol{\mu }_{\mathrm{s},h} ) \in{N}_{h}^{\mathrm{d}}\times{\mathbf N}_{h}^{\mathrm{s}}$|⁠, where |$j_{\mathrm{d}}^{\mathrm{n} c}$| and |$\mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}$| are given functions associated with the nonconformity that occurs at the interface, belonging to |$ {L}^{2}({\mathcal{I}}^{\mathrm{d}}_{h})$| and |${\mathbf L}^{2}({\mathcal{I}}^{\mathrm{s}}_{h})$|⁠, respectively. Similarly, |$j_{\mathrm{d}}^{\delta }$| and |$\mathbf{j}_{\mathrm{s}}^{\delta }$| are associated with the gap between the discrete interfaces |$\mathcal{I}_{h}^{\mathrm{s}}$| and |$\mathcal{I}_{h}^{\mathrm{d}}$|⁠, also belonging to |$ {L}^{2}({\mathcal{I}}^{\mathrm{d}}_{h})$| and |${\mathbf L}^{2}({\mathcal{I}}^{\mathrm{s}}_{h})$|⁠, respectively. We emphasize that the motivation to introduce these given functions is only to unify the well-posedness and error analysis in a compact way. In particular, as we will see below in Corollary 1, to show well-posedness, |$\mathrm{J}_{\mathrm{s}}$|⁠, |$\mathbf{J}_{\mathrm{d}}$|⁠, |$j_{\mathrm{d}}^{\mathrm{n} c}$|⁠, |$\mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}$|⁠, |$j_{\mathrm{d}}^{\delta }$| and |$\mathbf{j}_{\mathrm{s}}^{\delta }$| will be zero. On the other hand, as we will see in Section 5, |$\mathrm{J}_{\mathrm{s}}$| and |$\mathbf{J}_{\mathrm{d}}$| are associated with projection errors when proving the error bounds, and will appear in the terms on the right-hand side of (5.1a) and (5.1f). Similarly, |$j_{\mathrm{d}}^{\mathrm{n} c}$|⁠, |$\mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}$|⁠, |$j_{\mathrm{d}}^{\delta }$| and |$\mathbf{j}_{\mathrm{s}}^{\delta }$| will appear in the terms on the right-hand side of (5.4a) and (5.4b).

4.1 An energy argument

Before presenting the energy estimate, we proceed to deduce how the transmission conditions (3.3) connect |$\langle \widehat{{\boldsymbol \sigma }}_{\mathrm{s},h}{}\boldsymbol{\mathrm{n}}_{\mathrm{s}}, \widehat{\mathbf{u}}_{\mathrm{s},h} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}$| and |$\langle \widehat{\mathbf{u}}_{\mathrm{d},h} \cdot \boldsymbol{\mathrm{n}}_{\mathrm{d}},\widehat{p}_{\mathrm{d},h}\rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}$|⁠. To this end, we define |$\mathbb{T}=-\langle \widehat{{\boldsymbol \sigma }}_{\mathrm{s},h}{}\boldsymbol{\mathrm{n}}_{\mathrm{s}}, \widehat{\mathbf{u}}_{\mathrm{s},h} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}+\langle \widehat{\mathbf{u}}_{\mathrm{d},h} \cdot \boldsymbol{\mathrm{n}}_{\mathrm{d}},\widehat{p}_{\mathrm{d},h}\rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}$| and we write it in terms associated with the mismatch between |${\mathcal{I}}^{\mathrm{s}}_{h}$| and |${\mathcal{I}}^{\mathrm{d}}_{h}$|⁠. More precisely, we prove the following lemma.

In what follows, we define

and provide an upper bound for this energy term |$ \mathcal{Q}(\mathrm{L}_{\mathrm{s}, h}, \boldsymbol{\mathrm{u}}_{\star ,h}, \widehat{\mathbf{u}}_{\mathrm{s},h},p_{\mathrm{d},h},\widehat{p}_{\mathrm{d},h})$|⁠. This bound depends, in addition to the sources, on the norms of the approximations of the velocity |$\mathbf{u}_{\mathrm{s},h}$|⁠, pressures |$p_{\mathrm{s},h}$|⁠, and |$p_{\mathrm{d},h}$|⁠.

Our next goal is to provide an estimate for the |$ {L}^{2}$|-norm of |$\mathbf{u}_{\mathrm{s},h}$|⁠, |$p_{\mathrm{s},h}$| and |$p_{\mathrm{d},h}$|⁠. To bound |$\|\mathbf{u}_{\mathrm{s},h}\|_{\varOmega ^{\mathrm{s}}_{h}}$| and |$\|p_{\mathrm{d},h}\|_{\varOmega ^{\mathrm{d}}_{h}}$| we employ a duality argument, whereas for |$\|p_{\mathrm{s},h}\|_{\varOmega ^{\mathrm{s}}_{h}}$| we use an inf-sup condition.

4.2 A duality argument

In order to estimate |$\|\mathbf{u}_{\mathrm{s},h}\|_{\varOmega ^{\mathrm{s}}_{h}}+\|p_{\mathrm{d},h}\|_{\varOmega ^{\mathrm{d}}_{h}}$|⁠, we now proceed as in Solano et al. (2022), Manríquez et al. (2022), Gatica & Sequeira (2017) and incorporate a suitable auxiliary problem. More precisely, in what follows we consider the continuous problem (1.1a)-(1.1b)-(1.1c) with sources given by |$\boldsymbol{\mathrm{f}}_{\mathrm{s}} := \varTheta _{\mathrm{s}} \in \mathbf{L}^{2}(\varOmega _{\mathrm{s}})$| and |${\mathrm{f}}_{\mathrm{d}} := \varTheta _{\mathrm{d}} \in{L}^{2}(\varOmega _{\mathrm{d}})$|⁠, that is:

In addition, we proceed to carry out the same decomposition performed for the pressure of the continuous problem, that is |$\widetilde{\varphi }_{\mathrm{s}}=\widetilde{\alpha }_{\mathrm{s}}+\varphi _{\mathrm{s}}$|⁠, where |$\varphi _{\mathrm{s}} \in{L}^{2}_{0}(\varOmega ^{\mathrm{s}}_{h})$|⁠, and |$\widetilde{\alpha }_{\mathrm{s}}=\frac{-1}{|\varOmega _{\mathrm{s}}|}\int _{\varOmega _{\mathrm{s}}\setminus \varOmega ^{\mathrm{s}}_{h}}\varphi _{\mathrm{s}}$|⁠, is the constant similar to (2.1). Furthermore, suppose that elliptic regularity holds, that is,

We are now in position to establish the main result of this section.

4.3 Semi-aligned discrete interfaces

With the aim of improving the estimate given in Theorem 1, in this section, we consider a specific structure of nonconforming meshes, where the discrete interfaces satisfy the following: if |$\boldsymbol{v}_{\mathrm{d}}$| is a vertex in |${\mathcal{I}}_{h}^{\mathrm{d}}$|⁠, then |$\boldsymbol{\psi }(\boldsymbol{v}_{\mathrm{d}})$| is a vertex |$\boldsymbol{v}_{\mathrm{s}}$| in |${\mathcal{I}}_{h}^{\mathrm{s}}$| (see Fig. 2). In this way, we refer to |${\mathcal{I}}_{h}^{\mathrm{d}}$| and |${\mathcal{I}}_{h}^{\mathrm{s}}$| as semi-aligned discrete interfaces.

Fig. 2.

Example of semi-aligned discrete interfaces separated by a uniform gap of size |$\delta $|⁠.

Open in new tabDownload slide

Under this condition and the fact that in the error analysis |$ \mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}:= (\mathcal{P}^{\mathrm{d}}p_{\mathrm{d}}-{p}_{\mathrm{d}}) \circ \boldsymbol{\psi }^{-1} \boldsymbol{\mathrm{n}}_{\mathrm{d}}$| and |$ j_{\mathrm{d}}^{\mathrm{n} c} := (\boldsymbol{\mathcal{P}}^{\mathrm{s}}\mathbf{u}_{\mathrm{s}}-\mathbf{u}_{\mathrm{s}}) \circ \boldsymbol{\psi } \cdot \boldsymbol{\mathrm{n}}_{\mathrm{s}}$| (see equation (5.4) in Lemma 8), it is easy to see that |$ \,\langle{j_{\mathrm{d}}^{\mathrm{n} c}}, {\mu }_{\mathrm{d},h} \rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}$| in (4.12a) vanishes. In fact, since |$h_{\mathrm{s}}\leq h_{\mathrm{d}}$|⁠, the term |$ {\mu }_{\mathrm{d},h} \boldsymbol{\mathrm{n}}_{\mathrm{s}}\circ \boldsymbol{\psi }^{-1}$| belongs to |$ {\mathbf N}_{h}^{\mathrm{s}}$| and thus |$\langle j_{\mathrm{d}}^{\mathrm{n} c}, {\mu }_{\mathrm{d},h} \rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}= \langle (\boldsymbol{\mathcal{P}}^{\mathrm{s}}\mathbf{u}_{\mathrm{s}}-\mathbf{u}_{\mathrm{s}}), {\mu }_{\mathrm{d},h} \boldsymbol{\mathrm{n}}_{\mathrm{s}}\circ \boldsymbol{\psi }^{-1} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}=0$|⁠. Furthermore, in what follows we show that under this configuration |$\langle \mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}, \boldsymbol{\mathcal{P}}^{\mathrm{d}}(\boldsymbol{\phi }_{\mathrm{s}} \circ \boldsymbol{\psi })\circ \boldsymbol{\psi }^{-1} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}$| and |$\langle j_{\mathrm{d}}^{\mathrm{n} c}, {\mathcal{P}}^{\mathrm{d}}\varphi _{\mathrm{d}} \rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}$| also vanish. Indeed,

On the other hand, taking into account that |$h_{\mathrm{s}} \leq h_{\mathrm{d}}$|⁠, we have

since |${\mathcal{P}}^{\mathrm{d}}\varphi _{\mathrm{d}} \circ \boldsymbol{\psi }^{-1} \boldsymbol{\mathrm{n}}_{\mathrm{s}} \in{\mathbf M}_{h}^{\mathrm{s}}$|⁠. All these identities imply an improvement of the estimate of Lemma 5. In fact, we recall that |$\mathbb{B}_{3}\,= \ - \langle \mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}, \boldsymbol{\mathcal{P}}^{\mathrm{s}}\boldsymbol{\phi }_{\mathrm{s}} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}- \langle \mathbf{j}_{\mathrm{s}}^{\delta }, \boldsymbol{\mathcal{P}}^{\mathrm{s}}\boldsymbol{\phi }_{\mathrm{s}} \rangle _{{\mathcal{I}}^{\mathrm{s}}_{h}}$|⁠. Then, we have

which implies that

On the other hand, since |$\langle j_{\mathrm{d}}^{\mathrm{n} c}, {\mathcal{P}}^{\mathrm{d}}\varphi _{\mathrm{d}} \rangle _{{\mathcal{I}}^{\mathrm{d}}_{h}}=0$|⁠, by applying the continuous trace inequality and (4.19), it follows that

We now have the following result.

5. Error analysis

Our first goal in this section is to derive the error estimates of the proposed method. We employ the stability estimate deduced in previous sections. In what follows, we introduce the projection of the errors, namely |$\varepsilon ^{\mathrm{L}_{\mathrm{s}}}:= \varPi ^{\mathrm{s}}_{\mathrm{G}}\mathrm{L}_{\mathrm{s}}-\mathrm{L}_{\mathrm{s}, h}$|⁠, |$\varepsilon ^{\mathbf{u}_{\star }}:= \varPi ^{\star }_{{\mathbf V}} \mathbf{u}_{\star }-\mathbf{u}_{\star ,h}$|⁠, |$\varepsilon ^{\widehat{\mathbf{u}}_{\star }}:= \boldsymbol{\mathcal{P}}^{\star }_{}\mathbf{u}_{\star }-\widehat{\mathbf{u}}_{\star ,h}$|⁠, |$\varepsilon ^{{p}_{\star }}:= \varPi ^{\star }_{P}p_{\star }-p_{\star ,h}$|⁠, |$\varepsilon ^{{\hat{p}}_{\mathrm{d}}}:= \mathcal{P}^{\mathrm{d}}p_{\mathrm{d}}-\widehat{p}_{\mathrm{d},h}$|⁠, and |$\varepsilon ^{\alpha _{\mathrm{s}}}:= \varPi ^{\mathrm{s}}_{Q}\alpha _{\mathrm{s}}-\alpha _{\mathrm{s},h}$|⁠. In turn, the error of the projections are given by |$\mathrm{I}_{\mathrm{L}}^{\mathrm{s}}:= \mathrm{L}_{\mathrm{s}}{}-\varPi ^{\mathrm{s}}_{\mathrm{G}}\mathrm{L}_{\mathrm{s}}$|⁠, |$\mathrm{I}_{\mathbf{u}}^{\star }:= \mathbf{u}_{\star }-\varPi ^{\star }_{{\mathbf V}} \mathbf{u}_{\star }$|⁠, |$\mathrm{I}_{p}^{\star }:= p_{\star }-\varPi ^{\star }_{P}p_{\star }$|⁠, and |$\mathrm{I}_{\alpha }^{\mathrm{s}}:= \alpha _{\mathrm{s}}-\varPi ^{\mathrm{s}}_{Q}\alpha _{\mathrm{s}}$|⁠, and note that |$\mathrm{I}_{\alpha }^{\mathrm{s}}=0$|⁠.

We observe that the above equations are similar to those of the HDG scheme (3.1), where |$\mathrm{I}_{\mathrm{L}}^{\mathrm{s}}$|⁠, |$\mathrm{I}_{\mathbf{u}}^{\mathrm{d}}$|⁠, |$0$| and |$\boldsymbol 0$| play the role of |$\mathrm{J}_{\mathrm{s}}$|⁠, |$\mathbf{J}_{\mathrm{d}}$|⁠, |${\mathrm{f}}_{\mathrm{d}}$| and |$\boldsymbol{\mathrm{f}}_{\mathrm{s}}$|⁠, respectively. Moreover,

Hence, we consider the result of Theorem 1 applied to this context. More precisely, we notice that

and

and observe that the first terms in (5.8) and (5.9) can be bounded using the approximation properties of the |$ {L}^{2}$|-projection over |$ {N}_{h}^{\mathrm{d}}$| and |${\mathbf N}_{h}^{\mathrm{s}}$|⁠, respectively. In fact, we recall from Di Pietro & Ern (2012, Lemma 1.59) that there exist positive constants |$C_{\mathrm{n} \mathrm{s}}$| and |$C_{\mathrm{n} \mathrm{d}}$|⁠, independent of |$h$|⁠, such that

We stress that these constants take into account the nonconformity between the computational interfaces. Thus, we have

On the other hand, following the proof of Solano et al. (2022, Theorem 4.2) we deduce from (4.3), a scaling argument to bound the |$ {L}^{2}\big ({\mathcal{I}}_{h}^{\star }\big )-$|norm in terms of its |$ {L}^{2}(\varOmega _{h}^{\star })-$|norm, and Lemma 1, that

and

Next, by Assumptions (A), the fact that |$(\delta _{e}^{4}h_{e}^{-4}h^{{\beta }})$| and |$(\delta _{e}^{}h_{e}^{-1}h^{{\beta }})$| are bounded, with |$\star \in \{\mathrm{s}, \mathrm{d}\}$|⁠, it follows that

Then, bearing in mind the above, the estimates from Theorem 1 applied to (5.1) become

and

which, along with the properties of the HDG projectors (cf. (2.4)–(2.6)), we obtain the following result.

We now explain the consequences of this corollary for some particular cases. In what follows, in abuse of terminology, we will employ the nomenclature ‘hanging nodes’, even in the presence of a gap, when there is not a one-to-one correspondence between the nodes of the discrete interfaces.

• (C.1) No-gap and no hanging nodes. In this case, |$C_{\mathrm{n}\mathrm{s}}=C_{\mathrm{n} \mathrm{d}}=C_{\mathrm{g}}=0$| and our result shows optimal order of convergence of |$h^{k+1}$| for all the variables and order |$h^{k+2}$| for the projection of the errors |$\varepsilon ^{\mathbf{u}_{\mathrm{s}}}$| and |$\varepsilon ^{p_{\mathrm{d}}}$|⁠, as expected.

• (C.2) No-gap and hanging nodes. In this case, |$C_{\mathrm{g}}=0$|⁠, but |$C_{\mathrm{n}\mathrm{s}} \neq 0$| and |$C_{\mathrm{n} \mathrm{d}}\neq 0$|⁠. In this situation, since |$C_{\mathrm{g}}=0$|⁠, we can take |$\beta = 2$|⁠. Therefore, we obtain optimal order of convergence of |$h^{k+1}$| for the variables |$\mathbf{u}_{\mathrm{d}}$|⁠, |$\mathrm{L}_{\mathrm{s}}$| and |$p_{\mathrm{s}}$|⁠, but suboptimal order of |$h^{k+1/2}$| for the projection of the errors |$\varepsilon ^{\mathbf{u}_{\mathrm{s}}}$| and |$\varepsilon ^{p_{\mathrm{d}}}$| and also for the errors in |$\mathbf{u}_{\mathrm{s}}$| and |$p_{\mathrm{d}}$|⁠.

• (C.3) Gap |$\delta $| of order |$h^{2}$| and no hanging nodes. Here |$\gamma =1$|⁠, which implies that |$\beta =1-\epsilon $| for all |$\epsilon>0$|⁠, whereas the nonconformity constants |$C_{\mathrm{n} \mathrm{s}}$| and |$ C_{\mathrm{n} \mathrm{d}}$| vanish. This yields, for all the variables, order of convergence |$h^{k+1-\epsilon }$| for all |$\epsilon>0$|⁠.

• (C.4) Gap |$\delta $| of order |$h^{2}$| and hanging nodes. Again, |$\gamma =1$| and |$\beta =1-\epsilon $| for all |$\epsilon>0$|⁠, but now |$C_{\mathrm{n} \mathrm{s}} \neq 0$| and |$ C_{\mathrm{n} \mathrm{d}} \neq 0$|⁠. Therefore, an order of convergence of |$h^{k+1/2-\epsilon }$| for all |$\epsilon>0$| is attained for the variables |$\mathbf{u}_{\mathrm{d}}$|⁠, |$\mathrm{L}_{\mathrm{s}}$| and |$p_{\mathrm{s}}$| and |$h^{k+1/2}$| for all the other variables.

• (C.5) Gap |$\delta $| of order |$h^{7/4}$| and no hanging nodes. In this case |$\gamma =3/4$|⁠, |$\beta =1/2-\epsilon $| for all |$\epsilon>0$|⁠, and |$C_{\mathrm{n} \mathrm{s}}= C_{\mathrm{n} \mathrm{d}}=0$|⁠. Then, the order of convergence is |$h^{k+5/8-\epsilon }$| for all |$\epsilon>0$| for the variables |$\mathbf{u}_{\mathrm{d}}$|⁠, |$\mathrm{L}_{\mathrm{s}}$| and |$p_{\mathrm{s}}$| and |$h^{k+7/8}$| for the rest of the variables.

• (C.6) Gap |$\delta $| of order |$h^{7/4}$| and hanging nodes. Here |$\gamma =3/4$| and |$\beta =1/2-\epsilon $| for all |$\epsilon>0$|⁠, whereas |$C_{\mathrm{n} \mathrm{s}} \neq 0$| and |$ C_{\mathrm{n} \mathrm{d}} \neq 0$|⁠. Hence, an order of convergence of |$h^{k+1/4-\epsilon }$| for all |$\epsilon>0$| is attained for the variables |$\mathbf{u}_{\mathrm{d}}$|⁠, |$\mathrm{L}_{\mathrm{s}}$| and |$p_{\mathrm{s}}$| and |$h^{k+1/2}$| for the rest of the variables.

We observe that the introduction of the constant |$\beta $| as exponent in the first term of the right-hand side of the first equation of Corollary 2 slightly improves the theoretical convergence rate of the variables involved there, despite the presence of the nonconformity constants.

We end this section by considering the particular case where the discrete interfaces |${\mathcal{I}}_{h}^{\mathrm{d}}$| and |${\mathcal{I}}_{h}^{\mathrm{s}}$| satisfy the requirement set out in Section 4.3. In this case, Lemma 7 suggests an improvement of |$h^{3/2}$| in the term associated to |$\mathbf{j}_{\mathrm{s}}^{\mathrm{n} c}$|⁠, whereas the term associated to |$j_{\mathrm{d}}^{\mathrm{n} c}$| vanishes. Thus, the semi-aligned variant of Corollary 2 is established as follows.

Let us comment on the consequences of this result.

• (D.1) No-gap and hanging nodes. Here, |$C_{\mathrm{g}}=C_{\mathrm{n}\mathrm{s}}= 0$|⁠, but |$C_{\mathrm{n} \mathrm{d}}\neq 0$|⁠. Therefore, optimal order of convergence of |$h^{k+1}$| for all the variables and order |$h^{k+2}$| for the projection of the errors |$\varepsilon ^{\mathbf{u}_{\mathrm{s}}}$| and |$\varepsilon ^{p_{\mathrm{d}}}$|⁠, which improves the power result in (C.2).

• (D.2) Gap |$\delta =h^{2}$| and hanging nodes. Here, |$\gamma =1$| and |$\beta =1-\epsilon $| for all |$\epsilon>0$|⁠, but now |$C_{\mathrm{n} \mathrm{s}} = 0$| and |$ C_{\mathrm{n} \mathrm{d}} \neq 0$|⁠. Therefore, we improve the order of convergence stated in (C.4) since now |$h^{k+1-\epsilon }$| for all |$\epsilon>0$| is attained for the projection of the errors |$\varepsilon ^{\mathbf{u}_{\mathrm{s}}}$| and |$\varepsilon ^{p_{\mathrm{d}}}$| and also for the errors in |$\mathbf{u}_{\mathrm{s}}$| and |$p_{\mathrm{d}}$|⁠.

• (D.3) Gap |$\delta =h^{7/4}$| and hanging nodes. In this case, there is no improvement in the order of convergence compared to case (C.6).

6. Numerical results

We consider four numerical experiments with the aim of illustrating the convergence of our HDG method presented in Section 3 for the two-dimensional case. In all of them, we consider the computational domain |$\varOmega = \varOmega _{\mathrm{s}} \cup \varOmega _{\mathrm{d}} \cup \varSigma $|⁠, where |$\varOmega _{\mathrm{s}} = (0,1) \times (1/2, 1)$|⁠, |$\varOmega _{\mathrm{d}} = (0, 1) \times (0, 1/2)$|⁠, and |$\varSigma = (0,1)\times \{1/2\}$|⁠. In turn, we approach our numerical examples by two computational subdomains |$\varOmega ^{\mathrm{s}}_{h} = (0,1) \times (1/2+\delta , 1)$|⁠, |$\varOmega ^{\mathrm{d}}_{h} = (0, 1) \times (0, 1/2-\delta )$|⁠, i.e., two rectangular meshed subdomains separated by a flat interface centered at |$y=1/2.$| In addition, we define the manufactured exact solution:

Also, hereafter we take |$\kappa = \mathbb{I}$|⁠, |$\nu = 1$|⁠, |$\omega _{1} = 1$| and the stabilization parameter |$\tau \equiv 1$|⁠. Subsequently, we define the errors:

where |$\boldsymbol{E}$| denotes the local extrapolation defined in (2.2) and |$p_{\mathrm{s},h}$| is the discrete pressure of our HDG scheme that satisfies (3.1e). Next, the experimental convergence rates are set as

where |$e$| and |$\widetilde{e}$| denote errors computed on two consecutive meshes with |$N$| and |$\widetilde{N}$| elements, respectively.

6.1 No gap

In our first numerical experiment, we take |$\delta = 0$| and consider two different scenarios: one free of hanging nodes and one containing hanging nodes on the discrete interfaces. For the first scenario (case (C.1) above), the results in Table 1 confirm the theoretical rate of convergence for all variables provided by Corollary 2 (for |$k = 4$|⁠, the error |$\widehat{e}$| calculated in the last two rows is affected by round-off errors). For the second scenario, we take coarser meshes for |$\varOmega _{\mathrm{d}}$| such that each interface edge corresponds to two interface edges on the meshes for |$\varOmega _{\mathrm{s}}$|⁠, effectively introducing one hanging node per side. This corresponds to case (D.1). Table 2 shows optimal order of convergence for all variables in this case, observing superconvergence |$h^{k+2}$| in the numerical trace, which is the theoretical order of convergence provided by Corollary 3.

6.2 Gap of order |$h^{7/4}$|

According to Corollary 2, |$\gamma $| must be larger than |$3/4$|⁠. In this experiment, we want to observe the behavior of the errors when |$\gamma $| is equal to |$3/4$|⁠, which means |$\delta $| of order |$h^{7/4}$|⁠. Similarly to the previous example, we divide this experiment in two different scenarios. First, we suppose again that the meshes are free of hanging nodes (case (C.5)), i.e., there is a one-to-one face bijection between the two interfaces, which means that |$C_{\mathrm{n}\mathrm{s}}=C_{\mathrm{n} \mathrm{d}}=0$|⁠. The behavior of the errors reported in Table 3 is better than the prediction of Corollary 2. In fact, we observe optimal rates for all the variables. The orders |$\widehat{e}$| are not quite clear, but they are at least |$k+1$| and the average (considering the five meshes) is |$k+2$|⁠. For the second scenario, we add hanging nodes on the bottom mesh as before, introducing one per interface edge. This corresponds to case (D.3). In Table 4 we show the results for this case. We again observe optimal and superconvergent rates, which is better than what the theory predicted.

6.3 Gap of order |$h^{}$|

Now we consider two numerical examples taking a gap |$\delta =h$|⁠. Since |$\gamma =0$|⁠, this case is not covered by the analysis. First we present a case where the meshes are free of hanging nodes. According to Table 5, the experimental rates of convergence seem to oscillate around |$k+1$| for |$e(\mathrm{L})$|⁠, |$e(\mathbf{u})$| and |$e(p)$|⁠, while for |$\widehat{e}$| the orders are at least |$k+1$|⁠. We also consider hanging nodes at the interfaces as Fig. 3 shows, noting that the meshes used are not semi-aligned for this experiment. Similarly to the previous case, according to Table 6, the errors are oscillatory but they all seem to be of order |$k+1$|⁠.