Coupling Staggered-Grid and MPFA Finite Volume Methods for Free Flow/Porous-Medium Flow Problems

TL;DR

This paper introduces a coupled discretization method combining staggered-grid and MPFA finite volume techniques to accurately simulate free flow and anisotropic porous medium flow, ensuring stability and good agreement with existing methods.

Contribution

It presents a novel coupling scheme that integrates staggered-grid and MPFA finite volume methods for improved modeling of free flow and anisotropic porous media.

Findings

Stable and applicable in anisotropic conditions

Enhanced accuracy near porous medium boundaries

Good agreement with existing models

Abstract

A discretization is proposed for models coupling free flow with anisotropic porous medium flow. Our approach employs a staggered grid finite volume method for the Navier-Stokes equations in the free flow subdomain and a MPFA finite volume method to solve Darcy flow in the porous medium. After appropriate spatial refinement in the free flow domain, the degrees of freedom are conveniently located to allow for a natural coupling of the two discretization schemes. In turn, we automatically obtain a more accurate description of the flow field surrounding the porous medium. Numerical experiments highlight the stability and applicability of the scheme in the presence of anisotropy and show good agreement with existing methods, verifying our approach.

Figures

Equations

$$\frac{\partial \varrho }{\partial t}+\nabla \cdot (\varrho \mathbf{v})$$

$$\frac{\partial (\varrho \mathbf{v})}{\partial t}+\nabla \cdot \left(\varrho \mathbf{v}{\mathbf{v}}^{\mathrm{T}}-\mu (\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^{\mathrm{T}})+p\mathbf{I}\right)$$

$$\mathbf{v}$$

$$p$$

$$\frac{\partial \varrho }{\partial t}+\nabla \cdot \left(\varrho \mathbf{v}\right)$$

$$\mathbf{v}+\frac{1}{\mu }\mathbf{K}\left(\nabla p-\varrho \mathbf{g}\right)$$

$$\mathbf{v}\cdot \mathbf{n}$$

$$p$$

$${\mathbf{v}}^{\mathrm{ff}}\cdot \mathbf{n}$$

$${\left(\varrho \mathbf{v}{\mathbf{v}}^{\mathrm{T}}-\mu (\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^{\mathrm{T}})+p\mathbf{I}\right)}^{\mathrm{ff}}\mathbf{n}$$

$${\left(-\frac{\sqrt{\mathbf{t}\cdot \mathbf{Kt}}}{{\alpha }_{\mathrm{BF}}}(\nabla \mathbf{v})\mathbf{n}-\mathbf{v}\right)}^{\mathrm{ff}}\cdot \mathbf{t}$$

$${\int }_{K^{\ast }}\frac{\partial (\varrho \mathbf{v})}{\partial t}\mathrm{d}x+{\int }_{\partial K^{\ast }}(\varrho \mathbf{v}{\mathbf{v}}^{\mathrm{T}})\cdot \mathbf{n}\mathrm{d}\Gamma$$

$$\{\mathbf{v}\}{|}_{\mathcal{P}}$$

$$[[\mathbf{v}]]{|}_{\mathcal{P}}$$

$${\int }_{K_i^{\ast }}\frac{\partial (\varrho v_i)}{\partial t}\mathrm{d}x+{\int }_{\partial K_i^{\ast }}(\{\mathbf{v}\}(\varrho v_i{)}^{\mathrm{up}})\cdot \mathbf{n}\mathrm{d}\Gamma -{\int }_{\partial K_i^{\ast }}{\mu }^{\mathrm{avg}}[[\mathbf{v}]]\cdot \mathbf{n}\mathrm{d}\Gamma$$

$${\int }_{\partial K^{\ast }}f\mathrm{d}\Gamma$$

$${\int }_{\partial K^{\ast }}\mathbf{f}\cdot \mathbf{n}\mathrm{d}\Gamma$$

$${\int }_K\frac{\partial \varrho }{\partial t}\mathrm{d}x+{\int }_{\partial K}{\varrho }^{\mathrm{up}}\mathbf{v}\cdot \mathbf{n}\mathrm{d}\Gamma ={\int }_{K}q\mathrm{d}x.$$

$${\int }_K\frac{\partial \varrho }{\partial t}\mathrm{d}x+\sum _{\sigma \in {\mathcal{E}}_K}{\int }_{\sigma }\frac{{\varrho }^{\mathrm{up}}}{{\mu }^{\mathrm{up}}}\left(-\mathbf{K}\left(\nabla p-\varrho \mathbf{g}\right)\right)\cdot \mathbf{n}\mathrm{d}\Gamma ={\int }_{K}q\mathrm{d}x.$$

$${\int }_K\frac{\partial \varrho }{\partial t}\mathrm{d}x+\sum _{\sigma \in {\mathcal{E}}_K}\frac{{\varrho }^{\mathrm{up}}}{{\mu }^{\mathrm{up}}}{F}_{K,\sigma }=Q_K,\forall K\in \mathcal{T},$$

$$F_{K,{\sigma }^v}+F_{L,{\sigma }^v}=0.$$

$${\nabla }_{\mathcal{D}}^{K^v}p\cdot ({\mathbf{x}}_{{\sigma }_i^v}-{\mathbf{x}}_K)=p_{{\sigma }_i^v}-p_K.$$

$${\nabla }_{\mathcal{D}}^{K^v}p={\mathbb{D}}_{K^v}^{-T}\left[\begin{array}{c}{p}_{{\sigma }_1^v}-p_K\\ p_{{\sigma }_3^v}-p_K\end{array}\right],\text{ with }{\mathbb{D}}_{K^v}:=\left[\begin{array}{cc}{\mathbf{x}}_{{\sigma }_1^v}-{\mathbf{x}}_K& {\mathbf{x}}_{{\sigma }_3^v}-{\mathbf{x}}_K\end{array}\right].$$

$$F_{K,{\sigma }^v}:=-|{\sigma }^v|{\mathbf{n}}_{K,{\sigma }^v}^T{\mathbf{K}}_K{\nabla }_{\mathcal{D}}^{K^v}p+{\gamma }_{K,{\sigma }^v},$$

$$\mathsf{A}{\mathsf{p}}_{\sigma }=\mathsf{B}{\mathsf{p}}_{\mathsf{K}}+\Delta \mathsf{g}.$$

$$\mathsf{f}=\mathsf{C}{\mathsf{p}}_{\sigma }+\mathsf{D}{\mathsf{p}}_{\mathsf{K}}+\mathsf{g}.$$

$$\mathsf{f}=\underset{=:\mathsf{T}}{\underbrace{\left(\mathsf{C}{\mathsf{A}}^{-1}\mathsf{B}+\mathsf{D}\right)}}{\mathsf{p}}_{\mathsf{K}}+\mathsf{C}{\mathsf{A}}^{-1}\left(\Delta \mathsf{g}\right)+\mathsf{g}.$$

$$F_{K,{\sigma }_3^v}$$

$$\mathsf{A}{\mathsf{p}}_{\sigma }=\mathsf{B}{\mathsf{p}}_{\mathsf{K}}+\mathsf{W}{\mathsf{v}}_{\mathrm{ff}}+\Delta \mathsf{g}.$$

$$\mathsf{f}=\mathsf{T}{\mathsf{p}}_{\mathsf{K}}+\mathsf{C}{\mathsf{A}}^{-1}\left(\mathsf{W}{\mathsf{v}}_{\mathrm{ff}}+\Delta \mathsf{g}\right)+\mathsf{g}$$

$$\mathbf{K}=\mathbf{R}(\alpha )\left(\begin{array}{cc}\frac{1}{\beta }k& 0\\ 0& k\end{array}\right){\mathbf{R}}^{-1}(\alpha ),\text{with}\mathbf{R}(\alpha )=\left(\begin{array}{cc}\cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{array}\right),$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Coupled systems of free flow and flow through a porous medium can be found ubiquitously in various kinds of natural and industrial contexts, including soil water evaporation [1], fuel cell water management [2], food processing [3], and evaporative cooling for turbomachinery [4]. Despite every-growing computational capacities, a discrete numerical modeling these kind of systems, including the porous geometry, is only feasible for small-scale problems. For larger domain sizes, averaging techniques involving the concept of a representative elementary volume (REV) [5], are used in order to yield upscaled models for the description of the porous domain [6]. These models can then be coupled to the free-flow region either using a single-domain or a two-domain approach. For the former, both the porous medium and the free flow are described by a single set of equations, as first introduced by Brinkman [7]. The two domains are then discerned by a spatial variation of material parameters. On the other hand, the two-domain approach decomposes the problem into two disjoint subdomains. The free-flow region is then governed by the Navier-Stokes equations while Darcy’s or Forchheimer’s law is used in the porous medium subdomain [8, 9, 10, 11]. In order to maintain thermodynamic consistency, appropriate coupling conditions have to be formulated which enforce the conservation of mass, momentum and energy across the interface between the two domains [12]. The aim of this work is to couple two different discretization schemes, and we will focus on the two-domain approach since it provides this flexibility more readily. We will consider laminar single-phase flow in the following but the extension to compositional multi-phase flow in the porous domain [11] and the use of turbulence models in the free flow part [13] is possible.

Being an active area of research, various mathematical and numerical models for the coupling of (Navier-) Stokes and Darcy/Forchheimer have been developed during recent years. Examples using the same spatial discretization scheme for both domains range from the staggered-grid finite volume method [14, 15], finite element method [16] or the box-scheme, a vertex-centered finite volume method [11, 17]. Furthermore, combinations of colocated and staggered-grid schemes have been developed [18, 19, 13]. For more discretization approaches to the modeling of coupled free flow and porous medium flow, we refer the reader to [16] and references therein.

In this work, the free flow equations are discretized using the staggered-grid method, which forms a stable numerical method for such problems. In turn, no additional stabilization techniques are necessary, in contrast to colocated schemes [20], for example. In the porous medium, we employ a multi-point-flux-approximation method (MPFA) [21]. This method was developed to overcome the shortcomings of the classical two-point-flux approximation. In particular, the MPFA scheme does not require the grid to be $\mathbf{K}$-orthogonal, meaning that the grid cells need not be in line with the principal directions of the permeability tensor $\mathbf{K}$. This is especially important for skewed and unstructured grids or in case the principal directions of the permeability tensor are inclined, in the case of layered porous structures or faults, for example. The MPFA method has been applied previously for solving Brinkman’s equation [22, 7] for a coupled system of free flow and porous medium. The referenced works thus use the same discretization scheme for both subdomains. The novelty of this work is that we employ a staggered-grid method for the free-flow model (Navier-Stokes), and couple it with a MPFA method for the porous medium model (Darcy).

The paper is organized as follows. Section 2 introduces the coupled model by presenting the governing equations in the free flow and porous medium flow subdomains, respectively, as well as the coupling conditions at the interface. The discretization schemes are introduced in Section 3, including the newly proposed coupling at the interface between the two subdomains. Numerical results are presented in Section 4 with the use of two test cases. Finally, Section 5 focuses on the conclusions.

2 Governing Equations

In order to present the governing equations for the coupled model, we first introduce the assumptions on the geometry and the notational conventions. After these preliminary definitions, we continue with the equations governing free flow and porous medium flow, followed by the coupling conditions.

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$, be an open, connected Lipschitz domain with boundary $\partial\Omega$ and $d$-dimensional measure $|{\Omega}|$. Furthermore, let $\Omega^{\mathrm{pm}}$ and $\Omega^{\mathrm{ff}}$ be a disjoint partition of $\Omega$ representing the porous medium and free flow subdomains, respectively. The subdomain boundaries are given by the interface $\Gamma^{\mathrm{if}}:=\partial\Omega^{\mathrm{ff}}\cap\partial\Omega^{\mathrm{pm}}$ as well as the remainders $\Gamma^{\mathrm{pm}}:=\partial\Omega^{\mathrm{pm}}\setminus\Gamma^{\mathrm{if}}$ and $\Gamma^{\mathrm{ff}}:=\partial\Omega^{\mathrm{ff}}\setminus\Gamma^{\mathrm{if}}$. For brevity, the superscripts are often omitted when no ambiguity arises.

The external boundary $\partial\Omega=\Gamma^{\mathrm{ff}}\cup\Gamma^{\mathrm{pm}}$ is further decomposed such that $\Gamma^{\mathrm{ff}}=\Gamma_{\mathrm{v}}^{\mathrm{ff}}\cup\Gamma_{\mathrm{p}}^{\mathrm{ff}}$ and $\Gamma^{\mathrm{pm}}=\Gamma_{\mathrm{v}}^{\mathrm{pm}}\cup\Gamma_{\mathrm{p}}^{\mathrm{pm}}$ disjointly. Here, the subscript denotes whether the velocity or the pressure is prescribed as a boundary condition. To ensure unique solvability of the resulting system, we assume that $|{\Gamma_{\mathrm{p}}^{\mathrm{pm}}\cup\Gamma_{\mathrm{p}}^{\mathrm{ff}}}|>0$, i.e. that a pressure boundary condition is imposed on a subset of the boundary $\partial\Omega$ with positive measure.

Let $\mathbf{n}$ denote the unit normal vector on $\partial\Omega$ oriented outward with respect to $\Omega$. We abuse notation and let $\mathbf{n}$ moreover denote the unit normal vector on $\Gamma^{\mathrm{if}}$ oriented outward with respect to $\Omega^{\mathrm{ff}}$.

2.1 Free Flow

In our model, the Navier-Stokes equations govern the free flow in subdomain $\Omega^{\mathrm{ff}}$. These equations are given by:

[TABLE]

The unknown variables are the velocity $\mathbf{v}$ and the pressure $p$. Here, $\varrho$ and $\mu$ denote the potentially pressure-dependent density and viscosity, respectively, while $q$ is a source (or sink) term, and $\mathbf{g}$ describes the influence of gravity. $\mathbf{I}$ is the identity tensor in $\mathbb{R}^{d\times d}$. The boundary conditions are given by known quantities $\mathbf{v}_{\Gamma}$ and $p_{\Gamma}$, representing the velocity or pressure at the corresponding boundary.

2.2 Porous Medium Flow

The equations governing single-phase flow in the porous-medium $\Omega^{\mathrm{pm}}$ are given by

[TABLE]

Equation (2b) states that the momentum balance in the porous medium is given by Darcy’s law, i.e. the Darcy velocity is calculated as $\mathbf{v}=-\frac{1}{\mu}\mathbf{K}\left(\nabla p-\varrho\mathbf{g}\right)$, with $\mathbf{K}$ being the permeability tensor. Similar to the free-flow equations (1), $q$ denotes a source or sink term. Finally $v_{\Gamma}$ and $p_{\Gamma}$ are known quantities representing the normal flux and pressure on the corresponding boundaries, respectively.

2.3 Coupling Conditions

In order to derive a thermodynamically consistent formulation of the coupled problem, conservation of mass and momentum has to be guaranteed at the interface between the porous medium and the free-flow domain. We therefore impose the following interface conditions:

[TABLE]

The momentum transfer normal to the interface is given by (3b) [9]. Condition (3c) is the commonly used Beavers-Joseph-Saffman slip condition [23, 24]. Here, $\mathbf{t}$ denotes any unit vector from the tangent bundle of $\Gamma^{\mathrm{if}}$ and $\alpha_{\mathrm{BF}}$ is a parameter to be determined experimentally. We remark that this condition is technically a boundary condition for the free flow, not a coupling condition between the two flow regimes. Furthermore, it has been developed for free flow strictly parallel to the interface and might lose its validity for other flow configurations.

3 Discretization

This section is devoted to giving an outline of the numerical schemes used in the individual subdomains and the incorporation of the interface conditions (3). However, we first introduce some notational conventions concerning the partition of $\Omega$ in the following definition.

Definition 1** (Grid discretization).**

The tuple $\mathcal{D}:=(\mathcal{T},\mathcal{E},\mathcal{P},\mathcal{V})$ denotes the grid discretization, in which

1. (i)

$\mathcal{T}$* is the set of control volumes (cells) such that $\overline{\Omega}=\cup_{K\in\mathcal{T}}\overline{K}$. For each cell $K\in\mathcal{T}$, $|{K}|>0$ denotes the cell volume.* 2. (ii)

$\mathcal{E}$* is the set of faces such that each face $\sigma$ is a $(d-1)$-dimensional hyperplane with measure $|{\sigma}|>0$. For each cell $K\in\mathcal{T}$, $\mathcal{E}_{K}$ is the subset of $\mathcal{E}$ such that $\partial K=\cup_{\sigma\in\mathcal{E}_{K}}{\sigma}$. Furthermore, $\mathbf{x}_{\sigma}$ denotes the face evaluation points and $\mathbf{n}_{K,\sigma}$ the unit vector that is normal to $\sigma$ and outward to $K$.* 3. (iii)

$\mathcal{P}:=\{\mathbf{x}_{K}\}_{K\in\mathcal{T}}$* is the set of cell centers (not required to be the barycenters) such that $\mathbf{x}_{K}\in K$ and $K$ is star-shaped with respect to $\mathbf{x}_{K}$. For all $K\in\mathcal{T}$ and $\sigma\in\mathcal{E}_{K}$, let $d_{K,\sigma}$ denote the Euclidean distance between $\mathbf{x}_{K}$ and $\sigma$.* 4. (iv)

$\mathcal{V}$* is the set of vertices of the grid, corresponding to the corners of the cells.*

For ease of exposition, we assume $d=2$, however, this is not a limitation and the model can be readily extended for three dimensions.

3.1 Staggered grid scheme

A staggered-grid finite volume scheme, also known as MAC scheme [15], is used in the free-flow subdomain. Here, scalar quantities including pressure and density are stored on the cell centers $\mathcal{P}$ while the velocity degrees of freedom are located on the primary control volumes’ faces $\mathcal{E}$. The resulting scheme is stable, hence oscillation-free solutions are guaranteed without the need for additional stabilization techniques. This is in contrast with colocated schemes, in which all unknowns are defined at the same location [20]. $\mathcal{D}$ is a uniform Cartesian grid with mesh size $h$ in both directions.

We start with the momentum balance equations. For each face in $\mathcal{E}$, we construct a secondary control volume $K^{*}$ with boundary $\partial K^{*}$, as depicted in Figure 1. On each $K^{*}$, we integrate the momentum balance equation and apply Gauss divergence theorem to obtain

[TABLE]

The first and second components of this vector equation are considered separately. Due to the different locations at which the degrees of freedom are defined, we require interpolation operators in order to continue. Let us introduce the average ($\{\cdot\}$) and jump quantities ($\llbracket\cdot\rrbracket$) on cell centers ($\mathcal{P}$) and vertices ($\mathcal{V}$) of the primal grid for $\mathbf{v}=[v_{x},v_{y}]^{\text{T}}$ (see Figure 1):

[TABLE]

Here, the superscript $\{E,N,W,S\}$ refers to the closest degree of freedom East, North, West, or South of the evaluation point, see Figure 1.

In the following, the superscript “up” denotes the upwind quantity relative to the velocity $\mathbf{v}$. Moreover, we introduce $\mu^{\mathrm{avg}}$ such that in all cell centers $\mathcal{P}$, it denotes the corresponding viscosity whereas in each vertex $\mathcal{V}$, it is the viscosity averaged over the adjacent cells.

With the operators defined, we discretize the momentum balance equation for each secondary control volume $K_{i}^{*}$, in which the subscript $i\in\{x,y\}$ denotes whether the control volume surrounds a vertical or horizontal face, respectively. The discretized equation for component $i$ is then given by

[TABLE]

We emphasize that the boundary integrals are computed numerically using the following rules for a scalar-valued function $f$ and a vector-valued function $\mathbf{f}$:

[TABLE]

By definition, each $\mathbf{x}_{\sigma^{*}}$ will either be a cell center ($\mathcal{P}$) or a vertex ($\mathcal{V}$) of the grid.

Finally, the mass balance is evaluated on each cell of the grid, i.e. we compute for each $K\in\mathcal{T}$:

[TABLE]

The above equations fully define the staggered-grid discretization scheme for all internal faces of the grid $\mathcal{D}$. For incorporation of boundary conditions, we refer the reader to [20].

3.2 Cell-centered finite volume scheme

A cell-centered finite volume scheme is employed in the Darcy subdomain, i.e. the grid elements are used as control-volumes and the degrees of freedom are associated with the cell centers. Typically, the finite volume formulation is obtained by integrating the first equation of (2) over a control volume $K\in\mathcal{T}$ and by applying the Gauss divergence theorem:

[TABLE]

Replacing, the exact fluxes by an approximation, i.e. $F_{K,\sigma}\approx\int_{\sigma}\left(-\mathbf{K}_{K}\left(\nabla p-\varrho\mathbf{g}\right)\right)\cdot\mathbf{n}\,\mathrm{d}\Gamma$ (here $\mathbf{K}_{K}$ is the value of $\mathbf{K}$ associated with control volume $K$), yields

[TABLE]

where $F_{K,\sigma}$ is the discrete flux through face $\sigma$ flowing out of cell $K$, $Q_{K}:=\int_{K}q\,\mathrm{d}x$ is the integrated source/sink term, and $(\cdot)^{\mathrm{up}}$ denotes upwinding with respect to the sign of the flux $F_{K,\sigma}$.

Finite volume schemes primarily differ in the approximation of the term $(\mathbf{\mathbf{K}}_{K}\nabla p)\cdot\mathbf{n}$ (i.e. the choice of the fluxes $F_{K,\sigma}$). The widely used linear two-point flux approximation (TPFA), for example, is a simple but robust scheme. However, it is well-known that it is inconsistent on grids that are not $\mathbf{K}$-orthogonal (see e.g. [25]). In this work we consider anisotropic permeability tensors in the porous medium and $\mathbf{K}$-orthogonality of the grid can thus not be guaranteed. Therefore, we employ a multi-point flux approximation (MPFA) scheme for the formulation of the discrete fluxes, which has been presented in [21]. This particular scheme is termed MPFA-O and is only one among many methods that fall into the family of MPFA schemes ([26, 27]). Please note that we will omit the suffix “-O” throughout this document wherever it would affect the readability.

For the computation of the fluxes, a dual grid is created by connecting the barycenters of the cells with the barycenters of the faces ($d=2$) or the barycenters of the faces and edges ($d=3$). This divides each cell into sub-control volumes $K^{v}$ with $K\in\mathcal{T}$ and $v\in\mathcal{V}$. Analogously, each face is sub-divided into sub-control volume faces $\sigma^{v}$, see Figure 2. Expressions for the face fluxes $F_{K,\sigma^{v}}$ are obtained by introducing the face pressures $p_{\sigma^{v}}$. The location of these face pressures along the sub-control volume faces $\sigma^{v}$ is parameterized by $\xi$, $0\leq\xi<1$, and is the center of the original face $\sigma$ for $\xi=0$ and would be the position of the vertex $v$ for $\xi=1$. These face pressures are then eliminated by enforcing the continuity of fluxes across each sub-control volume face. I.e., for each face $\sigma^{v}$ between $K^{v}$ and $L^{v}$, we impose:

[TABLE]

We allow for piecewise constant $\mathbf{K}$ (denoted as $\mathbf{K}_{K}$ for each cell $K$) and construct discrete gradients $\nabla_{\mathcal{D}}^{K^{v}}p$, per sub-control volume $K^{v}$, depending on its two embedded sub-control volume faces. Let us consider $K^{v}$ in Figure 2 with faces $\sigma_{1}^{v}$ and $\sigma_{3}^{v}$. Here, the discrete gradients are constructed to be consistent such that the following holds for $i\in\{1,3\}$:

[TABLE]

Thus, a discrete gradient (for sub-control volume $K^{v}$) that fulfills the two conditions (12) is defined by

[TABLE]

This enables us to write the discrete flux across $\sigma^{v}$ between $K^{v}$ and $L^{v}$ as follows:

[TABLE]

where we introduced $\gamma_{K,\sigma^{v}}=\rho_{\sigma^{v}}|{\sigma^{v}}|\mathbf{n}_{K,\sigma^{v}}^{T}\mathbf{\mathbf{K}}_{K}\mathbf{g}$, with $\rho_{\sigma^{v}}=\frac{\rho_{K}+\rho_{L}}{2}$, to incorporate the effect of gravity.

To deduce a cell-centered scheme, the face pressures $p_{\sigma^{v}_{i}}$ are eliminated. This is done by enforcing flux continuity (11a) within each interaction volume and by solving a local system of equations.

We rewrite these conditions in matrix form, and introduce the sans serif font to denote the corresponding matrices and vectors. All local face pressures of an interaction region are collected in the vector $\mathsf{p_{\sigma}}$, cell pressures in the vector $\mathsf{p_{K}}$, and all terms related to gravity in the vector $\mathsf{g}$. Flux continuity then allows us to rewrite the face pressures in terms of the cell pressures:

[TABLE]

Here, the $\Delta$ represents the difference in contributions due to gravity over each face. Let $\mathsf{f}$ denote the vector of all fluxes across the sub-control volume faces of the interaction region. These can be expressed in matrix form using equation (14):

[TABLE]

With these introduced matrices the final expressions for the local sub-control volume face fluxes, related to the interaction region, read:

[TABLE]

The entries of the matrix $\mathsf{T}$ are often referred to as the transmissibilities.

3.3 Coupling

In this section, we consider the realization of the coupling conditions (3). As depicted in Figure 3, the grids are chosen to be non-matching at the interface such that each sub-control volume face coincides with a face $\sigma\in\mathcal{E}^{\mathrm{ff}}$ of the free-flow domain. In turn, a natural coupling arises between the staggered grid discretization and the MPFA method, due to the coinciding degrees of freedom (for $\xi=0.5$, see Figure 3). We emphasize that each cell in the porous-medium subdomain has two neighboring cells in the free-flow subdomain (two-dimensional setup).

We start with the flux continuity condition (3a). As depicted in Figure 3, let us consider a sub-control volume $K^{v}$ located at the interface such that $\sigma_{3}^{v}$ is located on $\Gamma^{\mathrm{if}}$. We then let the velocity from the free-flow domain determine the flux over $\sigma_{3}^{v}$:

[TABLE]

Collecting the right-hand side of (18) in the matrix-vector product $\mathsf{Wv_{\mathrm{ff}}}$, we can rewrite (15) for interaction regions that are located at the interface to the free-flow domain:

[TABLE]

In the situation shown in Figure 3, the face pressures $p^{\mathrm{pm}}_{\sigma^{v}_{1}},p^{\mathrm{pm}}_{\sigma^{v}_{2}},p^{\mathrm{pm}}_{\sigma^{v}_{3}}$ are thus dependent on the primary unknowns $p^{\mathrm{pm}}_{K},p^{\mathrm{pm}}_{L}$ of the porous-medium domain and the face velocities $v^{\mathrm{ff}}_{\sigma^{v}_{2}},v^{\mathrm{ff}}_{\sigma^{v}_{3}}$ of the free-flow domain (i.e. $p^{\mathrm{pm}}_{\sigma^{v}_{i}}=p^{\mathrm{pm}}_{\sigma^{v}_{i}}(p^{\mathrm{pm}}_{K},p^{\mathrm{pm}}_{L},v^{\mathrm{ff}}_{\sigma^{v}_{2}},v^{\mathrm{ff}}_{\sigma^{v}_{3}})$). Insertion of (19) in (16) leads to the expression

[TABLE]

for the fluxes across sub-control volume faces within interaction regions located at the interface to the free-flow domain.

The remaining coupling conditions are imposed as followed. The momentum balance (3b) is enforced using the reconstructed face pressures from (19). On the other hand, the Beavers-Joseph-Saffman condition (3c) is technically a boundary condition for the free-flow problem, as previously noted in Section 2, and is implemented accordingly.

Finally, we remark that in the case of compressible fluids, the density and viscosity are pressure-dependent. The upwind terms $\mu^{\mathrm{pm},\mathrm{up}},\varrho^{\mathrm{pm},\mathrm{up}},\varrho^{\mathrm{ff},\mathrm{up}}$ are then evaluated using the cell-pressure unknowns, i.e. for a face $\sigma$ between porous-medium cell $K$ and free-flow cell $L$, these terms therefore depend on the pressures $p^{\mathrm{pm}}_{K}$ and $p^{\mathrm{ff}}_{L}$.

4 Numerical results

All simulations are performed using the open-source simulator DuMu${}^{\text{x}}$ [28], which comes in the form of an additional DUNE module [29]. We employ a monolithic approach, where both sub-problems are assembled into one system of equations and use an implicit Euler method for the time discretization. Newton’s method is used to solve the non-linearities involved in the systems of equations. For all test cases, the compressible fluid “air” (see the DuMu${}^{\text{x}}$ documentation [28]) is used. We consider two-dimensional setups, however, the implementation is also able to handle three-dimensional domains.

4.1 Test case 1

The first test case is similar to the one that has been presented in [14]. However, the authors of [14] used a Navier-Stokes-Brinkman-type system for both domains and this system is discretized by using a staggered-grid scheme in both domains. For anisotropic permeability tensors, this requires the interpolation of all velocity components at grid faces.

In this work, we use a different approach, where a staggered-grid scheme is used to discretize the free-flow system and a cell-centered finite volume scheme to discretize the porous-medium system. Thus, no additional velocity degrees of freedom are needed for the porous-medium domain. However, for anisotropic permeability tensors or unstructured grids this requires more sophisticated cell-centered finite volume schemes, as for example the MPFA scheme that has been presented in Section 3.2.

Air is flowing through a two-dimensional channel which is partially blocked by a rectangular porous medium as shown in Figure 4. The first test case involves small Reynolds numbers ($Re\ll 1$) with respect to the average velocity in the narrow section in the channel above the porous medium. In this case, only the stationary solution (which is reached after a few time steps, starting from a resting fluid) is investigated. The top and bottom of the domain are considered as rigid, impermeable walls with $\mathbf{v}=0$ (including the wall part below the porous box). Flow is driven by a pressure difference between the left and the right boundary which is set to $\Delta p=10^{-6}\,\mathrm{Pa}$.

The permeability tensor in the porous medium is given as

[TABLE]

where $\alpha$ is the rotation angle, $\beta\geq 1$ the anisotropy ratio, and $k=10^{-6}\,\mathrm{m^{2}}$. In the following, the influence of the anisotropy on the total mass fluxes crossing the boundaries $\Gamma^{\mathrm{pm}}_{\mathrm{in}},\Gamma^{\mathrm{pm}}_{\mathrm{out}},\Gamma^{\mathrm{pm}}_{\mathrm{top}}$ is investigated for $\alpha\in\{-45^{\circ},-30^{\circ},0^{\circ},30^{\circ},45^{\circ}\}$ and $\beta\in\{10,100\}$. In the free-flow domain the same grid is used for the MPFA and the TPFA schemes, whereas in the porous medium the grid used for the MPFA scheme is coarsened such that we obtain a non-matching interface, as shown in Figure 3.

Figure 5 and 6 show the velocity magnitude and the pressure profile of the MPFA and TPFA schemes for $\beta=100$ and $\alpha\in\{-45^{\circ},45^{\circ}\}$. For the TPFA scheme, the results are only shown for $\alpha=45^{\circ}$ because the results for $\alpha=-45^{\circ}$ are identical. Due to the flow resistance imposed by the porous medium, most of the air passes this obstacle through the constricted section above the block, thus leading to the highest flow velocities there. While the gas passes the porous block virtually parallel for both $\alpha=45^{\circ}$ and $\alpha=0^{\circ}$ when applying the TPFA method, the effect of anisotropy is clearly visible in the MPFA results. Here, the flow follows the inclined principal direction of the permeability tensor, exiting or entering the porous domain at the top. Small regions of local recirculation can be found within the obstacle which are caused by the medium’s anisotropy and the closed wall at the bottom. For $\alpha=-45^{\circ}$, the upward flow in the box creates a recirculation at the right part of the bottom, where a small amount of gas is actually pulled from the free-flow channel into the porous medium. Analogously for $\alpha=45^{\circ}$, the downward flow causes a recirculation at the bottom left, where the gas cannot exit through the solid wall and thus leaves the domain towards the left, in opposition to the general flow field.

Due to the small pressure differences, we subtract a reference value of $p_{\mathrm{ref}}=10^{5}\,\mathrm{Pa}$ for improved visibility as shown on the right side of Figures 5 and 6. Almost the entire pressure drop along the channel is observed at the porous domain. Again, the influence of anisotropy is clearly visibly for the MPFA results in terms of an inclined pressure profile which also reflects the closed bottom of the domain. None of these effects can be captured by the TPFA method.

Figure 7 shows the total mass fluxes over the porous-medium free-flow interfaces. For $\alpha=0$, the TPFA and MPFA scheme result in almost the same solutions (small differences occur because a finer mesh is used for the TPFA scheme in the porous medium). For this angle, the TPFA scheme is also consistent and produces the correct results. However, it is observed that for all other angles the total mass fluxes significantly differ. The off-diagonal terms of $\mathbf{K}$ are not considered in the TPFA transmissibilities, which explains why the TPFA results are independent from the direction of rotation. Furthermore, the total fluxes over $\Gamma^{\mathrm{pm}}_{\mathrm{top}}$ are small for the TPFA scheme, in contrast to the MPFA scheme, where the total mass fluxes increase with increasing rotation angle due to the contribution from the non-parallel porous-medium flow as mentioned above.

4.2 Test case 2

The next test case investigates the solution behavior for higher Reynolds numbers. It uses a similar setting as the previous one with the difference that the channel is elongated in $x$-direction. The computational domains are given by $\Omega=[0,2.5]\times[0,0.25]$ m and $\Omega^{\mathrm{pm}}=[0.4,0.6]\times[0,0.2]$ m such that $\Omega^{\mathrm{ff}}=\Omega\setminus\Omega^{\mathrm{pm}}$. A pressure difference of $\Delta p=2\cdot 10^{-3}\,\mathrm{Pa}$ between the left and the right boundary results in $Re\approx 130$ in the channel right atop the porous block.

For this test case, it takes much longer until a stationary solution is reached, which is why the solutions are investigated at different time steps. In the following, we focus on the discussion of the case where $\alpha=0^{\circ}$ or $\alpha=45^{\circ}$ and $\beta=100$. Figure 8 shows the resulting velocity fields of the MPFA method at different times for $\alpha=45^{\circ}$. As before, the gas hits the front of the porous block and is forced mainly through the narrow channel section above it which acts as some sort of duct through which a jet of high velocity fluid streams in the open channel on the right side. As the velocity gradually increases over time, the formation of vortex structures within the free-flow channel can be observed after around $20\,\mathrm{s}$. Having reached equilibrium at around $1000\,\mathrm{s}$, the system features two stable countercurrent, larger vortices downstream the porous obstacle and one small recirculation zone in front of the block. Again, the flow within the obstacle is clearly influenced by the anisotropy, a feature not reproduced by the TPFA method as seen in Figure 9. Here, the flow within the block does not follow $\mathbf{K}$’s orientation. Instead, the fluid immediately strives for the top of the obstacle once it has entered it which is due to the strongly increased vertical permeability ($\beta=100$). Again, the TPFA scheme is not able to capture this relevant anisotropy effect as it only considers the main diagonals of $\mathbf{K}$. Figure 10 depicts the cumulative mass fluxes per unit depth across the porous medium’s boundaries (top row) and across the plane at the center of the narrow channel section ($x=0.5\,\mathrm{m}$, bottom row) for a case with $\alpha=0^{\circ}$ (left column) and one with $\alpha=45^{\circ}$ (right column). For the former, the results of the TPFA and MPFA are very similar. With increasing time, the total mass flux entering the obstacle’s front (red line) increases in accordance with the global velocity field (see Figure 8) until it approaches a constant value after approximately $100\,\mathrm{s}$. The same applies for the fluxes leaving the obstacle’s top (blue line), which approaches the same value as the incoming fluxes but with a different sign. This indicates that no mass flux occurs over the obstacle’s back which can also be seen from the black line. This line (and thus, the fluxes over the obstacle’s top) only deviates from zero at the very beginning of the simulation where it has to balance out the initial disequilibrium between the red and blue curve. Due to the imposed anisotropy ratio of $\beta=100$, the horizontal permeability is 100 times lower than the vertical one and fluid is immediately pushed towards the top, where it exits the domain again. Even though the fluid’s density is actually pressure-dependent, significant compressibility effects could not be observed for the given setup. The differences between TPFA and MPFA with respect to the blue and red curves (front and back of the obstacle) are most likely due to the different discretization width in the porous domain as explained before. In total, less gas seems to enter the porous block for the MPFA method. The plot on the lower left of Figure 10 shows the temporal evolution of the total mass flux through the channel above the obstacle. Both the TPFA and the MPFA scheme converge to the same result, the grid resolution of the free-flow domain is identical for both cases. Note that even with a constant flux through the constriction after $100\,\mathrm{s}$, the global velocity field downstream the obstacle still changes, including the formation of vortices as described before.

Considering the right column of Figure 10 ($\alpha=45^{\circ}$), the differences between the two schemes are significantly higher. For the MPFA method, we observe a considerable inflow across the obstacle’s top boundary, coming from the constricted free-flow channel and drawn by the downwards inclined flow field within the porous medium. At the same time, there is only a very limited inflow through the obstacle’s front which is due the obstacle’s anistropy.

This is in strong contrast to the TPFA method’s result, where air is still leaving the obstacle’s top as also described before and seen in Figure 9. The MPFA method results in a higher flux through the constricted channel (see bottom right of Figure 10) as there is less fluid entering the porous obstacle as explained above.

4.3 Test case 3

Another advantage of using MPFA in the porous medium domain is the ability to use unstructured grids while maintaining consistency of the scheme. The test case presented in this section considers a porous medium domain in which geometrical constraints favor the use of triangles (unstructured) over quadrilaterals for its discretization. This situation can arise, for example, in environmental applications considering the exchange processes between the atmosphere and the subsurface, where the latter can be composed of complex shaped geological layers. An illustration of the setup of this test case is given in Figure 11, while a detailed view on the discretization of the two compartments is depicted in Figure 12.

Two sets of simulations were performed using the rotation angles $\alpha=45^{\circ}$ and $\alpha=-45^{\circ}$, together with $\beta=10$, $k=10^{-6}\,\mathrm{m^{2}}$ and $\Delta p=1\cdot 10^{-6}\,\mathrm{Pa}$. The results using both MPFA and TPFA in the porous medium domain are shown in the Figures 13 and 14, respectively. Please note that unlike before, the velocity vectors are now scaled by magnitude. As expected, the solutions using TPFA exhibit a more distorted velocity field originating from the scheme being inconsistent on both unstructured grids and for anisotropic tensors. On the other hand, the solutions using MPFA provide velocity fields that follow the geometrical features of the porous medium. For $\alpha=-45^{\circ}$, two main regions in which the flux from the porous medium into the free-flow domain is concentrated, and where the maximum velocities occur, can be observed (see Figure 13). The direction of these maximum velocities coincides with the direction of highest permeability. In contrast to that, the inflow from the free-flow domain into the porous medium occurs in the direction of the lowest permeability and thus at smaller velocities. In the case of $\alpha=45^{\circ}$, the highest velocities are observed in the regions where an inflow from the free-flow domain into the porous medium occurs, again following the direction of the highest permeability (see Figure 14). With TPFA being used in the porous medium domain, the differences between the velocity fields obtained from the two angles turn out to be smaller in comparison to the solutions obtained with MPFA. Furthermore, the low velocity regions seem to be generally overestimated. This effect shows itself also in the integrated transfer flux across the interface, which results in $1.3\cdot 10^{-8}\,\mathrm{kg/s}$ for MPFA and $2.2\cdot 10^{-8}\,\mathrm{kg/s}$ for TPFA, thus, around $70\,\mathrm{\%}$ higher.

5 Conclusions

In this work, a discretization method is proposed for problems concerning free flow coupled to porous medium flow. The method combines the stability of the staggered grid finite volume method for the free-flow equations with the consistency of MPFA finite volume methods for flows in anisotropic porous media. We have shown how appropriate alignment of the grids results in a natural coupling between the two discretization schemes due to coinciding degrees of freedom.

The stability and consistency of the method are emphasized with the use of numerical experiments. Especially in the presence of anisotropy in the porous medium, a significant difference is observed with respect to the widely used, but inconsistent, TPFA finite volume method. We moreover emphasized the use of unstructured grids in the porous medium, allowing for computations on more general geometries. The use of unstructured grids in the free-flow domain [30] is a topic for future investigation.

Future work will moreover address the extension of the presented model to multi-phase flow simulations including compositional and non-isothermal effects, with a special focus on evaporative processes at the interface between the porous medium and the free-flow domain. This extension would make the model applicable to a large variety of applications, as e.g. the drying of soil due to evaporation [11] and subsequent soil salinization [31]. This drying process can lead to the creation of fractures within the porous medium. Therefore, we want to employ a discrete fracture model for the description of the porous medium domain, for which there is a preexisting implementation available in DuMu${}^{\text{x}}$ (see [32]). To increase efficiency, we will also use other finite volume schemes for the porous domain, as for example those that have been recently presented in [33].

Acknowledgements

We thank the German Research Foundation (DFG) for supporting this work by funding SFB 1313, Research Project A02.