2.1 Direct numerical simulations

The two-dimensional Stokes equations are solved with a no-slip condition imposed at the cylinder surfaces as well as on the vertical and bottom walls of the cavity. The equations are given by,

(2.1) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\boldsymbol{u} & = & \displaystyle 0,\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u} & = & \displaystyle 0,\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle \boldsymbol{u} & = & \displaystyle (0,0)\quad \text{on}~\unicode[STIX]{x1D6E4}_{n},~\unicode[STIX]{x1D6E4}_{c},\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle \boldsymbol{u} & = & \displaystyle (U_{w},0)\quad \text{at}~y=H,\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ is fluid viscosity, $\unicode[STIX]{x1D6E4}_{n}$ denotes the bottom and side boundaries of the cavity and $\unicode[STIX]{x1D6E4}_{c}$ denotes the boundary of the cylinders in the porous bed. The computations are performed with FreeFEM++ (Hecht Reference Hecht2012), using a triangular mesh and a Taylor–Hood finite element space (P2+P1) for velocity and pressure. We set mesh spacing to $\unicode[STIX]{x0394}s_{1}=0.125l$ at the outer boundaries (cavity walls) of the domain, and $\unicode[STIX]{x0394}s_{2}=0.050l$ at the surface of the cylinders. We have carried out a simulation of the same configuration with half the mesh spacing ( $\unicode[STIX]{x0394}s_{1}=0.063l$ and $\unicode[STIX]{x0394}s_{2}=0.025l$ ), and observed that the slip velocity changes by 0.6 %.

Figure 2.

Left frame (a) shows the streamwise velocity component, whereas the right frame (b) shows wall-normal velocity component at $x/H=-0.1$ . Solid black lines depict direct numerical simulations of the porous cavity problem. Blue lines with circular markers correspond to the continuum model of the porous bed coupled to Stokes solver above the bed via interface conditions at $y_{s}/H=0.01$ (red dashed line), which is equivalent to $y_{s}/l=0.1$ in pore scale. Insets show velocity profiles near the interface.

In figures 2 and 3, we present the velocity profiles obtained from DNS (direct numerical simulations) with solid black lines. Figure 2 shows the streamwise and wall-normal velocity profiles for the fixed streamwise position $x/H=-0.1$ . The insets show the detailed microscale fluctuations of the velocities near the interface and the rapid transition to the macroscale velocity in the free-fluid region. This transition occurs in a thin layer near the top row of cylinders.

Figure 3.

In the top frame (a), we compare slip velocity $u_{s}$ prediction from the continuum model with DNS. The bottom frame (b) compares the penetration velocity $v_{p}$ between the two approaches. Velocities are sampled at the interface $y=y_{s}$ .

Figure 3 shows the velocity along $x$ at a virtual free-fluid–porous interface placed at $y_{s}/l=0.1$ . If $y_{s}$ would have been an interface with a rigid wall, these graphs would show the no-slip condition, i.e. zero velocity for both wall-normal and streamwise velocity components. However, at the interface with a porous medium, there is a slip velocity $u_{s}$ and a penetration velocity $v_{p}$ . The underlying structure of porous medium manifests as microscale oscillations in slip and penetration velocities. Apart from the microscale oscillations, we can observe that both velocity components exhibit macroscale variations. The negative slip velocity, which is induced by the spanwise vortex above the interface, has a maximum value $u_{s}/U_{w}=-0.0121$ at the centre of the cavity. Although the slip velocity is small compared to the bulk flow, it may have a significant physical effect on the characteristics of the overlying fluid. For example Rosti et al. (Reference Rosti, Cortelezzi and Quadrio2015) recently showed that slip velocities below 3 % had a significant effect on the flow statistics in a turbulent channel. Additionally, Carotenuto & Minale (Reference Carotenuto and Minale2013) showed that precise predictions of slip velocity can be essential to get accurate viscosity measurements from rheology tests. The penetration velocity shows a sinusoidal behaviour; for $x<0$ , there is a net mass transport from the pore region to the free-fluid region, whereas for $x>0$ the net mass flow is in the opposite direction. Similarly, the net momentum transport into the porous region is in opposite directions whether $x>0$ or $x<0$ .

The values of $u_{s}$ and $v_{p}$ , which are essential to capture the momentum and mass transport across the interface, depend both on the flow in the pores and on the microscale geometry of the pores. In the next section, we introduce a fully non-empirical method to compute $u_{s}$ and $v_{p}$ with an error of $O(l/H)$ without resorting to DNS of the full domain.

2.2 Simulation of homogenized equations

We start by replacing the full DNS domain with two rectangular domains, where the free-fluid region $\unicode[STIX]{x1D6FA}_{f}$ and porous region $\unicode[STIX]{x1D6FA}_{p}$ are separated by the interface $\unicode[STIX]{x1D6E4}$ , as shown in figure 1(c). In the free-fluid part, we do not employ homogenization and therefore the flow is governed by the Stokes equations

(2.5) $$\begin{eqnarray}\displaystyle -\unicode[STIX]{x1D735}\hat{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}} & = & \displaystyle 0,\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}} & = & \displaystyle 0,\end{eqnarray}$$

where $\hat{\boldsymbol{u}}$ and $\hat{p}$ are flow and pressure fields in $\unicode[STIX]{x1D6FA}_{f}$ , respectively. Dirichlet conditions are imposed for $\hat{\boldsymbol{u}}$ on the vertical side walls and the top wall of the cavity (as for the full DNS). The boundary condition at the interface $\unicode[STIX]{x1D6E4}$ in contact with the porous region is

(2.7) $$\begin{eqnarray}\hat{\boldsymbol{u}}=(\hat{u} _{s},\hat{v}_{p}),\end{eqnarray}$$

where the slip velocity and the penetration velocity depend on the flow in the porous region. In the porous part $\unicode[STIX]{x1D6FA}_{p}$ , the flow is governed by the well-known Darcy’s law

(2.8) $$\begin{eqnarray}\hat{\boldsymbol{u}}=-\frac{\overline{\overline{\unicode[STIX]{x1D646}}}^{\text{itr}}}{\unicode[STIX]{x1D707}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\hat{p},\end{eqnarray}$$

and mass conservation

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}=0,\end{eqnarray}$$

where $\overline{\overline{\unicode[STIX]{x1D646}}}^{\text{itr}}$ is the interior permeability tensor. It is convenient to combine mass conservation with Darcy’s law to arrive at a single equation for the pore pressure, which – assuming that permeability tensor $\overline{\overline{\unicode[STIX]{x1D646}}}^{\text{itr}}$ is constant over space and isotropic – reads

(2.10) $$\begin{eqnarray}\unicode[STIX]{x0394}\hat{p}=0.\end{eqnarray}$$

We complement this equation with homogenous Neumann conditions on the side walls and the bottom wall, which corresponds to zero transpiration. At the interface $\unicode[STIX]{x1D6E4}$ , we impose a Dirichlet condition with the pressure obtained from (2.5)–(2.6). This continuous pressure condition is valid up to $O(l/H)$ under the theoretical assumptions employed in this paper, which will be discussed in following sections. Since the particular porous bed we are investigating is isotropic (see figure 1 a,b), the continuity of pressure at the interface is also implied by works of Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) and Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013).

Returning to the velocity boundary conditions (2.7), that are required for solving Stokes system (2.5)–(2.6), we simply state the conditions of order $O(l/H)$ that will be derived in the next sections (with the final result in (5.11)). The penetration velocity component $\hat{v}_{p}$ is given by

(2.11) $$\begin{eqnarray}\hat{v}_{p}=-\frac{K_{cyl}}{\unicode[STIX]{x1D707}}\unicode[STIX]{x2202}_{y}\hat{p}^{-},\end{eqnarray}$$

where $K_{cyl}$ is the isotropic permeability of the porous medium consisting of a regular array of circular cylinders. Note that, although pressure is continuous at $\unicode[STIX]{x1D6E4}$ in our case, the pressure gradient is not necessarily continuous; $\unicode[STIX]{x2202}_{y}\hat{p}^{-}$ in (2.11) denotes the pressure gradient when approaching the interface from the porous bed. The condition for $v_{p}$ can also be obtained from mass conservation for a thin rectangular control volume around $\unicode[STIX]{x1D6E4}$ with periodic streamwise velocity on the vertical sides. The condition for slip velocity $\hat{u} _{s}$ is

(2.12) $$\begin{eqnarray}\hat{u} _{s}=-\frac{K_{s}}{\unicode[STIX]{x1D707}}\unicode[STIX]{x2202}_{x}\hat{p}^{-}+L_{s}(\unicode[STIX]{x2202}_{y}\hat{u} +\unicode[STIX]{x2202}_{x}\hat{v}).\end{eqnarray}$$

This expression is similar to the condition obtained empirically by Beavers and Joseph, except that $K_{s}$ is the interface permeability (e.g. $K_{s}\neq K_{cyl}$ ), related to a semi-permeable transition layer between the porous medium and the free fluid. Another difference with the BJ condition is that the strain term $\unicode[STIX]{x2202}_{x}\hat{v}$ is included in addition to $\unicode[STIX]{x2202}_{y}\hat{u}$ . It has been argued that the term $\unicode[STIX]{x2202}_{x}\hat{v}$ should be present for curved boundaries (Jones Reference Jones1973), but to the authors’ knowledge, it has not been derived earlier for flat interfaces, although it has be conjectured to exist by Nield (Reference Nield2009). The constant $L_{s}$ is related to the slip length in the Navier boundary condition. The constants appearing in boundary conditions (2.11)–(2.12) are provided by microscale simulations in interface cells, described in following sections. In order to provide an overview of the applicability of the derived boundary conditions, we summarize here briefly the practical limits that will be determined both theoretically and numerically in the remaining part of this paper. These limits are; (i) moderate scale separation $l/H\leqslant 0.1$ ; (ii) restriction on the Reynolds number based on the seepage velocity ( $Re_{d}\leqslant 1$ ); and (iii) restriction on the Reynolds number based on the lid velocity $U_{w}$ ( $Re_{f}\leqslant \unicode[STIX]{x1D716}^{-1}$ ). The corresponding Reynolds numbers are defined later.

We solve the set of equations (2.5)–(2.12) using FreeFEM++ with mesh spacing $\unicode[STIX]{x0394}s=0.125l$ . Figure 2 (blue curve with circular symbols) compares the obtained velocity profiles over a vertical slice to the DNS results, where one can observe an excellent agreement between the two. We note that the effective macroscale behaviour is captured, while underlying oscillations arising from the small-scale characteristics of porous bed are not modelled. The consequence of using an $O(l/H)$ accurate model in the interior (Darcy’s law) is that the diffusion process from the free flow to the pore flow (which defines the transition layer of height ${\sim}l$ ) is not captured. However, from the perspective of the free fluid, the macroscopic effect of the porous bed is essentially the same using fully resolved DNS and the continuum model. In figure 3(a) (blue curve with circular symbols), we also observe a good agreement between DNS and the continuum model for slip velocity over a horizontal slice, despite that the latter approach does not resolve the microscale dynamics between and around the cylinders. The model is able to predict the maximum slip velocity at the centre of the cavity $\hat{u} _{s}/U_{w}=-0.0119$ . Figure 3(b) compares the predicted penetration velocity and DNS results. There one can observe that, although microscale oscillations dominate the DNS results, the macroscale sinusoidal behaviour is correctly captured by the model.

To close this section, we want to point out that the effective problem is computationally much cheaper than DNS. The number of degrees of freedom used for the DNS in § 2.1 for the region below the interface is approximately $2.0\times 10^{5}$ , whereas for the two-domain approach in § 2.2 it is $1.4\times 10^{4}$ . This difference arises from coarser mesh in the porous region, as well as the reduced number of variables (the model equation defines pressure only). In more complex and three-dimensional cases the difference can be significantly larger, and only averaged models might be computationally feasible to solve numerically.

3.1 Dimensionless Navier–Stokes equations

The free-fluid region and the porous region are characterized by different spatial and temporal scales. We choose to non-dimensionalize the Navier–Stokes equations using the characteristic scales of the porous medium. In appendix A, the reader can find detailed analysis of these scales. The characteristic velocity of the flow in the porous region $U^{d}$ is

(3.1) $$\begin{eqnarray}U^{d}\sim \frac{l^{2}\unicode[STIX]{x0394}P}{\unicode[STIX]{x1D707}H},\end{eqnarray}$$

where $\unicode[STIX]{x0394}P$ is the characteristic macroscopic global pressure, $\unicode[STIX]{x1D707}$ is the fluid viscosity and $H$ and $l$ are the macroscopic and microscopic length scales, respectively. Consequently, we use the following relationships between dimensional (denoted with ‘tilde’) and dimensionless variables

(3.2) $$\begin{eqnarray}\tilde{u} _{i}=U^{d}u_{i},\quad \tilde{p}=\unicode[STIX]{x0394}Pp,\quad \tilde{x}_{i}=lx_{i},\quad \text{and}\quad \tilde{t}=\frac{l}{U^{d}}t.\end{eqnarray}$$

Here, time is non-dimensionalized with the convection time scale at the microscale. In order to simplify the notation, we use $x_{1}$ and $x$ , and $x_{2}$ and $y$ interchangeably. We may now write the Navier–Stokes equations in the following dimensionless form

(3.3) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}^{2}Re_{d}(\unicode[STIX]{x2202}_{t}u_{i}+u_{j}u_{i,j}) & = & \displaystyle -p_{,i}+\unicode[STIX]{x1D716}u_{i,jj},\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle u_{i,i} & = & \displaystyle 0,\end{eqnarray}$$

where $Re_{d}=\unicode[STIX]{x1D70C}_{f}U^{d}H/\unicode[STIX]{x1D707}$ is the Darcy Reynolds number. The different order of the scale separation parameter $\unicode[STIX]{x1D716}$ in front of the terms provides an estimate of the relative magnitude of the terms within the porous region; pressure force plays a dominant role and the inertial force is a higher-order effect. This conclusion holds if $Re_{d}\leqslant 1$ , which is an assumption for the presented approach. Technically, equations (3.3)–(3.4) hold also in the free-fluid region, although the relative magnitude between the terms is not anymore characterized by the scale separation parameter; the scaling (3.2), except for pressure, is not suited for the free flow.

3.2 Decomposition of the flow field

We continue by choosing an interface $\unicode[STIX]{x1D6E4}$ at a vertical coordinate $x_{2}=y_{s}$ , which divides the fluid domain into a free-fluid region and a porous region. The accuracy of the final interface condition does not depend on $y_{s}$ (up to certain limits) as proven by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012). We confirm this statement numerically in the § 5.3.

We separate the flow above the interface (domain $\unicode[STIX]{x1D6FA}_{f}$ ) into a fast flow $(U_{i},P)$ and a perturbation $(u_{i}^{+},p^{+})$ ,

(3.5) $$\begin{eqnarray}u_{i}=U_{i}+u_{i}^{+},\quad p=P+p^{+}.\end{eqnarray}$$

The terms $(u_{i}^{+},p^{+})$ are generated by the porous medium and will – as shown below – be responsible for the induced slip and penetration velocities. The pressure and the velocity below the interface (domain $\unicode[STIX]{x1D6FA}_{p}$ ), denoted by

(3.6) $$\begin{eqnarray}u_{i}=u_{i}^{-},\quad p=p^{-},\end{eqnarray}$$

represent the slow flow and the pressure field in the pores. Table 1 summarizes the introduced quantities in $\unicode[STIX]{x1D6FA}_{f}$ and $\unicode[STIX]{x1D6FA}_{p}$ . By inserting the decomposition (3.5) and the quantities (3.6) into (3.3)–(3.4) and grouping the different terms, the equations governing the dynamics of the different quantities are obtained.

Table 1.

List of the defined quantities and their properties. For each dimensionless quantity we provide the corresponding dimensional scale, domain, illustrate if microscale and macroscale variations are present and also state the dimensionless order. Appendix A provides more details and discussion on the scaling.

3.2.1 Fast flow

The global pressure $P$ and the fast flow $U_{i}$ are governed by the Navier–Stokes equations with no-slip condition at $y_{s}$ , i.e.

(3.7) $$\begin{eqnarray}\displaystyle {\mathcal{A}}(U_{i},P,\unicode[STIX]{x1D716})= & \!\unicode[STIX]{x1D716}^{2}Re_{d}(\unicode[STIX]{x2202}_{t}U_{i}+U_{j}U_{i,j}),\quad & \displaystyle y\geqslant y_{s},\end{eqnarray}$$

(3.8) $$\begin{eqnarray}\displaystyle U_{i}= & \!0,\quad & \displaystyle y=y_{s}.\end{eqnarray}$$

Here, ${\mathcal{A}}$ is a linear Stokes operator defined by (B 1) in appendix B, where we have summarized definitions of various operators and tensors. The macroscale pressure associated with fast flow is assumed to have a magnitude of $\tilde{P}\sim \unicode[STIX]{x0394}P$ . The fast velocity $U_{i}$ has a characteristic free-flow velocity, $\tilde{U} _{i}\sim U^{f}$ , that is ‘faster’ than the Darcy velocity scale $U^{d}$ ; we assume

(3.9) $$\begin{eqnarray}U^{f}\sim \unicode[STIX]{x1D716}^{-1}U^{d}.\end{eqnarray}$$

Using the introduced non-dimensionalization (3.2) on the variable estimates, we arrive at a priori orders of the fast flow velocity and the global pressure

(3.10) $$\begin{eqnarray}U_{i}=O(\unicode[STIX]{x1D716}^{-1}),\quad P=O(1).\end{eqnarray}$$

The dimensionless fast flow field $U_{i}$ becomes very large for small $\unicode[STIX]{x1D716}$ , whereas the global pressure in $\unicode[STIX]{x1D6FA}_{f}$ – either externally imposed as in pressure-driven channel flow or induced by the fast flow as in lid-driven cavity – is of order one.

Note that the assumption (3.9) is not a universal one, and depends on the bulk Reynolds number; for example, for Stokes flow $U^{f}\sim \unicode[STIX]{x1D716}^{-2}U^{d}$ . Appendix A shows that (3.9) is obtained, when $Re_{f}=\unicode[STIX]{x1D70C}_{f}HU^{f}/\unicode[STIX]{x1D707}\sim \unicode[STIX]{x1D716}^{-1}$ . This a priori scale estimate simplifies the multiscale expansion outlined in § 4 and its consequence is of a theoretical nature; in practice, our derived boundary condition predicts very accurately the slip and penetration velocity for a wide range of parameters, as we demonstrate in § 5.3.

3.2.2 Perturbations and slow flow

The perturbations above the interface are governed by

(3.11) $$\begin{eqnarray}\displaystyle {\mathcal{A}}(u_{i}^{+},p^{+},\unicode[STIX]{x1D716})= & \!\unicode[STIX]{x1D716}^{2}Re_{d}(\unicode[STIX]{x2202}_{t}u_{i}^{+}+u_{j}^{+}u_{i,j}^{+}+U_{j}u_{i,j}^{+}+u_{j}^{+}U_{i,j}),\quad & \displaystyle y\geqslant y_{s},\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle u_{i}^{+}= & \!u_{i}^{-},\quad & \displaystyle y=y_{s}.\end{eqnarray}$$

In order to solve these equations, one has to know the flow below the interface ( $u_{i}^{-}$ ). The pressure and the slow velocity below the interface are governed by

(3.13) $$\begin{eqnarray}\displaystyle {\mathcal{A}}(u_{i}^{-},p^{-},\unicode[STIX]{x1D716})= & \!\unicode[STIX]{x1D716}^{2}Re_{d}(\unicode[STIX]{x2202}_{t}u_{i}^{-}+u_{j}^{-}u_{i,j}^{-}),\quad & \displaystyle y\leqslant y_{s},\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F4}_{ij}^{u^{-}}n_{j}= & \!\unicode[STIX]{x1D6F4}_{ij}^{u^{+}}n_{j}+\unicode[STIX]{x1D6F4}_{ij}^{U}n_{j},\quad & \displaystyle y=y_{s},\end{eqnarray}$$

where the stress tensors containing $u_{i}^{+}$ , $u_{i}^{-}$ and $U_{i}$ in (3.14) are defined by (B 7)–(B 9) in appendix B. Note that the decomposition of the flow introduced in § 3.2 is exact, since continuity of both velocity (3.12) and total stress (3.14) is imposed at the interior interface $\unicode[STIX]{x1D6E4}$ . In other words, if one would sum up (3.7), (3.11), (3.13) and the boundary conditions (3.8), (3.12), (3.14), the Navier–Stokes equations (3.4) defined in the full domain would be recovered.

The perturbation terms $(u_{i}^{+},p^{+})$ in $\unicode[STIX]{x1D6FA}_{f}$ are an effect of the porous medium. Therefore the perturbation velocity is estimated by the seepage velocity, in dimensional setting $\tilde{u} _{i}^{+}\sim U^{d}$ , and the pressure perturbation by the microscale pressure $\tilde{p}^{+}\sim \unicode[STIX]{x0394}p$ , where $\unicode[STIX]{x0394}p$ is the pressure induced by the seepage velocity, see appendix A. In the current non-dimensional setting, we have

(3.15) $$\begin{eqnarray}u_{i}^{+}=O(1),\quad p^{+}=O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

Comparing to (3.10), we observe $\Vert p^{+}\Vert \ll \Vert P\Vert$ and $\Vert u_{i}^{+}\Vert \ll \Vert U_{i}\Vert$ when $\unicode[STIX]{x1D716}\ll 1$ .

The slow velocity below the interface can be estimated as $\tilde{u} _{i}^{-}\sim U^{d}$ , since it is reasonable to assume that the flow inside the porous medium has the magnitude of the Darcy velocity. The pore pressure, on the other hand, can be estimated as $\tilde{p}^{-}\sim \unicode[STIX]{x0394}P$ , because the global macroscale pressure is present also in the porous medium. For dimensionless variables we then have

(3.16) $$\begin{eqnarray}u_{i}^{-}=O(1),\quad p^{-}=O(1).\end{eqnarray}$$

The dimensional estimates and dimensionless orders of the decomposed quantities are summarized in table 1.

4.1 $O(1)$ equation and its analytical solution in an interface cell

Collecting the terms with prefactor $1$ , we get the following system

(4.7) $$\begin{eqnarray}\displaystyle p^{-(0)}n_{i}= & \!Pn_{i},\quad & \displaystyle y=y_{s},\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle p_{,i_{1}}^{-(0)}= & \!0,\quad & \displaystyle y\leqslant y_{s}.\end{eqnarray}$$

We observe that the zeroth-order pressure in the porous region $p^{-(0)}$ is independent of the microscale coordinate $x_{i}$ . For our purpose, which is to derive the macroscale effective boundary condition, it is sufficient to solve this equation in an elongated cell near the vicinity of the interface (figure 4 b). The size of this cell

(4.9) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}_{cell}=\{y_{1}\leqslant y\leqslant y_{2},-{\textstyle \frac{1}{2}}\leqslant x\leqslant {\textstyle \frac{1}{2}}\}\end{eqnarray}$$

is chosen in such a way to capture only the microscale behaviour near the interface. Recall that we use $x$ to describe the streamwise coordinate $x_{1}$ and $y$ to describe the interface-normal coordinate $x_{2}$ . The solution to (4.7)–(4.8) below the interface in $\unicode[STIX]{x1D6FA}_{cell}$ is constant and equal to global pressure $P$ at the interface,

(4.10) $$\begin{eqnarray}p^{-(0)}=P|_{y_{s}},\quad y_{1}\leqslant y\leqslant y_{s}.\end{eqnarray}$$

We will use this result in § 5.1 to derive a macroscale pressure condition at the interface.

4.2 $O(\unicode[STIX]{x1D716})$ equation and its computational solution in an interface cell

Next, we collect the first-order terms with prefactor $\unicode[STIX]{x1D716}$ , which results in the Stokes equations for $(u_{i}^{\pm (0)},p^{\pm (1)})$ . Specifically, above the interface, we have

(4.11) $$\begin{eqnarray}{\mathcal{A}}_{1}(u_{i}^{+(0)},p^{+(1)},1)=0,\quad y\geqslant y_{s},\end{eqnarray}$$

where ${\mathcal{A}}_{1}$ denotes the Stokes operator, which contains derivatives with respect to microscale $x_{i}$ , as defined in (B 4). Below the interface, we have

(4.12) $$\begin{eqnarray}{\mathcal{A}}_{1}(u_{i}^{-(0)},p^{-(1)},1)=p_{,i_{0}}^{-(0)}(X_{i}),\quad y\leqslant y_{s}.\end{eqnarray}$$

We observe that the slow flow $u_{i}^{-(0)}$ is forced by the macroscale gradient of the pore pressure term $(p_{,i_{0}}^{-(0)})$ . In contrast, above the interface the equation for perturbation $u_{i}^{+(0)}$ is not driven by a lower-order pressure term; the $O(1)$ pressure above the surface is $P$ , which is contained in equations for $U_{i}$  (3.7)–(3.8). The same global macroscale pressure is driving the fast flow and the slow flow, but whereas this pressure is defined as $P$ above the interface, it is obtained as the leading-order expansion term below the interface. This is a consequence of the fact that the fast flow is not expanded, since it varies only with macroscale coordinate by definition.

The boundary conditions of $O(\unicode[STIX]{x1D716})$ equations at the interface $y_{s}$ are continuity of velocity

(4.13) $$\begin{eqnarray}u_{i}^{-(0)}=u_{i}^{+(0)},\end{eqnarray}$$

and jump in stress

(4.14) $$\begin{eqnarray}\unicode[STIX]{x1D6F4}_{ij}^{u^{-}(1)}n_{j}=\unicode[STIX]{x1D6F4}_{ij}^{u^{+}(1)}n_{j}+\unicode[STIX]{x1D61A}_{ij}n_{j}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D6F4}_{ij}^{u^{\pm }(1)}=-p^{\pm (1)}\unicode[STIX]{x1D6FF}_{ij}+(u_{i,j_{1}}^{\pm (0)}+u_{j,i_{1}}^{\pm (0)})$ is the stress tensor of the perturbation velocity and the slow velocity, whereas $\unicode[STIX]{x1D61A}_{ij}=U_{i,j}+U_{j,i}$ is the strain tensor of the fast flow.

The system of equations (4.11)–(4.14) is solved in the interface cell $\unicode[STIX]{x1D6FA}_{cell}$ shown in figure 4(b). To complete the problem formulation, boundary conditions are needed at the sides of the cell. Due to regularity of the porous structure, we impose periodic conditions on the sides, as shown in figure 4(c). At the bottom of the cell, we impose the interior solution, which is

(4.15) $$\begin{eqnarray}u_{i}^{-(0)}=-\unicode[STIX]{x1D612}_{ij}^{\text{itr}}p_{,j_{0}}^{-(0)},\quad y=y_{1},\end{eqnarray}$$

where $\unicode[STIX]{x1D612}_{ij}^{\text{itr}}$ is the classical interior permeability field. For the derivation of this expression and the corresponding microscale problem, the reader is referred to the book by Mei & Vernescu (Reference Mei and Vernescu2010). From the literature (Mikelić & Jäger Reference Mikelić and Jäger2000; Jäger & Mikelić Reference Jäger and Mikelić2009; Marciniak-Czochra & Mikelić Reference Marciniak-Czochra and Mikelić2012; Carraro et al. Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013) it is well known that the interface cell is exposed to a zero-stress condition at the infinity within the method of matched asymptotic expansion, therefore we impose zero-stress condition, i.e. $\unicode[STIX]{x1D6F4}_{ij}^{u^{+}(1)}n_{j}=0$ , at $y=y_{2}$ .

We may consider the Stokes equations (4.11)–(4.12) and the boundary conditions (4.13)–(4.15) as one linear problem with four unknowns $(u_{i}^{\pm (0)},p^{\pm (1)})$ . Due to linearity, we can construct the solutions for the velocity and pressure fields as the superposition of $p_{,i_{0}}^{-(0)}(X_{i})$ and $\unicode[STIX]{x1D61A}_{ij}(X_{i})$ , i.e.

(4.16) $$\begin{eqnarray}u_{i}^{\pm (0)}=-\unicode[STIX]{x1D612}_{ij}^{\pm }p_{,j_{0}}^{-(0)}+\unicode[STIX]{x1D613}_{ijk}^{\pm }\unicode[STIX]{x1D61A}_{jk}|_{y_{s}},\end{eqnarray}$$

and

(4.17) $$\begin{eqnarray}p^{\pm (1)}=-\unicode[STIX]{x1D608}_{j}^{\pm }p_{,j_{0}}^{-(0)}+\unicode[STIX]{x1D609}_{ij}^{\pm }\unicode[STIX]{x1D61A}_{ij}|_{y_{s}}.\end{eqnarray}$$

The average of $\unicode[STIX]{x1D612}_{ij}$ is a tensorial effective Darcy permeability and the average of $\unicode[STIX]{x1D613}_{ijk}$ is related to the tensorial version of the Navier-slip coefficient. Similarly, the tensors $\unicode[STIX]{x1D608}_{j}$ and $\unicode[STIX]{x1D609}_{ij}$ are transfer coefficients from the driving pressure gradient and fluid strain, respectively, to the perturbation pressure.

Figure 5.

Solutions of interface problems for the coefficients of the Darcy term ( $\unicode[STIX]{x1D719}=0.15$ and $y_{s}=0.1$ ). The frames (a–d) from left to right correspond to the flow fields $\unicode[STIX]{x1D612}_{11},\unicode[STIX]{x1D612}_{21},\unicode[STIX]{x1D612}_{12}$ and $\unicode[STIX]{x1D612}_{22}$ . The arrows indicate the direction of the constant unit volume forcing below the interface (horizontal dashed line). Rightmost frame (e) shows plane-averaged profiles; the streamwise component provides the interface permeability ( $\langle \unicode[STIX]{x1D612}_{11}\rangle =0.014$ above $y_{s}$ ), whereas the wall-normal component is constant and corresponds to the interior permeability ( $\langle \unicode[STIX]{x1D612}_{22}\rangle =0.0266$ ).

4.2.1 Microscale Stokes problems for Darcy term

By inserting the ansatzes (4.16)–(4.17) into (4.11)–(4.15), it follows that the tensors $\unicode[STIX]{x1D612}_{ij}^{\pm }$ and $\unicode[STIX]{x1D608}_{i}^{\pm }$ satisfy,

(4.18) $$\begin{eqnarray}\displaystyle {\mathcal{A}}_{1}(\unicode[STIX]{x1D612}_{ik}^{+},\unicode[STIX]{x1D608}_{k}^{+},1)= & \!0,\quad & \displaystyle y\geqslant y_{s},\end{eqnarray}$$

(4.19) $$\begin{eqnarray}\displaystyle {\mathcal{A}}_{1}(\unicode[STIX]{x1D612}_{ik}^{-},\unicode[STIX]{x1D608}_{k}^{-},1)= & \!-\unicode[STIX]{x1D6FF}_{ik},\quad & \displaystyle y\leqslant y_{s},\end{eqnarray}$$

with boundary conditions at the interface $y_{s}$ given by

(4.20) $$\begin{eqnarray}\unicode[STIX]{x1D612}_{ik}^{-}=\unicode[STIX]{x1D612}_{ik}^{+},\quad \unicode[STIX]{x1D6F4}_{ij}^{K^{+}}n_{j}=\unicode[STIX]{x1D6F4}_{ij}^{K^{-}}n_{j}.\end{eqnarray}$$

At the bottom boundary $y=y_{1}$ , we have $\unicode[STIX]{x1D612}_{ik}^{-}=\unicode[STIX]{x1D612}_{ik}^{\text{itr}}$ (note that fields $\unicode[STIX]{x1D612}_{ik}^{-}$ and $\unicode[STIX]{x1D612}_{ik}^{\text{itr}}$ are not coinciding, but have the same values at $y=y_{1}$ ). The field $\unicode[STIX]{x1D612}_{jk}^{\pm }$ represents the $j$ th velocity component of the $k$ th Stokes problem. Thus to determine every component of $\unicode[STIX]{x1D612}_{ij}^{\pm }$ and $\unicode[STIX]{x1D608}_{i}^{\pm }$ , 3 pairs of Stokes problems have to be solved coupled at the interface through continuity of velocity field and stress. Note that below the interface, the flow is driven by a unit forcing in one direction at a time. Therefore the physical interpretation of $\unicode[STIX]{x1D612}_{i1}$ for example, is the flow response to forcing in the horizontal direction below the interface, as shown in figure 5(a,b).

To obtain reliable results, the interface cell needs to extend sufficiently into the free fluid such that variations of $\unicode[STIX]{x1D612}_{ij}^{+}$ are small and sufficiently into the porous medium such that variations of $\unicode[STIX]{x1D612}_{ij}^{-}$ are periodic. We have investigated different heights of the interface cell, and have determined that height of $10l$ (containing 5 cylinders below the interface) is sufficient. We have computed the $\unicode[STIX]{x1D612}_{ij}$ and $\unicode[STIX]{x1D613}_{ijk}$ fields for finer meshes and found no numerical artefacts nor noticeable modification of the results. Figure 5 shows $\unicode[STIX]{x1D612}_{ij}^{\pm }$ fields and corresponding plane-averaged profiles near the tip of the solid structure for interface location $y_{s}=0.1$ .

4.2.2 Microscale Stokes problems for Navier-slip term

By inserting the ansatzes (4.16)–(4.17) into (4.11)–(4.15), the following equations for the tensors $\unicode[STIX]{x1D613}_{ijk}^{\pm }$ and $\unicode[STIX]{x1D609}_{ij}^{\pm }$ are obtained,

(4.21) $$\begin{eqnarray}\displaystyle {\mathcal{A}}_{1}(\unicode[STIX]{x1D613}_{ikl}^{+},\unicode[STIX]{x1D609}_{kl}^{+},1)= & \!0,\quad & \displaystyle y\geqslant y_{s},\end{eqnarray}$$

(4.22) $$\begin{eqnarray}\displaystyle {\mathcal{A}}_{1}(\unicode[STIX]{x1D613}_{ikl}^{-},\unicode[STIX]{x1D609}_{kl}^{-},1)= & \!0,\quad & \displaystyle y\leqslant y_{s},\end{eqnarray}$$

with boundary conditions at $y_{s}$ given by,

(4.23) $$\begin{eqnarray}\unicode[STIX]{x1D613}_{ikl}^{-}=\unicode[STIX]{x1D613}_{ikl}^{+},\quad \unicode[STIX]{x1D6F4}_{ij}^{L^{+}}n_{j}=\unicode[STIX]{x1D6F4}_{ij}^{L^{-}}n_{j}-\unicode[STIX]{x1D6FF}_{ik}n_{l}.\end{eqnarray}$$

At the bottom boundary $y=y_{1}$ , we have $\unicode[STIX]{x1D613}_{ikl}^{-}=0$ . The tensors transferring the stress of the free fluid to the perturbation velocity require the solution of $9$ pairs of coupled Stokes problems. The forcing for these equations is at the interface in the form of a stress condition. For example, the $\unicode[STIX]{x1D613}_{i12}$ component is the flow response to a unit tangential stress at the interface, whereas $\unicode[STIX]{x1D613}_{i22}$ is the response to unit normal stress at the interface. In general, for problems with flat interfaces described in a coordinate system aligned with the interface, only three pairs of problems are forced.

Figure 6.

Solutions of interface problem for the non-zero component ( $\unicode[STIX]{x1D613}_{i12}$ ) of the Navier-slip term ( $\unicode[STIX]{x1D719}=0.15$ and $y_{s}=0.1$ ). The left and centre frames correspond to $\unicode[STIX]{x1D613}_{112}$ and $\unicode[STIX]{x1D613}_{212}$ , respectively. The arrows indicate the direction of the constant unit boundary forcing at the interface location. Rightmost frame (c) shows plane-averaged profiles; the streamwise component provides the interface slip length ( $\langle \unicode[STIX]{x1D613}_{112}\rangle =0.15$ above $y_{s}$ ), whereas the average of wall-normal component is zero (i.e. $\langle \unicode[STIX]{x1D613}_{212}\rangle =0$ ), since the two-dimensional field is antisymmetric with respect to the centre axis.

Returning to our two-dimensional configuration with a flat interface and aligned coordinate system, $\unicode[STIX]{x1D613}_{ij1}$ are unforced problems, leading to trivial solutions for all components. Out of the forced problems $\unicode[STIX]{x1D613}_{ij2}$ , the components $\unicode[STIX]{x1D613}_{122}$ and $\unicode[STIX]{x1D613}_{222}$ are zero since the forcing is in a constrained direction, i.e. due to mass conservation, the motion in the vertical direction is zero when the no-slip condition is enforced at the bottom boundary. We are thus left with only one non-trivial problem ( $\unicode[STIX]{x1D613}_{i12}$ ), for which the flow fields are shown in figure 6, along with corresponding plane-averaged profiles.

5.1 Condition for pore pressure

Using decomposition (3.5) and the scaling for pressure perturbation (3.15), we can write the pressure field above the interface as

(5.1) $$\begin{eqnarray}p=P+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

Taking the average of the expression above in a microscale volume of size $l^{3}$ – see definition (B 17) – above the interface gives,

(5.2) $$\begin{eqnarray}\langle p\rangle =P+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

Here, we have no angle brackets around $P$ , since it is independent of the microscale coordinate. The multiscale expanded pressure field (4.6) below the interface, on the other hand, is

(5.3) $$\begin{eqnarray}p^{-}=p^{-(0)}+O(\unicode[STIX]{x1D716}),\end{eqnarray}$$

where its volume-averaged form is

(5.4) $$\begin{eqnarray}\langle p^{-}\rangle =p^{-(0)}+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

According to solution of the $O(1)$ -problem (4.10), one can state that the pressure field in the whole interface cell is constant and equal to the macroscale pressure $P$ . Inserting the solution of the $O(1)$ -problem (4.10) in (5.4) and then equating it to expression (5.2), we obtain the following at the interface,

(5.5) $$\begin{eqnarray}\langle p^{-}\rangle =\langle p\rangle +O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

For brevity, we denote averaged dimensional quantities with a ‘hat’ (e.g. $\hat{p}^{-}=\unicode[STIX]{x0394}P\langle p^{-}\rangle$ ), which gives the pressure interface condition in its final dimensional form,

(5.6) $$\begin{eqnarray}\hat{p}^{-}=\hat{p}.\end{eqnarray}$$

Working with the chosen estimates (see table 1), one obtains pressure continuity up to $O(\unicode[STIX]{x1D716})$ for any anisotropic porous bed. We point out that from (4.17), one may formulate a pressure condition valid to $O(\unicode[STIX]{x1D716}^{2})$ . This is however out of the scope for this work. We note that this result is different compared to works by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) and Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013). This is a direct consequence of the theoretical assumption $Re_{f}\sim \unicode[STIX]{x1D716}^{-1}$ , which leads to the $O(1)$ -problem for pressure being trivial.

5.2 Velocity boundary condition for free fluid

The solution to the $O(\unicode[STIX]{x1D716})$ -problem (see § 4.2) is obtained numerically by computing $\unicode[STIX]{x1D612}_{ij}$ and $\unicode[STIX]{x1D613}_{ijk}$ . One may then proceed to construct the fully resolved flow field with error $O(\unicode[STIX]{x1D716})$ near the interface. First, we write the velocity above the interface as

(5.7) $$\begin{eqnarray}u_{i}^{+}=u_{i}^{+(0)}+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

The macroscale term $U_{i}$ does not appear in the expression above, because it is constant in the interface cell and this constant has be zero due to the boundary condition (3.8). Inserting the above expression into (4.16) gives

(5.8) $$\begin{eqnarray}u_{i}^{+}=-\unicode[STIX]{x1D612}_{ij}^{+}p_{,j_{0}}^{-(0)}+\unicode[STIX]{x1D613}_{ijk}^{+}(U_{j,k}|_{y_{s}}+U_{k,j}|_{y_{s}})+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

We average out the microscale oscillations by forming the volume average above the interface,

(5.9) $$\begin{eqnarray}\langle u_{i}^{+}\rangle =-\langle \unicode[STIX]{x1D612}_{ij}^{+}\rangle \,p_{,j_{0}}^{-(0)}+\langle \unicode[STIX]{x1D613}_{ijk}^{+}\rangle \,(U_{j,k}|_{y_{s}}+U_{k,j}|_{y_{s}})+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

Here, we have no angle brackets around pressure and velocity gradients since these quantities are independent of the microscale coordinate.

Now, by inserting the approximations of the velocity strain (C 3) and the pressure gradient (C 4) – see appendix C – into (5.9), we obtain

(5.10) $$\begin{eqnarray}\langle u_{i}\rangle =-{\mathcal{K}}_{ij}\frac{1}{\unicode[STIX]{x1D716}}\langle p^{-}\rangle _{,j}+{\mathcal{L}}_{ijk}(\langle u_{j}\rangle _{,k}+\langle u_{k}\rangle _{,j})+O(\unicode[STIX]{x1D716}),\end{eqnarray}$$

where for convenience, we have denoted ${\mathcal{K}}_{ij}=\langle \unicode[STIX]{x1D612}_{ij}^{+}\rangle$ and ${\mathcal{L}}_{ijk}=\langle \unicode[STIX]{x1D613}_{ijk}^{+}\rangle$ . The coefficients ${\mathcal{K}}_{ij}$ and ${\mathcal{L}}_{ijk}$ are evaluated over $l^{3}$ cube ( $l^{2}$ in 2D setting) at the top of the finite interface cell. Finally, we have used that $u_{i}=U_{i}+u_{i}^{+}=u_{i}^{+}$ at the interface.

In order to return to the boundary conditions (2.11)–(2.12) used for the lid-driven cavity problem in § 2, we revert the boundary condition to dimensional quantities,

(5.11) $$\begin{eqnarray}\hat{u} _{i}=-\left({\mathcal{K}}_{ij}\frac{l^{2}}{\unicode[STIX]{x1D707}}\right)\hat{p}_{,j}^{-}+({\mathcal{L}}_{ijk}l)(\hat{u} _{j,k}+\hat{u} _{k,j}).\end{eqnarray}$$

The coefficients, based on solutions of the interface cell problems, for the cavity flow are

(5.12) $$\begin{eqnarray}K_{cyl}=l^{2}\,{\mathcal{K}}_{22},\quad K_{s}=l^{2}\,{\mathcal{K}}_{11}\quad \text{and}\quad L_{s}=l\,{\mathcal{L}}_{112}.\end{eqnarray}$$

All other components of tensors $\unicode[STIX]{x1D612}_{ij}$ and $\unicode[STIX]{x1D613}_{ijk}$ are zero for the porous bed with regular circular cylinders, and therefore there is no velocity shear term for the penetration velocity (2.11). Actual values used in § 2 are ${\mathcal{K}}_{11}=0.0312$ , ${\mathcal{K}}_{22}=0.0986$ and ${\mathcal{L}}_{112}=0.1783$ .

Equation (5.11) is the final expression of the velocity boundary condition for a rigid porous bed. It can be used together with Navier–Stokes equations in any domain of interest, in order to take into account the effects of the porous medium, without resolving the microscale flow within the porous bed. We emphasize that the ‘minus’ notation for pressure means that the pressure gradient in the boundary condition is the gradient of the pore pressure. In the final subsection, we test the robustness of the derived boundary condition by varying solid volume fraction, scale separation parameter and interface location.

Table 2.

Cavity-slip velocity values at maximum predicted by model $\hat{u} _{s}$ normalized with DNS result $\bar{u}_{s}$ . The DNS results are plane averaged over a microscale length. Darcy contribution $\hat{u} _{sK}=-K_{s}/\unicode[STIX]{x1D707}\,\unicode[STIX]{x2202}_{x}\hat{p}^{-}$ and Navier-slip contribution $\hat{u} _{sL}=L_{s}(\unicode[STIX]{x2202}_{y}\hat{u} +\unicode[STIX]{x2202}_{x}\hat{v})$ are listed separately.

5.3 Accuracy of slip prediction and robustness to interface location

We now return to the example of the lid-driven cavity with a porous bed in order to illustrate more quantitatively that the proposed boundary condition yields accurate and robust slip velocity predictions. More specifically, we carry out a parametric study and report predictions of maximum slip velocity $\hat{u} _{s}$ at the centre of the cavity. In order to do a fair comparison, we surface average the DNS results $\bar{u}_{s}=\langle u_{s}\rangle _{S}$ at the interface. We also assess the contributions from two different terms in the derived boundary condition (5.11). Table 2 shows that for the range of parameters considered, the contribution at the interface from the Navier-slip term is always at least an order of magnitude larger than the contribution from the Darcy term. This is a consequence of (5.11), where $K_{s}\sim l^{2}$ and $L_{s}\sim l$ , and therefore $K_{s}\ll L_{s}$ for fine microstructures. This result is in agreement with previous work. It was first suggested by Saffman (Reference Saffman1971) that the Darcy term is of higher order and can be neglected, the same result was later rigorously proved by Mikelić & Jäger (Reference Mikelić and Jäger2000). One may therefore obtain a good approximation with only the Navier-slip term, as first suggested by Saffman (Reference Saffman1971) and later rigorously shown by Mikelić & Jäger (Reference Mikelić and Jäger2000). Including the Darcy term however yields – consistently – a smaller error, and therefore also a robust velocity boundary condition with respect to the interface location and different pore geometries. Additionally, the Darcy term is the only contribution appearing in the interface-normal direction, which is essential to capture the momentum transfer from and to the porous region.

Note that although the Darcy term is much smaller than the Navier-slip term, both terms appear in the $O(\unicode[STIX]{x1D716})$ equation (see § 4.2). This is a consequence of the estimate that perturbation velocity is of the same order as the Darcy velocity, i.e.  $u_{i}^{\pm }\sim U^{d}$ and that the fast flow velocity scales as (3.9). An alternative approach would be to assume $U^{f}\sim \unicode[STIX]{x1D716}^{-2}U^{d}$ , which would essentially result in three velocity scales to allow for the slip velocity to be faster than the Darcy flow but slower than the free flow. In such an approach – which would more sophisticated than the current one – the Darcy term can appear in higher-order equation than the equation for which the Navier-slip term appears. Nevertheless, one can observe that our direct numerical simulations are in good agreement with simulations of the homogenized model despite that we strictly speaking do not satisfy (3.9) for Stokes flow. Thus our method can be considered as practical parameter-free framework for computing the coefficients of a generalized condition proposed by Beavers and Joseph; a condition which has been employed under a variety of different flow conditions by experimentalists. To sum up, we have theoretically assumed that (i) $\unicode[STIX]{x1D716}\ll 1$ ; (ii) $Re_{d}\leqslant 1$ and (iii) $Re_{f}\sim \unicode[STIX]{x1D716}^{-1}$ . Based on numerical tests in this section for $Re_{f}=0$ and up to $\unicode[STIX]{x1D716}=0.1$ , we have determined the practical limits of the derived interface condition by relaxing conditions (i) and (iii), as summarized in § 2.2.