A framework for computing effective boundary conditions at the interface between free fluid and a porous medium

Abstract

Interfacial boundary conditions determined from empirical or ad hoc models remain the standard approach to model fluid flows over porous media, even in situations where the topology of the porous medium is known. We propose a non-empirical and accurate method to compute the effective boundary conditions at the interface between a porous surface and an overlying flow. Using a multiscale expansion (homogenization) approach, we derive a tensorial generalized version of the empirical condition suggested by Beavers & Joseph (J. Fluid Mech., vol. 30 (01), 1967, pp. 197–207). The components of the tensors determining the effective slip velocity at the interface are obtained by solving a set of Stokes equations in a small computational domain near the interface containing both free flow and porous medium. Using the lid-driven cavity flow with a porous bed, we demonstrate that the derived boundary condition is accurate and robust by comparing an effective model to direct numerical simulations. Finally, we provide an open source code that solves the microscale problems and computes the velocity boundary condition without free parameters over any porous bed.

1 Introduction

Surfaces found in nature are generally non-smooth with complex hierarchical structural features (Liu & Jiang 2011). The purposes of these surfaces vary greatly, ranging from camouflage and insulation to less obvious functions, such as passively interacting with surrounding fluid to reduce drag or noise (Abdulbari et al. 2013). These functions manifest as effective macroscale properties – for example, permeability, elasticity, slip and optical transparency – while their origin is the small-scale features of the surface. Therefore, to understand the hydrodynamic function of such complex surfaces, a systematic multiscale approach is required. In a bottom-up strategy, the microscale fluid–structure physics of the coating material is analysed first; the effective porosity, elasticity or slip are then induced naturally by upscaling the microscale features.

Volume-averaging and homogenization techniques (Davit et al. 2013) enable a bottom-up strategy by deriving the effective equations governing the macroscale coating dynamics, which contains parameters arising from microscale features. Whereas these techniques are routinely applied for homogeneous materials (e.g. the interior of a material), their application to inhomogeneous regions (e.g. near interfaces) has not reached the same level of maturity. One example, which is also the focus of the present work, is the interface between an overlying flow and a rigid porous surface. Recent works (Ochoa-Tapia & Whitaker 1995; Mikelić & Jäger 2000; Auriault 2010a ; Minale 2014) have treated the inhomogeneous interface problem theoretically with upscaling techniques. Ochoa-Tapia & Whitaker (1995) used a volume-averaging technique to derive a shear-stress jump condition. Later on, Valdés-Parada et al. (2013) used the same technique to analyse both stress and velocity jump across the interface. Interestingly, they identified a fixed location of the interface that yields best results when imposing a velocity jump. This is in contrast to both the theoretical findings by Marciniak-Czochra & Mikelić (2012) and the numerical results presented in this paper, which show that the accuracy of the velocity jump condition is independent of the interface location. Recently, Minale (2014) rederived the boundary conditions of Ochoa-Tapia & Whitaker (1995), elucidating how the stress from the free fluid is partitioned between the porous skeleton and the porous flow.

Volume-averaging techniques induce closure problems that need to be resolved using scale estimates. Homogenization techniques, on the other hand, begin with scale estimates and an expansion in small parameter $\unicode[STIX]{x1D716}=l/H$ , defining the scale separation between microscale $l$ and macroscale $H$ . With a homogenization approach, one obtains equations at different orders of $\unicode[STIX]{x1D716}$ and a decoupling of different quantities and thus also simpler closure problems. Mikelić & Jäger (2000) used homogenization and the method of matched asymptotic expansions to show that the Saffman (1971) version of the empirical boundary condition by Beavers & Joseph (1967) (called the BJ condition hereafter) is mathematically justified and its slip parameter can be computed by solving microscale problems in an interface unit cell. Auriault (2010a ) also used a homogenization technique to derive a BJ type of boundary condition valid for pressure-driven flows; he obtained however the condition at a different order compared to Mikelić & Jäger (2000), as seen in the discussions by Jäger & Mikelić (2010) and Auriault (2010b ). More recently, Carraro et al. (2015) repeated the procedure of Mikelić & Jäger (2000) to determine the boundary condition of the penetration (wall-normal) velocity component.

The boundary conditions derived using upscaling techniques have remained at a proof-of-concept level and only demonstrated on canonical one-dimensional flows. The theoretical progress has not yet resulted in a method that can, in a straightforward manner, be applied by practitioners and engineers. The reason is to some extent the non-trivial mathematical aspects – such as closure problems. Another reason is that the focus has been on mathematically justifying empirical boundary conditions, rather than presenting a step-by-step method for computing interfacial conditions. Therefore, investigations of practical interest of flows over porous media continue using empirical conditions or conditions with free unknown parameters (Han, Ganatos & Weinbaum 2005; Le Bars & Worster 2006; Rosti, Cortelezzi & Quadrio 2015; Zampogna & Bottaro 2016). Although these conditions provide physical models of the flow over porous media, they are based on lumping all unknown effects into few scalar parameters. This approach requires the support of empirical data (Zampogna & Bottaro 2016) or extensive computations to cover a large interval of parameters (Rosti et al. 2015).

In this work, we provide practitioners the framework to compute accurate interfacial velocity boundary conditions, instead of empirically determining them. We derive the interface boundary condition for slip velocity using homogenization and present the relevant Stokes equations to be solved in a microscale interface unit cell. Our main contribution is to provide a set of simple and numerically feasible microscale problems, which once solved, allows for a robust non-empirical effective interface condition. Our interface condition can be considered as a generalized version of the BJ condition, since it depends on interface permeability tensor and on the interface velocity strain rate tensor.

This paper is organised as follows. In § 2, using the lid-driven cavity with a porous bed as an example, we compare the velocity field computed from a direct numerical simulation to the field obtained by solving the homogenized (averaged) equations with the interface condition that is proposed in this paper. After that, in §§ 3–5, we derive the interface boundary condition. More specifically, in § 3 we decompose the physical domain into a porous part and free-fluid part and define an interface between the two domains. We then introduce the equations governing the microscale fluid flow in each part as well as their coupling through continuity of velocity and stress at the interface. In § 4, we use multiscale expansion to derive the relevant Stokes equations to be solved in a microscale interface unit cell in order determine the effective boundary condition. In § 5 we derive the interface conditions by employing a homogenization (averaging) technique and relate the obtained results back to the example presented in § 2. Finally, in § 6, we conclude the paper.

2 Direct numerical simulations versus continuum model

The purpose of this section is to compare two approaches to describe the flow in a lid-driven cavity with a homogenous bed of solid cylinders. In the first approach, represented in figure 1(a), we solve the Stokes equations over all spatial scales. This is possible for simplified geometries such as this one, but clearly for more complex (possibly three-dimensional) and much denser porous beds it is not feasible to solve Stokes equations with complete microscale resolution. In the second approach represented in figure 1(c), we reduce the degrees of freedom of the flow in the porous bed by homogenization.

Figure 1. Left frame (a) shows the lid-driven cavity domain with a porous bed. The top wall is driven with velocity $U_{w}$ . Centre frame (b) shows a magnified view of the porous bed. Cylinders are spaced apart by distance $l$ and $\unicode[STIX]{x1D6E4}_{c}$ defines the boundary of the cylinders. Right frame (c) shows a two-domain description of the same cavity problem in a homogenized setting. The parameters defining the problem are the volume fraction $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}r^{2}/l^{2}=0.02$ , scale separation $l/H=0.1$ and porous-bed depth of $d\approx 0.5H$ .

The cavity has a length $H$ and a depth of $(H+d)$ , where the porous bed is confined to $-d<y<0$ and $-H/2<x<H/2$ . Coordinate $y=0$ corresponds to the tangent plane of the top row of cylinders and the depth of the porous bed is $d\approx H/2$ . The top wall of the cavity is driven by a constant streamwise velocity $U_{w}$ , which is sufficiently slow to render fluid inertia negligible inside the cavity. Figure 1(b) shows the geometry of the porous bed, which consists of a lattice of cylinders with diameter $D=2r$ and with the spacing $l$ . For the particular example discussed in this section, the microscale length is $l/H=0.1$ and the cylinder volume fraction is $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}r^{2}/l^{2}=0.02$ .

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 2012), 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. (2015) recently showed that slip velocities below 3 % had a significant effect on the flow statistics in a turbulent channel. Additionally, Carotenuto & Minale (2013) 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ć (2012) and Carraro et al. (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 1973), but to the authors’ knowledge, it has not been derived earlier for flat interfaces, although it has be conjectured to exist by Nield (2009). 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 Governing equations and flow decomposition

While the derived effective boundary condition in this paper has been applied to the steady cavity flow over regular array of cylinders, the boundary condition is more general and can be applied to more complex three-dimensional flows, as schematically shown in figure 4(a). We relax the assumptions of a steady and non-inertial flow in the cavity, and start from incompressible Navier–Stokes equations. The length $H$ now corresponds to an appropriate macroscopic length scale of the flow, satisfying $\unicode[STIX]{x1D716}=l/H\ll 1$ .

Figure 4. Left frame (a) shows a schematic of a flow over porous bed with regular cylinders, where an interface $\unicode[STIX]{x1D6E4}$ between free fluid and porous bed has been introduced. The dashed rectangular domain corresponds to the interface cell used to compute the effective macroscale boundary condition. The coordinates and the boundary conditions imposed for solving cell problems are shown in frames (b) and (c), respectively.

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ć (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 Multiscale expansion

We now turn our attention to the multiscale analysis of the flow near the interface and construct an approximate description of it within an interface cell (defined below). To carry out the multiscale expansion, we introduce the macroscale and microscale coordinates

(4.1) $$\begin{eqnarray}X_{i}=\frac{\tilde{x}_{i}}{H}\quad \text{and}\quad x_{i}=\frac{\tilde{x}_{i}}{l},\end{eqnarray}$$

respectively. These coordinates are appropriate to describe the macroscopic and microscopic variations and are related to each other by $X_{i}=\unicode[STIX]{x1D716}x_{i}$ . In the new coordinates, there are two derivatives appearing due to the chain rule

(4.2) $$\begin{eqnarray}()_{,i}=()_{,i_{1}}+\unicode[STIX]{x1D716}()_{,i_{0}},\end{eqnarray}$$

where $()_{,i_{0}}$ denotes the derivative with respect to $X_{i}$ and $()_{,i_{1}}$ with respect to $x_{i}$ .

The fast flow $U_{i}$ and the global pressure $P$ do not depend on microscale coordinate, i.e. $U_{i,j_{1}}=0$ and $P_{,i_{1}}=0$ . This is a direct consequence of the definition of fast flow problem (3.7)–(3.8) and is valid for $\unicode[STIX]{x1D716}\ll 1$ . For $u_{i}^{\pm }$ and $p^{\pm }$ , which depend on both coordinates, we carry out the multiscale expansion as explained by Mei & Vernescu (2010). The perturbation velocity and the pressure above the interface ( $y\geqslant y_{s}$ ) are expanded as

(4.3) $$\begin{eqnarray}\displaystyle u^{+}(X_{i},x_{i}) & = & \displaystyle u_{i}^{+(0)}(X_{i},x_{i})+\unicode[STIX]{x1D716}u_{i}^{+(1)}(X_{i},x_{i})+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(4.4) $$\begin{eqnarray}\displaystyle p^{+}(X_{i},x_{i}) & = & \displaystyle \unicode[STIX]{x1D716}p^{+(1)}(X_{i},x_{i})+\unicode[STIX]{x1D716}^{2}p^{+(2)}(X_{i},x_{i})+O(\unicode[STIX]{x1D716}^{3}).\end{eqnarray}$$

The pressure expansion starts with the $O(\unicode[STIX]{x1D716})$ term, since $p^{+}=O(\unicode[STIX]{x1D716})$ . Below the interface ( $y\leqslant y_{s}$ ) the slow flow and the pore pressure are expanded as,

(4.5) $$\begin{eqnarray}\displaystyle u^{-}(X_{i},x_{i}) & = & \displaystyle u_{i}^{-(0)}(X_{i},x_{i})+\unicode[STIX]{x1D716}u_{i}^{-(1)}(X_{i},x_{i})+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle p^{-}(X_{i},x_{i}) & = & \displaystyle p^{-(0)}(X_{i},x_{i})+\unicode[STIX]{x1D716}p^{-(1)}(X_{i},x_{i})+O(\unicode[STIX]{x1D716}^{2}).\end{eqnarray}$$

We insert expansions (4.3)–(4.6) into the corresponding equations (3.11)–(3.14), and collect the terms at different orders. In the following subsections, we introduce and solve equation systems appearing at first two orders ( $O(1)$ and $O(\unicode[STIX]{x1D716})$ ).

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 (2010). From the literature (Mikelić & Jäger 2000; Jäger & Mikelić 2009; Marciniak-Czochra & Mikelić 2012; Carraro et al. 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 Effective interface conditions

This section provides the final forms of the effective boundary conditions by averaging the microscale solutions provided in the previous section. Since the interface cell is located at the boundary between the free-fluid and porous regions, the conditions for the free fluid should be evaluated using the solution above the interface, while conditions for the porous region should be evaluated using the solution below the interface. In particular, to solve for the pore pressure Laplace equation (2.10) below the interface, one needs a boundary condition for the pressure at the interface; this can be obtained by averaging the $O(1)$ -problem (see § 4.1). To solve the velocity of the free flow above the interface, one needs a condition for the velocity at the interface. We recall however that the $O(1)$ -problem given by (4.7)–(4.8) does not contain velocity. Therefore, one has to investigate solution of the $O(\unicode[STIX]{x1D716})$ -problem (4.16) to determine a boundary condition for the velocity. This is a consequence of the pore velocity viscous term being higher order compared to the pressure gradient term, see $\unicode[STIX]{x1D716}$ prefactors in (3.3).

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ć (2012) and Carraro et al. (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 (1971) that the Darcy term is of higher order and can be neglected, the same result was later rigorously proved by Mikelić & Jäger (2000). One may therefore obtain a good approximation with only the Navier-slip term, as first suggested by Saffman (1971) and later rigorously shown by Mikelić & Jäger (2000). 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.

6 Conclusions

We have presented a framework to construct a reduced homogenized model of the flow above and through a porous medium consisting of regular solid structures of general shape. The main contribution of the present paper is to provide the foundation and the tools to compute effective boundary conditions completely free from data fitting. The approach that we adopted can be summarized by the following four steps. First, the governing equations are made non-dimensional using scale estimates arising from flow in the porous domain. Second, the governing equations describing the fully resolved flow are separated at a virtual interface, and decomposed above the interface into equations for the fast flow and for perturbations. Third, multiscale expansion according to Mei & Vernescu (2010) is employed on perturbation and pore equations. Finally and after solving $O(1)$ - and $O(\unicode[STIX]{x1D716})$ -problems, we construct the interface conditions for pore pressure and fluid velocity using volume averages.

This procedure results in macroscopic description of the flow over a porous domain with an error $O(\unicode[STIX]{x1D716})$ . Specifically, using our a priori scaling estimates, the leading-order conditions are a pressure continuity condition and a generalized tensorial BJ boundary condition. The proposed velocity condition depends on the interface permeability and on the velocity strain, while the BJ condition contains interior permeability and velocity derivative of one component in one direction only. To the authors’ knowledge, such a general formulation has not been derived and validated before. Moreover, in order to obtain the constants of the effective boundary conditions, we derive a number of Stokes problems that need be solved numerically in small interface unit cells. Solvers for microscale problems have been released as an open source software (Lācis & Bagheri 2016), along with the solver used for the lid-driven cavity flow.

This work is also among the first to validate non-empirical boundary conditions on two-dimensional flows with DNS, where penetration velocity, slip velocity and pressure conditions have to be predicted to solve the coupled two-domain problem. The present boundary condition has been tested in the lid-driven cavity flow for a range of volume fractions from $\unicode[STIX]{x1D719}=0.02$ to $0.45$ , scale separation parameters from $\unicode[STIX]{x1D716}=0.02$ to $0.1$ and interface locations from $y_{s}=0.1$ to $y_{s}=0.5$ . When the homogenized model results are compared to DNS, the slip velocity predictions have been found to be robust and to give accurate predictions for all investigated parameters. We hope that, with this work and the release of the associated software, we can provide the numerical fluid dynamics community the tools to model flows over existing non-smooth surfaces as well as to design surfaces to modify fluid flow characteristics.