Validation and calibration of coupled porous-medium and free-flow problems using pore-scale resolved models

TL;DR

This study validates and calibrates coupled free-flow and porous-medium models using pore-scale simulations, highlighting the importance of interface conditions and effective parameters for accurate predictions in complex geometries.

Contribution

It introduces a validation framework for coupled Stokes-Darcy models using pore-scale resolved simulations across different geometries and interface conditions.

Findings

Coupled model results are sensitive to interface location and parameters.

Effective parameters vary with geometrical configuration.

Validation improves understanding of interface condition impacts.

Abstract

The correct choice of interface conditions and effective parameters for coupled macroscale free-flow and porous-medium models is crucial for a complete mathematical description of the problem under consideration and for accurate numerical simulation of applications. We consider single-fluid-phase systems described by the Stokes-Darcy model. Different sets of coupling conditions for this model are available. However, the choice of these conditions and effective model parameters is often arbitrary. We use large scale lattice Boltzmann simulations to validate coupling conditions by comparison of the macroscale simulations against pore-scale resolved models. We analyse two settings (lid driven cavity over a porous bed and infiltration problem) with different geometrical configurations (channelised and staggered distributions of solid grains) and different sets of interface conditions. Effective parameters for the macroscale models are computed numerically for each geometrical configuration. Numerical simulation results demonstrate the sensitivity of the coupled Stokes-Darcy problem to the location of the sharp fluid-porous interface, the effective model parameters and the interface conditions.

Figures

Equations

$$\nabla \cdot \mathrm{v}=0\text{in}{\Omega }_{\mathrm{ff}}.$$

$$-\nabla \cdot \mathsf{T}(\mathrm{v},p)-\rho \mathrm{g}=0\text{in}{\Omega }_{\mathrm{ff}},$$

$$\begin{array}{rl}\mathrm{v}& ={\overline{\mathrm{v}}}_{\mathrm{ff}}\text{on}{\Gamma }_{D,\mathrm{ff}},\\ \mathsf{T}(\mathrm{v},p)\cdot {\mathrm{n}}_{\mathrm{ff}}& ={\overline{\mathrm{h}}}_{\mathrm{ff}}\text{on}{\Gamma }_{N,\mathrm{ff}},\\ \partial {\mathrm{v}}_{\mathrm{ff}}/\partial {\mathrm{n}}_{\mathrm{ff}}& =0\text{on}{\Gamma }_{\text{out},\mathrm{ff}},\end{array}$$

$$\nabla \cdot \mathrm{v}=q\text{in}{\Omega }_{\mathrm{pm}},$$

$$\mathrm{v}=-\frac{\mathsf{K}}{\mu }\left(\nabla p-\rho \mathrm{g}\right).$$

$$\begin{array}{rl}p& ={\overline{p}}_{\mathrm{pm}}\text{on}{\Gamma }_{D,\mathrm{pm}},\\ \mathrm{v}\cdot {\mathrm{n}}_{\mathrm{pm}}& ={\overline{v}}_{\mathrm{pm}}\text{on}{\Gamma }_{N,\mathrm{pm}},\end{array}$$

$${\mathrm{v}}_{\mathrm{ff}}\cdot \mathrm{n}={\mathrm{v}}_{\mathrm{pm}}\cdot \mathrm{n}\text{on}\Gamma ,$$

$$-{\mathrm{n}}_{\mathrm{ff}}\cdot \mathsf{T}\left({\mathrm{v}}_{\mathrm{ff}},p_{\mathrm{ff}}\right)\cdot {\mathrm{n}}_{\mathrm{ff}}=p_{\mathrm{pm}}\text{on}\Gamma ,$$

$$\left({\mathrm{v}}_{\mathrm{ff}}-{\mathrm{v}}_{\mathrm{pm}}\right)\cdot {\tau }_j+\frac{\sqrt{\mathsf{K}}}{{\alpha }_{\text{BJ}}}\left(\nabla {\mathrm{v}}_{\mathrm{ff}}\cdot {\mathrm{n}}_{\mathrm{ff}}\right)\cdot {\tau }_j=0\text{on}\Gamma ,$$

$${\mathrm{v}}_{\mathrm{ff}}\cdot {\tau }_j+\frac{\sqrt{\mathsf{K}}}{{\alpha }_{\mathrm{BJ}}}\left(\nabla {\mathrm{v}}_{\mathrm{ff}}\cdot {\mathrm{n}}_{\mathrm{ff}}\right)\cdot {\tau }_j=0\text{on}\Gamma ,$$

$${\mathrm{v}}_{\mathrm{ff}}\cdot {\tau }_j+\frac{\sqrt{\mathsf{K}}}{{\alpha }_{\text{BJ}}}{\mathrm{n}}_{\mathrm{ff}}\cdot \mathsf{D}\left({\mathrm{v}}_{\mathrm{ff}}\right)\cdot {\tau }_j=0\text{on}\Gamma .$$

$$\begin{array}{rl}{v}_{\mathrm{ff},2}& =0,\\ p_{\text{ff}}+\mu C_{\omega }^{\text{bl}}\frac{\partial v_{\text{ff},1}}{\partial x_2}& =p_{\text{pm}},\\ v_{\text{ff},1}-\epsilon C_1^{\text{bl}}\frac{\partial v_{\text{ff},1}}{\partial x_2}& =0,\end{array}$$

$$\begin{array}{rl}{v}_{\mathrm{ff},2}& =v_{\mathrm{pm},2},\\ v_{\mathrm{ff},1}& =M_1^{\text{bl}}\frac{\partial p_{\mathrm{pm}}}{\partial x_2},\\ p_{\mathrm{pm}}& =0,\end{array}$$

$$\begin{array}{rl}-\mu \Delta {\text{v}}^{\epsilon }+\nabla p^{\epsilon }& =0\text{in }{\Omega }^{\epsilon },\\ \nabla \cdot {\text{v}}^{\epsilon }& =0\text{in }{\Omega }^{\epsilon },\end{array}$$

$${\mathrm{v}}^{\epsilon }=0\text{on }\partial {\Omega }^{\epsilon }\setminus \partial \Omega ,$$

$$\begin{array}{rl}{\text{v}}^{\epsilon }& =(0,0)\text{on }\partial \Omega \setminus \{x_2=H\},\\ {\text{v}}^{\epsilon }& =(1,0)\text{on }\{x_2=H\},\end{array}$$

$$\begin{array}{rl}& {\mathrm{v}}^{\epsilon }=(0,0)\text{on }\{x_1=0\}\cup \{x_1=L\},\\ & {\text{v}}^{\epsilon }=(0,-\sin (\pi x_1))\text{on }\{x_2=H\},\\ & \left(\nabla {\text{v}}^{\epsilon }-p^{\epsilon }\text{I}\right)\mathrm{n}=-p_b\mathrm{n}\text{on }\{x_2=h\},\end{array}$$

$$\begin{array}{rl}{\text{v}}^{\epsilon }(\text{x})& ={\epsilon }^2{\text{v}}_0(\text{x},\text{y})+{\epsilon }^3{\text{v}}_1(\text{x},\text{y})+\dots \\ p^{\epsilon }(\text{x})& =p_0(\text{x},\text{y})+\epsilon p_1(\text{x},\text{y})+\dots \end{array}$$

$$\begin{array}{r}\tilde{\mathsf{K}}=({\tilde{k}}_{ij}{)}_{i,j=1,2}={\int }_{Y_{\text{f}}}{\text{w}}_i^j\text{d}\text{y},\end{array}$$

$$\begin{array}{rl}-{\Delta }_{\text{y}}{\text{w}}^j+{\nabla }_{\text{y}}{\pi }^j& ={\text{e}}_j\text{in Yf},\\ {\nabla }_{\text{y}}\cdot {\text{w}}^j& =0\text{in Yf},\\ {\text{w}}^j& =0\text{on ∂Yf∖∂Y},\\ \{{\text{w}}^j,{\pi }^j\}& \text{ is 1-periodic in }\text{y}.\end{array}$$

$$\rho (\mathrm{x},t)=\sum _i{f}_i(\mathrm{x},t),\rho \mathrm{v}(\mathrm{x},t)=\sum _i{\mathrm{c}}_i{f}_i(\mathrm{x},t).$$

$$f_i(\mathrm{x}+{\mathrm{c}}_i\Delta t,t+\Delta t)=f_i(\mathrm{x},t)+{\Theta }_i(\mathrm{x},t)$$

$${\Theta }_i(\mathrm{x},t)=-\frac{f_i(\mathrm{x},t)-f_i(\mathrm{x},t{)}^{\text{eq}}}{\tau }\Delta t,$$

$$f_i^{\text{eq}}(\mathrm{x},t)=w_i\rho \left(1+\frac{\mathrm{v}\cdot {\mathrm{c}}_i}{c_{\text{s}}^2}+\frac{(\mathrm{v}\cdot {\mathrm{c}}_i{)}^2}{2c_{\text{s}}^4}-\frac{\mathrm{v}\cdot \mathrm{v}}{2c_{\text{s}}^2}\right)$$

$$\nu =c_{\text{s}}^2\left(\tau -\frac{\Delta t}{2}\right).$$

$$Re=\frac{v_{\text{lid}}\mathcal{L}}{\nu }$$

$$\tilde{\mathsf{K}}=\left(\begin{array}{cc}5.671\cdot 10^{-4}& 0\\ 0& 5.671\cdot 10^{-4}\end{array}\right).$$

$$\begin{array}{rl}& {\mathrm{v}}_{\mathrm{ff}}=(1,0)\text{on }\{x_2=1\},\\ & {\mathrm{v}}_{\mathrm{ff}}=(0,0)\text{on }(\{0\}\cup \{1\})\times (0,1),\\ & {\mathrm{v}}_{\mathrm{pm}}\cdot {\mathrm{n}}_{\mathrm{pm}}=0\text{on }\{x_2=-0.5\}\\ & {\mathrm{v}}_{\mathrm{pm}}\cdot {\mathrm{n}}_{\mathrm{pm}}=0\text{on }(\{0\}\cup \{1\})\times (-0.5,0).\end{array}$$

$$C_1^{\text{bl}}=-5.41\cdot 10^{-2},C_{\omega }^{\text{bl}}=0.$$

$$\tilde{\mathsf{K}}=\left(\begin{array}{cc}1.325\cdot 10^{-4}& 0\\ 0& 1.325\cdot 10^{-4}\end{array}\right),$$

Full text

1 Introduction

Coupled flow systems containing a porous-medium domain and a free-flow region appear in various environmental and technical applications, such as surface-water/ groundwater flow, industrial filtration and water management in fuel cells. The interaction between the flow regions is dominated by the interface driven processes and, due to their complexity, modelling such coupled flow systems is a challenging task.

Different spatial scales can be employed to investigate coupled free-flow and porous-medium systems. At the microscale (pore scale) the pore structure is fully resolved (Fig. 1, left) and the flow in the entire fluid domain (free-flow region and pore space in the porous medium) is described by the Navier–Stokes equations with the no-slip condition at the boundary of the solid inclusions. However, computing the microscale flow field is infeasible for practical applications as it requires detailed information about porous-medium morphology and topography which is usually unknown. Even if information on the pore structure is available, e.g., obtained from tomographic analysis Blunt_17 ; Wildenschild_Sheppard_13 or known in advance for artificially produced composite materials (GeoDict, www.math2market.com), performing microscale simulations is computationally highly demanding.

From the macroscale perspective, the whole system is described as two different continuum flow domains (free flow, porous medium) separated by an equi-dimensional transition region or a lower-dimensional (co-dimension one) sharp interface. Different mathematical models are usually used in the two flow regions. Thus, either coupling conditions at the sharp fluid-porous interface or equations that apply at the transition zone between the two flow systems are needed Beavers_Joseph_67 ; Goyeau_Lhuillier_etal_03 ; Nield_09 ; OchoaTapia_Whitaker_95 ; Saffman . Correct specification of these conditions is essential for a complete and accurate mathematical description of flow and transport processes in compositional systems as well as for an accurate numerical simulation of applications Discacciati_Miglio_Quarteroni_02 ; Discacciati_Quarteroni_09 ; Jackson_Rybak_etal_12 ; Jaeger_Mikelic_09 ; Layton_Schieweck_Yotov_03 . The Navier–Stokes equations, which in case of low Reynolds numbers can be simplified to the Stokes equations, are typically applied to describe fluid flow in the free-flow domain, while Darcy’s law or its extension is used for the flow through the porous medium. Depending on the flow regime, other models can be applied, e.g., shallow water equations or kinematic wave approximation of the Saint-Venant equation for surface flows Dawson_08 ; Sochala_Ern_Piperno_09 ; Reuter_etal_19 and the Brinkman equation, the Forchheimer equation, the Richards equation or the full two-phase porous-medium flow models for subsurface flows Brinkman_47 ; Cimolin_Discacciati_13 ; Mosthaf_Baber_etal_11 ; Rybak_etal_15 . Different sets of coupling conditions are used in the literature at the fluid-porous interface Beavers_Joseph_67 ; OchoaTapia_Whitaker_95 ; Carraro_etal_15 ; Nield_09 ; Angot_etal_17 . Some authors distinguish between two qualitatively different flow directions: near parallel flow, where the velocity of the free fluid is much larger than the filtration velocity in the porous medium and near normal flow, where velocities in the two regions have comparable magnitudes. However, these conditions cannot be applied to a general flow situation.

The macroscopic level of flow system description is the most common one used for numerical simulations of coupled systems in practical applications Arbogast_Brunson_07 ; Cimolin_Discacciati_13 ; HanspalWaghodeNassehiWakeman ; Iliev_Laptev_04 ; Mosthaf_Baber_etal_11 ; Riviere . Microscale (pore-scale resolved) numerical simulations are typically used to validate the macroscale models.

In this paper, we consider the Stokes–Darcy problem and analyse different coupling conditions at the fluid-porous interface. Most often, the conservation of mass, the balance of normal forces and the Beavers–Joseph–Saffman condition for the tangential velocity are used both for numerical analysis and simulations, e.g., Beavers_Joseph_67 ; Saffman ; Goyeau_Lhuillier_etal_03 ; Nield_09 ; Discacciati_Quarteroni_09 ; GiraultRiviere ; Layton_Schieweck_Yotov_03 ; Jaeger_Mikelic_09 . The Beavers–Joseph–Saffman condition that describes a tangential velocity slip at the fluid-porous interface was obtained for flows parallel to the porous layer. Nevertheless, this condition is often applied to other flow directions, e.g., HanspalWaghodeNassehiWakeman . Although other interface conditions are available for the Stokes–Darcy problem, they are either valid only for specific flow situations Carraro_etal_15 or contain parameters which are difficult to determine ValdesParada_etal_09 . There were also several attempts to validate the coupled models and the interface conditions, e.g., Lacis_Bagheri_17 ; Zampogna_Bottaro_16 ; Bars_Worster_06 ; Yang_etal_19 . However, to our knowledge no systematic study has been published. In this paper, we consider different interface conditions for the coupled Stokes–Darcy problem. We validate and calibrate the macroscale coupled model using pore-scale resolved numerical simulations. These fully resolved numerical simulations are performed using the lattice Boltzmann method, a numerical method for simulating fluid dynamics that is especially suited for parallel computations with high spatial and temporal resolutions in complex geometries. It thus provides first-principle insights into the physical flow processes in both domains and the flow behaviour at the interface region.

The paper is organised as follows. In Section 2 we provide the geometrical setting and the flow system description. The macroscale coupled model with different sets of interface conditions is presented in Section 3. The calculation of effective parameters for the macroscale model by means of homogenisation is discussed in Section 4 and the lattice Boltzmann method is described in Section 5. Numerical simulation results for model validation and calibration are given in Section 6. Finally, the conclusion is presented in Section 7.

2 Geometrical setting

We consider a free-flow region $\Omega_{\mathrm{ff}}\subset\mathbb{R}^{2}$ and a porous medium $\Omega_{\text{pm}}\subset\mathbb{R}^{2}$ with periodic arrangement of the solid obstacles (Fig. 1, left). The fluid is considered to be present in a single phase and to fully saturate the porous medium. Additionally, it is assumed to be incompressible and to have a constant viscosity. The flow is considered as being laminar and supposing a single chemical species, no compositional effects need to be modelled. The temperature is assumed constant, such that it is not required to model any energy balance equation.

At the microscale, we solve the Stokes equations in the free-flow domain and in the pore space of the porous medium, where no-slip boundary conditions are applied on the solid-fluid boundary. At the macroscale, the system is described as two different continuum flow domains separated by the sharp interface $\Gamma$ (Fig. 1, right). We solve the Stokes and Darcy equations in the two flow domains and impose different interface conditions on the fluid-porous interface (Sect. 3). These interface conditions are validated for different porous-medium morphologies and different flow directions.

3 Macroscale coupled model formulation

Fluid flow in the free-flow region is described by the Stokes equations, the flow through the porous medium is described by Darcy’s law. Different sets of interface conditions are imposed on the sharp fluid-porous interface. Effective parameters such as permeability, porosity and boundary layer constants are needed to describe the porous-medium system and the interface conditions. These parameters are calculated in Section 6 for the geometrical configurations studied in the paper.

3.1 Free-flow model

The fluid is considered to be incompressible and therefore the mass conservation equation becomes

[TABLE]

The Stokes equation describes the fluid flow

[TABLE]

where $\rho$ is the fluid density, ${\boldsymbol{\mathrm{v}}}$ is the fluid velocity, $p$ is the pressure, ${\boldsymbol{\mathrm{g}}}$ is the gravitational acceleration, $\boldsymbol{\mathsf{T}}({\boldsymbol{\mathrm{v}}},p)=2\mu\boldsymbol{\mathsf{D}}\left({\boldsymbol{\mathrm{v}}}\right)-p{\boldsymbol{\mathsf{I}}}$ is the stress tensor, $\mu$ is the dynamic viscosity, $\displaystyle\boldsymbol{\mathsf{D}}\left({\boldsymbol{\mathrm{v}}}\right)=\tfrac{1}{2}\left(\nabla{\boldsymbol{\mathrm{v}}}+\left(\nabla{\boldsymbol{\mathrm{v}}}\right)^{{}^{\mathrm{T}}}\right)$ is the rate of strain tensor and ${\boldsymbol{\mathsf{I}}}$ is the identity tensor.

On the external boundary of the free-flow domain $\Gamma_{\mathrm{ff}}=\partial\Omega_{\mathrm{ff}}\setminus\Gamma$, the following boundary conditions are imposed

[TABLE]

where ${\boldsymbol{\mathrm{n}}}_{\mathrm{ff}}$ is the unit outward normal vector from domain $\Omega_{\mathrm{ff}}$ on its boundary and $\Gamma_{\mathrm{ff}}=\Gamma_{D,{\mathrm{ff}}}\cup\Gamma_{N,{\mathrm{ff}}}\cup\Gamma_{\text{out},{\mathrm{ff}}}$.

3.2 Porous-medium model

We consider Darcy’s flow equations in the porous-medium domain

[TABLE]

where Darcy’s velocity is given by

[TABLE]

Here ${\boldsymbol{\mathrm{v}}}$ is the fluid velocity through the porous medium, $p$ is the fluid pressure, $\boldsymbol{\mathsf{K}}$ is the intrinsic permeability of the porous medium, $\mu$ is the dynamic viscosity of the fluid, $\rho$ is the fluid density, $q$ is the source term and ${\boldsymbol{\mathrm{g}}}$ is the gravitational acceleration. The permeability tensor $\boldsymbol{\mathsf{K}}$ is symmetric, positive definite and bounded. Its values are computed numerically for the considered settings in Section 6. In this paper, we neglect gravitational effects setting ${\boldsymbol{\mathrm{g}}}={\boldsymbol{\mathrm{0}}}$.

We prescribe the following boundary conditions on the external boundary of the porous-medium domain $\Gamma_{\mathrm{pm}}=\partial\Omega_{\mathrm{pm}}\setminus\Gamma$:

[TABLE]

where ${\boldsymbol{\mathrm{n}}}_{\mathrm{pm}}$ is the unit outward normal vector from domain $\Omega_{\mathrm{pm}}$ on its boundary, $\Gamma_{\mathrm{pm}}=\Gamma_{D,{\mathrm{pm}}}\cup\Gamma_{N,{\mathrm{pm}}}$, $\Gamma_{D,{\mathrm{pm}}}\cap\Gamma_{N,{\mathrm{pm}}}=\emptyset$, and $\Gamma_{D,{\mathrm{pm}}}\neq\emptyset$.

3.3 Interface conditions

In addition to the boundary conditions prescribed on the external boundary of the coupled domain, a set of interface conditions has to be defined on $\Gamma$ in order to obtain a closed problem formulation. Different sets of interface conditions have been proposed in the literature depending on the flow direction.

3.3.1 Flow parallel to the interface

When the free flow is almost parallel to the porous medium, the conservation of mass, the balance of normal forces and a variation of the Beavers–Joseph condition is typically applied.

The conservation of mass across the interface $\Gamma$ is given by

[TABLE]

where ${\boldsymbol{\mathrm{n}}}$ is the unit normal vector at the fluid-porous interface $\Gamma$ pointing outward from the free-flow domain $\Omega_{\mathrm{ff}}$ (Fig. 1, right).

The balance of normal forces at the interface $\Gamma$ is

[TABLE]

where ${\boldsymbol{\mathrm{n}}}_{\mathrm{ff}}={\boldsymbol{\mathrm{n}}}=-{\boldsymbol{\mathrm{n}}}_{\mathrm{pm}}$.

There exist several possibilities for the interface condition on the tangential component of the velocity Nield_09 . The Beavers–Joseph interface condition Beavers_Joseph_67 for the difference between the tangential component of the free-flow velocity and the porous-medium velocity can be written as

[TABLE]

where $\alpha_{\text{BJ}}>0$ is the Beavers–Joseph parameter, ${\boldsymbol{\mathrm{\tau}}}_{j}$ is a unit vector tangential to the interface, $j=1,\ldots,d-1$, and $d$ is the number of space dimensions.

Saffman Saffman simplified the Beavers–Joseph condition (9) as follows

[TABLE]

arguing that the tangential velocity of a fluid through a porous medium is small and thus can be neglected.

Jones Jones_73 considered the symmetrised rate of strain tensor $\boldsymbol{\mathsf{D}}$ in the Beavers–Joseph–Saffman interface condition

[TABLE]

Besides the classical interface conditions (7), (8) and (9)–(11), which have been formulated based on heuristic arguments and were justified experimentally, in Jaeger_etal_01 the following interface conditions on $\Gamma$ were proposed

[TABLE]

where ${\boldsymbol{\mathrm{v}}}=(v_{1},v_{2})$ and $\varepsilon$ is the ratio of the length scales (Sect. 4). These coupling conditions are derived mathematically using homogenisation theory and boundary layer theory for flows parallel to the fluid-porous interface $\Gamma$. The first condition in (12) differs from the classical conservation of mass condition (7). However, it seems to be reasonable for flows parallel to the porous layer, because in this case the normal component of the Darcy velocity is of order $\varepsilon^{2}$ (see Eq. (18)) and thus can be neglected. The second condition in (12) links pressure fields in the porous domain and the free-flow domain, and is a variation of the classical balance of normal forces (8). The last equation in (12) is a modification of the Beavers–Joseph–Saffman law Beavers_Joseph_67 ; Saffman with $\displaystyle\alpha_{\mathrm{BJ}}^{-1}\boldsymbol{\mathsf{K}}^{1/2}=-\varepsilon C_{1}^{\text{bl}}{\boldsymbol{\mathsf{I}}}$. The effective coefficients $\displaystyle C_{1}^{\text{bl}}$ and $\displaystyle C_{\omega}^{\text{bl}}$ are determined through an auxiliary boundary layer problem. These coefficients and the parameter $\varepsilon$ are provided for different geometrical configurations in Section 6.

3.3.2 Flow perpendicular to the interface

For flows almost perpendicular to the porous medium, the following effective interface laws on $\Gamma$ are derived in Carraro_etal_15 using the theory of homogenisation and boundary layers

[TABLE]

where $M_{1}^{\text{bl}}$ is a boundary layer corrector constant. However, these conditions are valid only for a very specific setting (given geometrical configuration and boundary conditions) and cannot be applied to a general infiltration problem. There is a lack of physically consistent interface conditions for flows perpendicular to the porous layer and moreover for arbitrary flow directions.

4 Effective properties of the coupled system

To obtain effective properties for coupled macroscale models (effective permeability and boundary layer constants), we use the theory of homogenisation and boundary layers Carraro_etal_15 ; Hornung_97 ; Jaeger_Mikelic_09 . Therefore, we consider periodic porous media and assume the separation of length scales $\varepsilon=\ell/\mathcal{L}\ll 1$, where $\ell$ is the characteristic pore size and $\mathcal{L}$ is the length of the domain of interest (Fig. 2).

We consider a coupled domain $\Omega=\Omega_{\mathrm{ff}}\cup\Omega_{\mathrm{pm}}$, consisting of a free-flow region $\Omega_{\mathrm{ff}}\subset\mathbb{R}^{2}$ and a porous medium $\Omega_{\text{pm}}\subset\mathbb{R}^{2}$ with periodic arrangement of the solid obstacles (Fig. 1). The flow region $\Omega^{\varepsilon}=\Omega_{\text{ff}}\cup\Gamma\cup\Omega_{\text{pm}}^{\varepsilon}$ consists of the free-flow domain $\Omega_{\mathrm{ff}}$, the permeable interface $\Gamma$ and the pore space of the porous-medium domain $\Omega_{\text{pm}}^{\varepsilon}$. The porous-medium geometry is defined by a periodic repetition of the scaled unit cell $Y=(0,1)^{2}$ in $\Omega_{\mathrm{pm}}$ (Fig. 2).

4.1 Microscale model

The fluid flow in $\Omega^{\varepsilon}$ is governed by the non-dimensional steady Stokes equations

[TABLE]

completed with the no-slip condition on the fluid-solid interface

[TABLE]

and the appropriate boundary conditions on the external boundary $\partial\Omega$.

For the lid driven cavity over the porous bed described in Section 6.1, we consider

[TABLE]

and for the infiltration problem validated in Section 6.2, we have

[TABLE]

where $p_{b}=100$.

4.2 Homogenisation

We obtain the corresponding macroscale formulation of the coupled system by studying the behaviour of the solutions to the microscopic problem (14)–(16) or (14), (15) and (17) when $\varepsilon\rightarrow 0$. In the limit, the equations in the free-flow region $\Omega_{\mathrm{ff}}$ remain unchanged and we need to homogenise the porous region. To derive the limit problem in $\Omega^{\varepsilon}_{\text{pm}}$, we use the idea of homogenisation with two-scale asymptotic expansions (Hornung_97, , chap. 1 and 3):

[TABLE]

where ${\boldsymbol{\mathrm{y}}}={\boldsymbol{\mathrm{x}}}/\varepsilon$ is the fast variable. We follow the classical procedure of homogenisation to derive Darcy’s law (5) as the upscaled problem in $\Omega_{\text{pm}}$ completed by the appropriate boundary conditions.

The permeability tensor $\boldsymbol{\mathsf{\widetilde{K}}}$ is given by

[TABLE]

where $\textbf{w}^{j}=(w_{1}^{j},w_{2}^{j})$ and $\pi^{j}$ are the solutions to the following cell problems for $j=1,2$:

[TABLE]

The unit cell $Y$ is presented in Figure 2. To obtain the effective permeability $\boldsymbol{\mathsf{K}}$, the tensor $\boldsymbol{\mathsf{\widetilde{K}}}$ has to be scaled with $\varepsilon^{2}$.

5 Lattice Boltzmann method

The lattice Boltzmann method (LBM) is a modern approach for simulating fluid dynamics. While the LBM historically evolved from lattice gas automata Wolf-Gladrow_00 , it can be regarded as a special discretisation of the continuous Boltzmann equation He_97 that describes the dynamics of a gas on a mesoscopic scale. On a macroscopic scale, however, the Boltzmann equation leads to the equations of fluid dynamics, i.e., mass, momentum and energy conservation. Therefore, by solving the Boltzmann equation for certain scenarios, one can also obtain the corresponding solutions to the Navier–Stokes equations Krueger_LBM .

5.1 Advantages

In comparison to conventional methods for simulating fluid dynamics, such as the finite difference, finite volume or finite element method, the LBM’s main advantages lay in its ease of implementation and its suitability for complex flow scenarios. As the most compute-intensive operations are local, i.e., restricted to neighbouring lattice cells, the LBM is inherently applicable for massively parallel computing. Furthermore, the LBM is well-suited for simulating particulate flows Rettinger_17 ; Kuron_16 , multi-phase flows Huang_15 and flows through complicated geometries such as porous media Fattahi_etal_16 ; Fattahi_etal_16b .

5.2 Description

Our introduction to the LBM follows Krueger_LBM ; the interested reader is referred to this reference and the references therein. The LBM’s most fundamental quantities are the so-called particle distribution functions $f_{i}({\boldsymbol{\mathrm{x}}},t)$ that represent the density of virtual particles with velocity ${\boldsymbol{\mathrm{c}}}_{i}=(c_{ix},c_{iy})$ at position ${\boldsymbol{\mathrm{x}}}$ and time $t$. The index $i\in\left[0;q\right)$ in $f_{i}$ refers to the corresponding index of velocity ${\boldsymbol{\mathrm{c}}}_{i}$ in a discrete set of velocities $\{{\boldsymbol{\mathrm{c}}}_{i},w_{i}\}$ with weighting coefficients $w_{i}$.

Velocity sets are denoted as D$d$Q$q$ with $d$ being the number of spatial dimensions and $q$ referring to the velocity set’s number of velocities. Velocity sets are a trade-off between accuracy and memory consumption. Different velocity sets and their corresponding coefficients can be found in the LBM literature such as Krueger_LBM . Here we limit ourselves to two spatial dimensions such that we use the well-established D2Q9 velocity set. In each velocity set, the constant $c_{s}$ determines the relation between the fluid pressure $p$ and the fluid density $\rho$ as $p=c_{s}^{2}\rho$. Due to this relation, $c_{s}$ is called the lattice speed of sound. In most velocity sets used to solve the Navier–Stokes equations, it can be expressed as $c_{s}=\sqrt{1/3}\Delta x/\Delta t$.

The distribution function $f_{i}$ is only defined at discrete points ${\boldsymbol{\mathrm{x}}}$ in a square lattice with lattice spacing $\Delta x$. Since the LBM’s origins lay in lattice gas automata, these discrete points are also often referred to as lattice cells. Similarly, $f_{i}$ is only defined at discrete times $t$ with time step $\Delta t$. As it is common in the LBM literature, we use lattice units which is an artificial set of units that is scaled such that the relations $\Delta x=1$ and $\Delta t=1$ hold.

In kinetic theory, one obtains the macroscopic mass and momentum density by taking the moments of the Boltzmann equation’s distribution functions. Similarly, one can find the mass density $\rho$ and the momentum density $\rho{\boldsymbol{\mathrm{v}}}$ with fluid velocity ${\boldsymbol{\mathrm{v}}}$ by

[TABLE]

The lattice Boltzmann equation

[TABLE]

describes the advection of particles $f_{i}({\boldsymbol{\mathrm{x}}},t)$ that travel with velocity ${\boldsymbol{\mathrm{c}}}_{i}$ to a neighbouring lattice cell ${\boldsymbol{\mathrm{x}}}+{\boldsymbol{\mathrm{c}}}_{i}\Delta t$ in the next time step $t+\Delta t$ (Fig. 3, left). After having moved from one lattice cell to another, the particles are affected by a collision operator $\Theta_{i}({\boldsymbol{\mathrm{x}}},t)$ that models particle collisions by redistributing particles among $f_{i}$ (Fig. 3, right).

The simplest collision operator that can be used for solving the Navier–Stokes equations is the Bhatnagar–Gross–Krook (BGK) operator. It linearly relaxes particle distributions towards an equilibrium $f_{i}^{\text{eq}}$ by

[TABLE]

with a relaxation rate determined by the relaxation time $\tau$. In contrast to other collision operators, the BGK operator is based on only one relaxation time and it is therefore also called single-relaxation-time (SRT) collision operator. Due to its simplicity, the BGK collision operator is often the standard choice for LBM simulations. However, in Pan_06 , it was found to be subject to a non-physical dependence between the fluid viscosity and the boundary locations, which is especially problematic in scenarios with large boundary areas such as flow through porous media. Therefore, for the LBM simulations in this paper, we use the more accurate two-relaxation-time (TRT) collision operator Ginzburg_08a ; Ginzburg_08b .

The particle distribution functions’ equilibria

[TABLE]

are derived by a Hermite polynomial expansion of the continuous Maxwell-Boltzmann distribution Krueger_LBM .

The moments of $f_{i}^{\text{eq}}$ are identical to those of the corresponding $f_{i}$, i.e., $\Sigma_{i}f_{i}^{\text{eq}}=\Sigma_{i}f_{i}=\rho$ and $\Sigma_{i}{\boldsymbol{\mathrm{c}}}_{i}f_{i}^{\text{eq}}=\Sigma_{i}{\boldsymbol{\mathrm{c}}}_{i}f_{i}=\rho{\boldsymbol{\mathrm{v}}}$. Since $f_{i}^{\text{eq}}$ only depends on the local quantities density $\rho$ and fluid velocity ${\boldsymbol{\mathrm{v}}}$ (Eq. (21)), global information exchange is not needed and the LBM is therefore well-suited for parallel computing.

Using the Chapman–Enskog analysis Krueger_LBM , one can show that the lattice Boltzmann equation (22) has a macroscopic behaviour that converges asymptotically to the solutions of the Navier–Stokes equations. The fluid’s kinematic shear viscosity $\nu$ is then a function of the relaxation time $\tau$ and obtained by

[TABLE]

6 Model validation and calibration

In this section, we validate and calibrate coupled Stokes–Darcy models using pore-scale resolved simulations for different geometrical configurations and different flow regimes.

We consider the free-flow region $\Omega_{\mathrm{ff}}=[0,1]\times[0,1]$ and the porous medium $\Omega_{\mathrm{pm}}=[0,1]\times[-0.5,0]$ and study two test cases (lid driven cavity over a porous bed and infiltration problem). The porous medium is set to be periodic for both test cases such that we can compute the effective parameters easily using homogenisation theory. Moreover, we analyse two different periodic geometries.

The porosity of the medium is taken to be $\phi=0.4$. Therefore, the radius of solid inclusions is given by $r=\varepsilon\sqrt{(1-\phi)/\pi}$. We solve the unit cell problems (20) and compute the effective permeability tensor for different geometrical settings according to equation (19). To compute the boundary layer constants, we solve the Stokes problem in a perforated vertical stripe containing unit cells (see e.g. Carraro_etal_15 ). For numerical simulations we use the software package FreeFEM++ with P2/P1 finite elements MR3043640 . The dynamic viscosity of the fluid is scaled for both test cases to $\mu=1$.

The macroscale problem is discretised using the finite volume method on staggered grids VersteegMalalasekra . The computational domains $\Omega_{\mathrm{ff}}$ and $\Omega_{\mathrm{pm}}$ are partitioned into equal blocks with grid size $h=1/800$, and the grids are conforming at the interface.

We apply large scale simulations using the LBM described in Section 5 to compute the microscale (pore-scale) velocity field and compare the macroscale simulations with different sets of interface conditions against these pore-scale resolved models.

The LBM simulations are implemented and performed with the open-source software-framework waLBerla Godenschwager_2013 ; Schornbaum_18 (www.walberla.net). We use the $D2Q9$ velocity set and the TRT collision operator. We set the TRT collision operator’s even order relaxation time $\tau=1$ and choose its magic parameter $\Lambda=3/16$ Ginzburg_08a ; Ginzburg_08b . With equation (25), the fluid’s kinematic viscosity in lattice units becomes $\nu=1/6$. We consider the simulation to have reached steady state when the relative temporal change in the absolute value of the fluid’s maximum velocity in $x_{2}$-direction is below $10^{-4}$.

In the following sections, we analyse the sensitivity of the coupled Stokes–Darcy models to the location of the fluid-porous interface $\Gamma$, the choice of the effective parameters and the interface conditions.

6.1 Lid driven cavity over the porous bed

For the lid driven cavity problem, we consider two geometrical configurations, a channelised one schematically presented in Figure 4 and a staggered arrangement of the solid inclusions presented in Figure 13. In both cases there are 40 solid inclusions in $x_{1}$-direction and 20 inclusions in $x_{2}$-direction in the porous domain $\Omega_{\mathrm{pm}}$. We define the Reynolds number

[TABLE]

using the $x_{1}$-component of the velocity at the lid $v_{\text{lid}}$, the length of the domain $\mathcal{L}$ in $x_{1}$-direction (Fig. 2) and the kinematic fluid viscosity $\nu$.

6.1.1 Channelised geometrical configuration

In this section, we consider the channelised geometrical configuration (Fig. 4). For this geometry, we solve the cell problems (20), compute the permeability tensor according to (19) and obtain

[TABLE]

The effective permeability tensor of the porous medium is given by $\boldsymbol{\mathsf{K}}=\varepsilon^{2}\boldsymbol{\mathsf{\widetilde{K}}}$, where the actual scale ratio is $\varepsilon=1/40$.

The corresponding macroscale model (1), (2), (4), (5) is closed by the following boundary conditions on the external boundary

[TABLE]

In the LBM simulation, the whole domain is resolved by $2400\times 3600$ lattice cells such that a solid inclusion has a diameter of about 52 lattice cells. We apply simple bounce back boundary conditions Krueger_LBM at the solid inclusions and at the external boundaries of the domain with respect to (28). In order to approximate Stokes flow, we set $Re=0.01$.

Two sets of interface conditions are analysed for the macroscale problem: (i) the classical conditions (7), (8), (11) and (ii) the interface conditions derived using homogenisation (12).

(i) Classical interface conditions

In this section, we study the sensitivity of the coupled macroscale Stokes–Darcy (SD) model (1), (2), (4), (5) with the classical interface conditions (7), (8), (11) to the location of the sharp fluid-porous interface and to the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$.

We evaluate the cross-sections at $x_{1}=0.4875$ (intersecting solid inclusions) and $x_{1}=0.5$ (between solid inclusions), where the flow is almost parallel to the interface and near the boundary at $x_{1}=0.025$, where the flow is almost perpendicular. The locations of the sharp fluid-porous interface for the macroscale model are $x_{2}=-0.01$, $x_{2}=0$ (no shift) and $x_{2}=0.01$.

Velocity profiles ($x_{1}$- and $x_{2}$-components) in and close to the horizontal centre of the domain are presented in Figures 6, 6. The Beavers–Joseph parameter is taken $\alpha_{\mathrm{BJ}}=0.5$ and different interface locations are considered: $(0,1)\times\{-0.01\}$ (shift below, $s=b$), $(0,1)\times\{0\}$ (no shift, $s=0$) and $(0,1)\times\{0.01\}$ (shift above, $s=t$). The profiles differ in the near-interface region and the best fit of the macroscale model to the pore-scale resolved simulations is obtained for $s=0$. This means that the sharp fluid-porous interface $\Gamma$ should be placed directly on top of the first row of solid inclusions (Fig. 4).

Figure 7 provides velocity profiles for different values of the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$. The best fit is obtained for $\alpha_{\mathrm{BJ}}=0.5$. This will be also confirmed later using interface conditions (12) obtained by homogenisation (Fig. 11–12). In the literature, $\alpha_{\mathrm{BJ}}=1$ is usually used, which is not the optimal choice.

Velocity profiles near the left boundary are provided in Figures 9, 9. Here, the location of the sharp interface and also the value of the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$ has almost no effect on the $x_{1}$-component of the velocity (Fig. 9). This is due to the fact that the velocity is almost perpendicular to the porous layer and the magnitude of the horizontal velocity is very small. However, for the vertical component of the velocity, the sharp interface location is important (Fig. 9).

(ii) Interface conditions based on

homogenisation

In this section, we study the sensitivity of the coupled macroscale Stokes–Darcy model (1), (2), (4), (5) with the interface conditions (12) proposed in Jaeger_etal_01 ; Jaeger_Mikelic_09 .

To apply these interface conditions we need to compute the boundary layer constants for our geometry (Fig. 4). For these constants, we obtain the following values

[TABLE]

In Figures 11, 11 we provide velocity profiles for the lid driven cavity problem in and close to the horizontal centre of the domain. For the macroscale model we observed that the interface conditions derived by homogenisation (profile: SD (homogen.)) fit the best to the classical coupling conditions for $\alpha_{\mathrm{BJ}}=0.5$. For such flow problems, we advise to use condition (12) instead of trying to find the optimal value for $\alpha_{\mathrm{BJ}}$.

Since the Beavers–Joseph–Saffman condition is the interface condition for the horizontal component of the velocity, the vertical component is only poorly sensitive to a change in the value of the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$. We observed that the interface condition for the $x_{2}$-component of the velocity, given by the first equation in (12), is nevertheless reasonable for such flow problems.

To demonstrate the sensitivity of the coupled Stokes–Darcy model to the choice of the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$ in (11), we make cross-sections of the horizontal component of the velocity at the fluid-porous interface ($x_{2}=0$). The numerical simulation results for different values of $\alpha_{\mathrm{BJ}}$ and also for the interface conditions derived by homogenisation (12) are presented in Figure 12.

For the macroscale simulations, the free-flow velocity at the interface is plotted. The averaged pore-scale simulations (Fig. 12, profile: LBM (avg.)) suggest that the choice $\alpha_{\mathrm{BJ}}=1$ which is commonly used in the literature is not optimal one. Again, as shown in Figures 7,  11,  11, we observed that $\alpha_{\mathrm{BJ}}=0.5$ as well as the interface conditions (12) provide a much better fit.

6.1.2 Staggered geometrical configuration

Since the Beavers–Joseph parameter $\alpha_{\mathrm{BJ}}$ contains the information about the geometry of the interfacial region, we consider another geometrical configuration for the porous bed. In this case, the porous medium consists of solid inclusions that are arranged in a periodic but staggered manner (Fig. 13). The porosity is kept the same as for the channelised setting $\phi=0.4$. The number of solid inclusions is also the same: 40 inclusions in $x_{1}$-direction and 20 inclusions in $x_{2}$-direction.

For this porous-medium geometry, we obtain the effective permeability tensor $\boldsymbol{\mathsf{K}}=\varepsilon^{2}\boldsymbol{\mathsf{\widetilde{K}}}$ in the same way as before with

[TABLE]

where $\varepsilon=1/20$ and the unit cell $Y$ contains several solid inclusions (Fig. 13, right).

We analysed the macroscale problem (1), (2), (4), (5) with the classical interface conditions (7), (8), (11) and the appropriate boundary conditions at the external boundary.

Figure 14 provides velocity profiles of the horizontal component of the velocity in the middle of the domain. Again, as we observed for the channelised geometry, the macroscale simulation results with $\alpha_{\mathrm{BJ}}=0.5$ fit much better to the pore-scale resolved simulations than the results with $\alpha_{\mathrm{BJ}}=1$ that is commonly used in the literature.

6.2 Infiltration

For the infiltration problem we consider only the staggered geometrical configuration of the porous medium (Fig. 13) to avoid a channel flow in the vertical direction between columns of solid inclusions (Fig. 4). Therefore, the permeability tensor is identical with the one given in (30).

The boundary conditions for the macroscale model given by equations (1), (2), (4), (5), (7), (8), (11) for the infiltration problem read

[TABLE]

Here, we define the Reynolds number

[TABLE]

where $v_{\text{in,max}}$ is the absolute maximum value of the inflow velocity ($x_{2}$-component). Similar as for the lid driven cavity problem (Sect. 6.1), in the LBM simulation, the whole domain is resolved by $2400\times 3600$ lattice cells. The solid inclusions have a diameter of around 52 lattice cells. To resemble (31), we use simple bounce back boundary conditions at the domain-boundaries in $x_{1}$-direction, at the inflow and at the solid inclusions. At the bottom of the domain in $x_{2}$-direction, we impose pressure boundary conditions based on the anti-bounce-back approach Krueger_LBM and set the lattice density to $\rho=1$. Again, $Re=0.01$ is chosen in order to approximate Stokes flow.

We make cross-sections for the velocity and pressure in the horizontal centre of the domain at $x_{1}=0.5$ and consider different locations of the sharp fluid-porous interface $\Gamma$ for the macroscale model. As in Section 6.1, the sharp interface is located at $x_{2}=-0.01$ ($s=b$), $x_{2}=0$ (no shift, $s=0$) and $x_{2}=0.01$ ($s=t$).

In Figure 16 we present the vertical component of the velocity computed for the macroscale model (SD) and the LBM simulations. We normalised the velocity by its maximal absolute value $v_{\text{in,max}}=0.1$. Although we used for the macroscale simulations the interface conditions (7), (8), (11), which are valid for parallel flows, we get a quite good fit to the results of the pore-scale resolved models. This is due to the fact that the $x_{1}$-component of the velocity is almost zero. Therefore, we do not provide the profiles for the horizontal velocity component.

For the infiltration problem, we observed that the pressure field is very sensitive to the location of the sharp interface (Fig. 16). Again, we obtain the best fit of the macroscale to the pore-scale model (Fig. 16, $s=0$) for the interface location presented in Figure 13. However, the differences between the pressure field obtained from the LBM simulations and the one from the best fitting interface location of the macroscopic model are more significant compared to the differences we observed for the velocity profiles. Here, it is important to recall that the permeability tensor for the macroscale model is computed using the finite element framework FreeFEM++. We suspect that one cause of the more notable differences in the pressure profiles could be slight differences in the representation of the solid inclusions in the computational grids of LBM and FreeFEM++.

7 Conclusion

In this paper, we validated two sets of interface conditions for the coupled Stokes–Darcy problem: (i) the classical conservation of mass across the interface, the balance of normal forces and the Beavers–Joseph–Saffman interface condition, and (ii) the interface conditions derived by means of the homogenisation and boundary layer theory Jaeger_Mikelic_09 ; Jaeger_etal_01 . We considered different geometrical configurations (channelised and staggered arrangement of solid grains) and different flow problems (lid driven cavity over a porous bed and infiltration problem). All effective model parameters are computed numerically for the geometrical settings examined in the paper. To be able to compute these parameters using homogenisation theory, we considered periodic porous media. The interface conditions are validated and calibrated numerically by comparing macroscale simulation results against pore-scale simulations. As a pore-scale resolved model we used large scale lattice Boltzmann simulations.

Numerical simulation results demonstrate sensitivity of the coupled Stokes–Darcy problem to the choice of the effective model parameters, the location of the sharp interface and the interface conditions. We have tested a number of characteristic flow regimes in low Reynolds number cases and in all of them $\alpha_{\mathrm{BJ}}=0.5$ is proved to be the best candidate for the Beavers–Joseph parameter. We observed that the commonly used parameter $\alpha_{\mathrm{BJ}}=1$ in the Beavers–Joseph–Saffman interface condition (11) is often not the optimal choice. Moreover, the sharp fluid-porous interface should be located directly on the top of the first row of solid inclusions.

Acknowledgements.

The authors thank Ivan Yotov, Jim Magiera and Christoph Rettinger for the valuable discussions.

I. Rybak and E. Eggenweiler thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) for supporting this work by funding SFB 1313, Project Number 327154368.

C. Schwarzmeier and U. Rüde thank the DFG for supporting project RU 422/27 and the Bundesministerium für Bildung und Forschung (BMBF, Federal Ministry of Education and Research) for supporting project SKAMPY (01IH15003A). The authors are grateful to the Regionales Rechenzentrum Erlangen (www.rrze.de) for providing access to supercomputing facilities.