2.1 Tangential interface velocity and slip length

The tangential velocity condition in the TR model is provided by the standard slip condition, which for 3-D textured surfaces is

(2.1)$$\begin{eqnarray}(u_{x},u_{y})=\boldsymbol{u}_{t}=\frac{\unicode[STIX]{x1D647}}{\unicode[STIX]{x1D707}}\boldsymbol{\cdot }\unicode[STIX]{x1D749}\quad \text{on }z=z_{i},\end{eqnarray}$$

where $\unicode[STIX]{x1D749}=\unicode[STIX]{x1D707}(\unicode[STIX]{x2202}_{z}u_{x}+\unicode[STIX]{x2202}_{x}u_{z},\unicode[STIX]{x2202}_{z}u_{y}+\unicode[STIX]{x2202}_{y}u_{z})$ and $\unicode[STIX]{x1D647}=(\unicode[STIX]{x1D613}_{xx},\unicode[STIX]{x1D613}_{xy};\unicode[STIX]{x1D613}_{xy},\unicode[STIX]{x1D613}_{yy})$ is the symmetric positive definite (Kamrin & Stone Reference Kamrin and Stone2011) surface slip length tensor. Here, $\boldsymbol{u}_{t}$ is the tangential velocity vector, and the subscript $t$ is used interchangeably with the $x$ and $y$ components.

Figure 3.

(a) Shows a flow domain with a generic free flow. The flat interface with the vertical coordinate $z_{i}$ above the surface texture is depicted using a transparent plane. The red rectangle is the interface cell. (b) Shows the interface cell with a bottom coordinate $\hat{z}_{b}(x,y)$ – describing the surface texture – and a top coordinate $\hat{z}_{t}$. The tangential shear stress is decomposed in unit forcing terms along the $x$ (c) and the $y$ (d) axes.

We consider a patterned wall and a vortical flow over it, as illustrated in figure 3(a). The scale separation assumption (A2) allows us to introduce two different spatial coordinates: $x_{i}$ and $\hat{x}_{i}$. The former is used to describe spatial variations over large length scales ($x_{i}\sim H$). The latter is used to describe microscopic variations over much smaller roughness scale ($\hat{x}_{i}\sim \ell$). The effective boundary condition (2.1) is a macroscopic condition; the microscopic features of the texture are embedded in an averaged sense in the slip length tensor $\unicode[STIX]{x1D647}$.

To determine $\unicode[STIX]{x1D647}$, we consider a small volume near the surface of the texture with a cross-section $\ell \times \ell$. This volume contains a representative surface structure, see figure 3(a,b). Within this interface cell, the scale separation assumption (A2) allows us to treat the shear stress from the free fluid $\unicode[STIX]{x1D749}$ as a spatially constant external parameter. Due to the creeping flow assumption (A1), the equations governing the flow response to the free-fluid shear stress are the Stokes equations,

(2.2)$$\begin{eqnarray}\displaystyle \displaystyle -\unicode[STIX]{x1D735}\hat{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}} & = & \displaystyle -\unicode[STIX]{x1D6FF}(\hat{z}-\hat{z}_{i})\unicode[STIX]{x1D749},\end{eqnarray}$$

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

where $\unicode[STIX]{x1D6FF}$ is the Dirac delta function with $\hat{z}$ as an argument and consequently it has the unit of inverse metres (m$^{-1}$). This set of equations is equivalent to a two-domain description satisfying velocity continuity and stress jump at the interface (appendix B). Additionally, equations (2.2)–(2.3) are the same as previously used and derived by Luchini, Manzo & Pozzi (Reference Luchini, Manzo and Pozzi1991), Kamrin, Bazant & Stone (Reference Kamrin, Bazant and Stone2010) and Luchini (Reference Luchini2013). The imposed boundary conditions are no slip and no penetration on the solid structure ($\boldsymbol{u}=0$ on $\hat{z}=\hat{z}_{b}$). We impose periodic conditions at the vertical faces of the interface cell (due to the assumption (A3)). At the top surface of the cell, we impose a zero-stress condition ($\unicode[STIX]{x1D72E}\boldsymbol{\cdot }\boldsymbol{n}=0$ on $\hat{z}=\hat{z}_{t}$) to keep the shear stress at the interface as the only driving force of the problem.

The linearity assumption (A1) allows us to write the solution as a product between a response operator $\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}$ and the free-fluid shear stress,

(2.4)$$\begin{eqnarray}\hat{\boldsymbol{u}}=\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\unicode[STIX]{x1D749}.\end{eqnarray}$$

This can be expanded as

(2.5a)$$\begin{eqnarray}\displaystyle \hat{\boldsymbol{u}} & = & \displaystyle (\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{x})\unicode[STIX]{x1D70F}_{x}+(\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{y})\unicode[STIX]{x1D70F}_{y}\end{eqnarray}$$

(2.5b)$$\begin{eqnarray}\displaystyle & = & \displaystyle \hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}\frac{\unicode[STIX]{x1D70F}_{x}}{\unicode[STIX]{x1D707}}+\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D707}},\end{eqnarray}$$

(2.6a,b)$$\begin{eqnarray}\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}=\unicode[STIX]{x1D707}\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{x}\quad \text{and}\quad \hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}=\unicode[STIX]{x1D707}\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{y}.\end{eqnarray}$$

Thus the velocity fields $\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}$ and $\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}$ are solutions of the following two fundamental problems,

(F P s x)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\hat{p}^{(\unicode[STIX]{x1D70F}x)}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}=-\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FF}(\hat{z}-\hat{z}_{i})\boldsymbol{e}_{x},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}=0;\end{eqnarray}$$

(F P s y)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\hat{p}^{(\unicode[STIX]{x1D70F}y)}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}=-\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FF}(\hat{z}-\hat{z}_{i})\boldsymbol{e}_{y},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}=0.\end{eqnarray}$$

The fundamental problems are forced in the $x$ and $y$ directions, respectively, with unit shear at the plane $z_{i}$. The boundary conditions are the same as for equations (2.2)–(2.3). The fundamental problems are named with two capital letters FP complemented with two small letters denoting the forcing direction. For example, equation (FPsx) reads ‘Fundamental Problem forced by shear in $\boldsymbol{x}$-direction’.

Taking the surface average of expression (2.5b) at the interface, we obtain

(2.7)$$\begin{eqnarray}\boldsymbol{u}=\langle \hat{\boldsymbol{u}}\rangle _{i}=\langle \hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}\rangle _{i}\frac{\unicode[STIX]{x1D70F}_{x}}{\unicode[STIX]{x1D707}}+\langle \hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}\rangle _{i}\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D707}}\quad \text{on }z=z_{i}.\end{eqnarray}$$

No average is carried out for the free-fluid shear stress, because it is constant within the interface cell (A2). The surface average of an arbitrary quantity $\hat{a}$ is defined as

(2.8)$$\begin{eqnarray}\langle \hat{a}\rangle _{i}=\langle \hat{a}\rangle _{s}(\hat{z}_{i}),\quad \text{where }\langle \hat{a}\rangle _{s}(\hat{z})=\frac{1}{\ell ^{2}}\int _{0}^{\ell }\int _{0}^{\ell }\hat{a}(\hat{x},{\hat{y}},\hat{z})\,\text{d}\hat{x}\,\text{d}{\hat{y}}.\end{eqnarray}$$

By comparing the surface-averaged velocity in the interface cell (2.7) with the slip boundary condition (2.1), we observe that the components of the slip length tensor can be obtained as

(2.9a-d)$$\begin{eqnarray}\unicode[STIX]{x1D613}_{xx}=\langle \hat{u} _{x}^{(\unicode[STIX]{x1D70F}x)}\rangle _{i},\quad \unicode[STIX]{x1D613}_{yx}=\langle \hat{u} _{y}^{(\unicode[STIX]{x1D70F}x)}\rangle _{i},\quad \unicode[STIX]{x1D613}_{xy}=\langle \hat{u} _{x}^{(\unicode[STIX]{x1D70F}y)}\rangle _{i},\quad \unicode[STIX]{x1D613}_{yy}=\langle \hat{u} _{y}^{(\unicode[STIX]{x1D70F}y)}\rangle _{i}.\end{eqnarray}$$

In terms of the response operator, the slip tensor becomes $\unicode[STIX]{x1D647}=\unicode[STIX]{x1D707}\langle \hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\rangle _{i}$. Note that the dimension of the vector fields $\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}x)}$ and $\hat{\boldsymbol{u}}^{(\unicode[STIX]{x1D70F}y)}$ is the same as for velocity per shear, which gives a unit of metres (m). The units of the different fields introduced throughout the paper are summarised in table 2.

2.2 Interface normal velocity and transpiration length

We begin with a simple motivation for the transpiration velocity based on the principle of mass conservation. We consider a two-dimensional rough surface and define a control volume (CV) below the interface $z=z_{i}$ as shown in figure 4. We assume that there is a slip velocity variation from $u_{x,1}$ at the left side of the CV to $u_{x,2}$ at the right side of the CV. The mass fluxes at the left and the right boundaries of the CV are proportional to the slip velocities at the interface (a direct consequence of (A1)). Consequently, mass conservation requires a non-zero transpiration velocity $u_{z}$ at the interface. If the slip velocity is increasing $u_{x,2}>u_{x,1}$, the generated transpiration velocity is, therefore, negative.

Figure 4.

Illustration of the transpiration velocity as a consequence from mass conservation owing to the variation of slip velocity along the interface. Control volume (CV) below the interface ($\hat{z}=\hat{z}_{i}$) is denoted with a shaded (blue) region. The quantities $u_{x,1}$ and $u_{x,2}$ denote slip velocities at the left and the right boundaries of the CV, respectively.

More generally, the interface normal velocity condition in the TR model for a 3-D textured surface is provided by a linear law relating the normal velocity with the tangential variation of slip velocity,

(2.10)$$\begin{eqnarray}u_{n}=u_{z}=-\unicode[STIX]{x1D648}\,\boldsymbol{ : }\,\unicode[STIX]{x1D735}_{2}\boldsymbol{u}_{t}\quad \text{on }z=z_{i},\end{eqnarray}$$

where $\unicode[STIX]{x1D648}=(\unicode[STIX]{x1D614}_{xx},\unicode[STIX]{x1D614}_{xy};\unicode[STIX]{x1D614}_{xy},\unicode[STIX]{x1D614}_{yy})$ is the transpiration length tensor – exhibiting the same symmetry properties as the slip length tensor – and $\unicode[STIX]{x1D735}_{2}=(\unicode[STIX]{x2202}_{x},\unicode[STIX]{x2202}_{y})$ is gradient operator containing the two tangential directions. The double dot product between two second-rank tensors $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ is defined as $\unicode[STIX]{x1D63C}\boldsymbol{ : }\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D608}_{ij}\unicode[STIX]{x1D609}_{ij}$, where summation over repeating indices is implied. We use the normal $n$ subscript interchangeably with the $z$ component. The proposed expression, motivated by the conservation of mass, also emerges from a formal multi-scale expansion (see § 5, appendix A and Sudhakar et al. Reference Sudhakar, Lācis, Pasche and Bagheri2019). A similar multi-scale expansion has been recently used by Bottaro (Reference Bottaro2019) to confirm the present transpiration velocity condition.

To determine $\unicode[STIX]{x1D648}$, we make use of mass conservation in a 3-D setting. We define a CV with size $(x_{2}-x_{1})\times (y_{2}-y_{1})\times (\hat{z}_{i}-\hat{z}_{b})$ over a number of texture elements as shown in figure 5. By definition, no flux goes through the impermeable bottom surface. Therefore mass conservation requires

(2.11)$$\begin{eqnarray}Q_{1}+Q_{2}+Q_{3}+Q_{4}+Q_{5}=0,\end{eqnarray}$$

where $Q_{i}$ is the volumetric flux through faces of the CV (figure 5). The flux through the vertical faces ($i=1,\ldots ,4$) of the CV can be evaluated as

(2.12)$$\begin{eqnarray}Q_{i}=\int _{S_{i}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S=\int _{\hat{z}_{b}}^{\hat{z}_{i}}\int _{s_{1}}^{s_{2}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\hat{z}\,\text{d}s,\end{eqnarray}$$

where $\boldsymbol{u}$ is the effective velocity field at the CV face, $\boldsymbol{n}$ is the unit normal vector of the surface and $s$ is either $x$ or $y$, depending on the orientation of the surface. Note that the integral in the wall-normal direction is carried out over the microscale $\hat{z}$, because macroscopically the textured surface is infinitesimal and variation in depth does not exist.

Figure 5.

Control volume for deriving the transpiration length tensor $\unicode[STIX]{x1D648}$. All the possible volumetric fluxes are indicated with thick arrows.

Next, we use equation (2.8) in conjunction with the solution of fundamental problems (FPsx), (FPsy) to rewrite equation (2.12) as

(2.13)$$\begin{eqnarray}\displaystyle Q_{i} & = & \displaystyle \int _{\hat{z}_{b}}^{\hat{z}_{i}}\int _{s_{1}}^{s_{2}}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\hat{z}\,\text{d}s=\int _{\hat{z}_{b}}^{\hat{z}_{i}}\int _{s_{1}}^{s_{2}}\langle \hat{\boldsymbol{u}}\rangle _{s}(x,y,\hat{z})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\hat{z}\,\text{d}s\nonumber\\ \displaystyle & = & \displaystyle \int _{\hat{z}_{b}}^{\hat{z}_{i}}\int _{s_{1}}^{s_{2}}[\langle \hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\rangle _{s}(\hat{z})\boldsymbol{\cdot }\unicode[STIX]{x1D749}(x,y)]\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}\hat{z}\,\text{d}s=\left(\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\int _{s_{1}}^{s_{2}}\unicode[STIX]{x1D749}\,\text{d}s\right)\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$

where we have defined the response tensor $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}$ as

(2.14)$$\begin{eqnarray}\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}=\int _{\hat{z}_{b}}^{\hat{z}_{i}}\langle \hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}\rangle _{s}\,\text{d}\hat{z}=\frac{1}{\ell ^{2}}\int _{0}^{\ell }\int _{0}^{\ell }\int _{z_{b}}^{z_{i}}\hat{\unicode[STIX]{x1D64D}}_{\unicode[STIX]{x1D70F}}(\hat{x},{\hat{y}},\hat{z})\,\text{d}\hat{x}\,\text{d}{\hat{y}}\,\text{d}\hat{z}.\end{eqnarray}$$

The flux through the top wall is expressed as

(2.15)$$\begin{eqnarray}Q_{5}=\int _{x_{1}}^{x_{2}}\int _{y_{1}}^{y_{2}}u_{z}\,\text{d}x\,\text{d}y.\end{eqnarray}$$

Inserting the expressions for the fluxes through the CV faces into the mass conservation identity (2.11) we obtain

(2.16)$$\begin{eqnarray}\displaystyle \int _{x_{1}}^{x_{2}}\int _{y_{1}}^{y_{2}}u_{z}\,\text{d}x\,\text{d}y & = & \displaystyle -\left(\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\int _{y_{1}}^{y_{2}}[\unicode[STIX]{x1D749}(x_{2},y)-\unicode[STIX]{x1D749}(x_{1},y)]\,\text{d}y\right)\boldsymbol{\cdot }\boldsymbol{e}_{x}\nonumber\\ \displaystyle & & \displaystyle -\left(\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\int _{x_{1}}^{x_{2}}[\unicode[STIX]{x1D749}(x,y_{2})-\unicode[STIX]{x1D749}(x,y_{1})]\,\text{d}x\right)\boldsymbol{\cdot }\boldsymbol{e}_{y}.\end{eqnarray}$$

To proceed towards the effective boundary condition (2.10), we take an infinitesimal CV limit, which gives us

(2.17)$$\begin{eqnarray}\displaystyle u_{z}\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y & = & \displaystyle -[\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }(\unicode[STIX]{x1D749}(x+\unicode[STIX]{x0394}x,y)-\unicode[STIX]{x1D749}(x,y))\unicode[STIX]{x0394}y]\boldsymbol{\cdot }\boldsymbol{e}_{x}\nonumber\\ \displaystyle & & \displaystyle -\,[\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }(\unicode[STIX]{x1D749}(x,y+\unicode[STIX]{x0394}y)-\unicode[STIX]{x1D749}(x,y))\unicode[STIX]{x0394}x]\boldsymbol{\cdot }\boldsymbol{e}_{y},\end{eqnarray}$$

were we have $\unicode[STIX]{x0394}y=y_{2}-y_{1}$, $\unicode[STIX]{x0394}x=x_{2}-x_{1}$ and $x_{1}=x$ and $y_{1}=y$. Dividing both sides by $\unicode[STIX]{x0394}x\unicode[STIX]{x0394}y$ and using the definition of a derivative, we obtain

(2.18)$$\begin{eqnarray}u_{z}=-(\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{x}\unicode[STIX]{x1D749})\boldsymbol{\cdot }\boldsymbol{e}_{x}-(\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\unicode[STIX]{x2202}_{y}\unicode[STIX]{x1D749})\boldsymbol{\cdot }\boldsymbol{e}_{y}.\end{eqnarray}$$

This expression can be rewritten using double dot product as

(2.19)$$\begin{eqnarray}u_{z}=-\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{ : }\unicode[STIX]{x1D735}_{2}\unicode[STIX]{x1D749}.\end{eqnarray}$$

To obtain the transpiration length tensor, we express the tangential shear stress from (2.1) and insert the result into (2.19). Comparing the final result with (2.10) yields

(2.20)$$\begin{eqnarray}\unicode[STIX]{x1D648}=\unicode[STIX]{x1D707}\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\unicode[STIX]{x1D647}^{-1}.\end{eqnarray}$$

Recall that the tensor $\unicode[STIX]{x1D64D}_{\unicode[STIX]{x1D70F}}$ can be obtained as a post-processing step from the fundamental problems (FPsx), (FPsy) using the volume integral (2.14).

It is interesting to note that the velocity conditions (2.1), (2.10) can be written in a more compact form,

(2.21)$$\begin{eqnarray}\left(\begin{array}{@{}c@{}}u_{x}\\ u_{y}\\ u_{z}\end{array}\right)=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D613}_{xx} & 0 & 0\\ 0 & \unicode[STIX]{x1D613}_{xx} & 0\\ 0 & 0 & \unicode[STIX]{x1D614}_{xx}\end{array}\right)\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}\unicode[STIX]{x2202}_{z}u_{x}+\unicode[STIX]{x2202}_{x}u_{z}\\ \unicode[STIX]{x2202}_{z}u_{y}+\unicode[STIX]{x2202}_{y}u_{z}\\ \unicode[STIX]{x2202}_{z}u_{z}\end{array}\right),\end{eqnarray}$$

valid for an incompressible flow over isotropic geometries or for incompressible two-dimensional flows. The upper left $2\times 2$ block corresponds to the slip length tensor $\unicode[STIX]{x1D647}$ introduced before, while the lower right element $\unicode[STIX]{x1D614}_{xx}$ is the first term of the transpiration length tensor $\unicode[STIX]{x1D648}$. The equivalence with the previous formulation is obtained using the continuity equation, i.e. $\unicode[STIX]{x2202}_{z}u_{z}=-\unicode[STIX]{x2202}_{x}u_{x}-\unicode[STIX]{x2202}_{y}u_{y}$. The form (2.21) can be useful in practice, for example, if boundary conditions are imposed weakly in the finite element method. A similar set of boundary conditions – obtained by neglecting $\unicode[STIX]{x2202}_{x}u_{z}$ and $\unicode[STIX]{x2202}_{y}u_{z}$ – has been numerically investigated by Gómez de Segura et al. (Reference Gómez de Segura, Fairhall, MacDonald, Chung and García-Mayoral2018). In their work, the focus was on elucidating the turbulent flow response to the boundary condition (2.21) where all the coefficients for the slip and the transpiration lengths could take different values; however, connection to a surface texture geometry has never been made.

2.3 Numerical validation of velocity conditions

We consider a lid-driven cavity whose bottom surface is made of a texture with the characteristic length scale $\ell$ (figure 6a). The macroscopic length scale $H$ corresponds to the cavity width and depth. The scale separation parameter is equal to $\unicode[STIX]{x1D716}=\ell /H=0.1$. A no-slip condition is applied on all surfaces except the top wall, which moves with a prescribed velocity $(U_{0},0)$. Details about the numerical solver can be found in § C.1.

Figure 6.

Lid-driven cavity with a textured bottom. (a) Shows the computational domain used for the resolved simulations. The bottom surface consists of ten regular cavities. The domain for effective model simulation is shown in (b).

The moving upper wall generates a clockwise rotating vortex. This vortex imposes a negative shear on the rough surface. It also induces a downward mass flux at the right half of the cavity and an upward mass flux at the left half of the cavity. Near the textured surface, one can observe velocity fluctuations with a wavelength corresponding to the texture size $\ell$. To obtain macroscopic flow fields from geometry-resolving direct numerical simulation (DNS), we average out the microscale oscillations by creating an ensemble of $50$ DNS simulations. The ensemble consists of configurations in which the textured surface at the bottom of the cavity is incrementally shifted in the $x$ direction. The shift between neighbouring simulations is $\ell /50$. We believe that the ensemble average is the most appropriate way to obtain macroscopic variations due to the boundaries of the cavity. For periodic configurations, moving average or Fourier filtering could be equally accurate choices. The tangential and the transpiration velocities from the ensemble-averaged DNS along $z=0.2\ell$ are shown with black lines in figures 7(a) and 7(b), respectively.

Figure 7.

Tangential (a) and normal (b) velocities along the interface between the flow in the lid-driven cavity and the rough bottom. Dashed vertical lines show the streamwise locations where the DNS and model predictions are compared in table 1.

For comparison, we set up an effective simulation of the problem with the domain and boundary conditions shown in figure 6(b). We observed that in this particular shear-driven flow configuration the $\unicode[STIX]{x2202}_{x}u_{z}$ term in the shear stress at the wall is very small compared to $\unicode[STIX]{x2202}_{z}u_{x}$. Therefore, in the effective simulation we neglect $\unicode[STIX]{x2202}_{x}u_{z}$. We define the interface location at the previously selected coordinate $z_{i}=0.2\ell$. The coefficients for the boundary conditions (2.1), (2.10) – the slip and the transpiration lengths – are computed as described in §§ 2.1 and 2.2, using a FreeFEM++ open-source code (Lācis & Bagheri Reference Lācis and Bagheri2016 --2019). At this specific interface location, we have $\unicode[STIX]{x1D613}_{xx}=0.218\ell$ and $\unicode[STIX]{x1D614}_{xx}=0.110\ell$. The velocities at $z=0.2\ell$ from the effective simulation are compared to the ensemble-averaged DNS in figure 7(a,b). It is clear that the effective boundary conditions accurately predict the ensemble average (or the macroscopic variation) of both velocity components.

Table 1.

The slip length $\unicode[STIX]{x1D613}_{xx}$ and the transpiration length $\unicode[STIX]{x1D614}_{xx}$ for a range of interface locations $z_{i}$ above the textured surface. The effective tangential velocity $u_{x}$ is sampled at $(0.5H,z_{i})$. The effective transpiration velocity $u_{z}$ is sampled at $(0.25H,z_{i})$. The model predictions are finally normalised using the ensemble averaged results from DNS $\bar{u}_{x}$ and $\bar{u}_{z}$.

The results obtained using the TR model are reasonably accurate for different interface locations. To show this, we repeat the previous effective computation for a range of interface locations $z_{i}=(0.0,0.1,0.3,0.4,0.5)\ell$. For each interface coordinate, the effective coefficients $\unicode[STIX]{x1D613}_{xx}$ and $\unicode[STIX]{x1D614}_{xx}$ are recomputed using the fundamental problems (FPsx)–(FPsy), see table 1. To quantitatively present the TR model predictions, we select two streamwise positions at the interface plane, shown with vertical dashed lines in figure 7. In table 1, we present the model predictions of $u_{x}$ sampled at point $(0.5H,z_{i})$ for all interface locations. As the interface moves upwards – further away from the solid structures – the value of the predicted slip velocity increases due to a larger distance over which the viscous friction can bring the velocity to the no-slip value at the wall. This effect is correctly captured by the model through the linear increase of the slip length $\unicode[STIX]{x1D613}_{xx}$ (table 1). In other words, the information about the interface location is provided to the effective model through the adjustment of the coefficients. The interface normal velocity component $u_{z}$ and transpiration length $\unicode[STIX]{x1D614}_{xx}$ exhibit the same behaviour.

In the last two columns of table 1, we present the ratio between model predictions and the ensemble-averaged DNS results. We observe that, for all interface locations, the relative error ($1-u/\bar{u}$) is below $7\,\%$ for the slip velocity. The relative error for the transpiration velocity is below $14\,\%$. There is a trend of an increasing error as the interface is moved upwards, except for the interface location $z_{i}=0$. Despite the trend of increasing error with interface location, the TR model has a remarkably good accuracy taking into account that the transpiration velocity varies over two orders of magnitude (see 5th column of table 1). It has to be mentioned that the convergence of the ensemble averaging for the interface location $z_{i}=0$ is very slow and we used $250$ simulations to obtain a smooth result. The shift of the textured wall between each simulation is then $\ell /250$. We have carried out similar numerical computations on equilateral triangular surface texture and obtained the same behaviour as reported above.

This investigation shows that it is possible to adjust the interface height over distances $O(\ell )$ without a significant loss of accuracy. A similar conclusion about the interface location has already been reached numerically by Lācis & Bagheri (Reference Lācis and Bagheri2016) and theoretically by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) for the slip velocity alone. However, as a ‘rule-of-thumb’, we suggest placing the interface as close to the solid structure as possible without intersecting the solids.

3.1 Velocity boundary conditions

For a porous surface, the tangential velocity boundary condition is identical to the textured surface, i.e. the slip condition (2.1). The interface normal velocity condition, on the other hand, becomes

(3.1)$$\begin{eqnarray}u_{n}=u_{n}^{-}-\unicode[STIX]{x1D648}\boldsymbol{ : }\unicode[STIX]{x1D735}_{2}\boldsymbol{u}_{t}\quad \text{on }z=z_{i},\end{eqnarray}$$

where $u_{n}^{-}$ is the interface normal Darcy velocity, satisfying

(3.2)$$\begin{eqnarray}u_{n}^{-}=\left(-\frac{\unicode[STIX]{x1D646}}{\unicode[STIX]{x1D707}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p^{-}\right)\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

This term is induced by mass conservation between the free fluid and the porous domain.

The tensors $\unicode[STIX]{x1D648}$ and $\unicode[STIX]{x1D647}$ are determined by solving the fundamental problems (FPsx)–(FPsy). The only difference from the textured surface is that the solid structures within the interface cell represent the porous material. Consequently, all the elements in $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D648}$ can be obtained through expressions (2.9) and (2.20), respectively. The interior permeability tensor ($\unicode[STIX]{x1D646}$) of the porous medium is computed through a set of Stokes equations in a bulk unit cell (Whitaker Reference Whitaker1998; Mei & Vernescu Reference Mei and Vernescu2010).

3.2 Pressure boundary condition

The pressure boundary condition for a general 3-D porous surface, obeying assumptions (A1)–(A3), is obtained through a balance between the normal free-fluid stress and the stress from the porous material, i.e.

(3.3)$$\begin{eqnarray}-p+2\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{z}u_{z}=-p^{-}+\boldsymbol{f}^{(1)}\boldsymbol{\cdot }\boldsymbol{u}^{-}+\boldsymbol{f}^{(2)}\boldsymbol{\cdot }\boldsymbol{u}_{t}\quad \text{on }z=z_{i}.\end{eqnarray}$$

The normal stress from the porous material consists of the pore pressure $p^{-}$ and two friction coefficients, $\boldsymbol{f}^{(1)}$ and $\boldsymbol{f}^{(2)}$. The coefficient $\boldsymbol{f}^{(1)}$ describes the interface normal resistance that the Darcy flow $\boldsymbol{u}^{-}$ must overcome to transport mass and momentum across and along the interface. The coefficient $\boldsymbol{f}^{(2)}$ provides the interface normal force due to the slip velocity near the interface. It exists only for anisotropic and misaligned surface geometries, similarly to the stabilisation parameter (see (1.2)) derived by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) and Carraro et al. (Reference Carraro, Marušić-Paloka and Mikelić2018). Misalignment in this context is a difference between the principal axes of the pore geometry and the coordinate axes.

The friction coefficients are again determined by considering the interface cell (figure 8b). The difference from the textured surface is the existence of the pore pressure resulting in an additional forcing term in the governing equations of the interface cell, yielding

(3.4)$$\begin{eqnarray}\displaystyle \displaystyle -\unicode[STIX]{x1D735}\hat{p}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}} & = & \displaystyle -\unicode[STIX]{x1D6FF}(\hat{z}-\hat{z}_{i})\unicode[STIX]{x1D749}+H(\hat{z}_{i}-\hat{z})\unicode[STIX]{x1D735}p^{-},\end{eqnarray}$$

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

Darcy’s law is valid only in the porous material, therefore the pressure gradient forcing is considered only below the interface. A 1-D Heaviside step function $H(\hat{z})$ is used to distinguish between regions above and below the interface. Boundary conditions for the interface cell are the same as for the equations (2.2)–(2.3) except at the bottom of the domain, where we impose the interior solution corresponding to the Darcy flow due to the same pressure gradient $\unicode[STIX]{x1D735}p^{-}$ (Whitaker Reference Whitaker1998; Mei & Vernescu Reference Mei and Vernescu2010; Lācis & Bagheri Reference Lācis and Bagheri2016).

Then, we use the linearity assumption (A1) and write the pressure as

(3.6)$$\begin{eqnarray}\hat{p}=\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\unicode[STIX]{x1D749}-\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p^{-}.\end{eqnarray}$$

Here, $\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}$ and $\hat{\boldsymbol{r}}_{p}$ are the response operators related to the shear stress $\unicode[STIX]{x1D749}$ and the pressure gradient $\unicode[STIX]{x1D735}p^{-}$, respectively. This expression is expanded as

(3.7)$$\begin{eqnarray}\displaystyle \hat{p} & = & \displaystyle \hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }(\unicode[STIX]{x1D70F}_{x}\boldsymbol{e}_{x}+\unicode[STIX]{x1D70F}_{y}\boldsymbol{e}_{y})-\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }(\unicode[STIX]{x2202}_{x}p^{-}\boldsymbol{e}_{x}+\unicode[STIX]{x2202}_{y}p^{-}\boldsymbol{e}_{y}+\unicode[STIX]{x2202}_{z}p^{-}\boldsymbol{e}_{z})\nonumber\\ \displaystyle & = & \displaystyle (\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{x})\unicode[STIX]{x1D70F}_{x}+(\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{y})\unicode[STIX]{x1D70F}_{y}-(\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{x})\unicode[STIX]{x2202}_{x}p-(\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{y})\unicode[STIX]{x2202}_{y}p-(\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{z})\unicode[STIX]{x2202}_{z}p\nonumber\\ \displaystyle & = & \displaystyle \hat{p}^{(\unicode[STIX]{x1D70F}x)}\frac{\unicode[STIX]{x1D70F}_{x}}{\unicode[STIX]{x1D707}}+\hat{p}^{(\unicode[STIX]{x1D70F}y)}\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D707}}-\hat{p}^{(px)}\frac{\unicode[STIX]{x2202}_{x}p}{\unicode[STIX]{x1D707}}-\hat{p}^{(py)}\frac{\unicode[STIX]{x2202}_{y}p}{\unicode[STIX]{x1D707}}-\hat{p}^{(pz)}\frac{\unicode[STIX]{x2202}_{z}p}{\unicode[STIX]{x1D707}},\end{eqnarray}$$

where we have defined

(3.8a,b)$$\begin{eqnarray}\hat{p}^{(\unicode[STIX]{x1D70F}x)}=\unicode[STIX]{x1D707}\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{x},\quad \hat{p}^{(\unicode[STIX]{x1D70F}y)}=\unicode[STIX]{x1D707}\hat{\boldsymbol{r}}_{\unicode[STIX]{x1D70F}}\boldsymbol{\cdot }\boldsymbol{e}_{y}\end{eqnarray}$$

and

(3.9a-c)$$\begin{eqnarray}\hat{p}^{(px)}=\unicode[STIX]{x1D707}\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{x},\quad \hat{p}^{(py)}=\unicode[STIX]{x1D707}\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{y},\quad \hat{p}^{(pz)}=\unicode[STIX]{x1D707}\hat{\boldsymbol{r}}_{p}\boldsymbol{\cdot }\boldsymbol{e}_{z}.\end{eqnarray}$$

Note that $\hat{p}^{(\unicode[STIX]{x1D70F}x)}$ and $\hat{p}^{(\unicode[STIX]{x1D70F}y)}$ are the pressure fields appearing in the fundamental problems (FPsx), (FPsy); they are the pressure responses to the interface shear forcing in the $x$ and $y$ directions, respectively. Furthermore, $\hat{p}^{(px)}$, $\hat{p}^{(py)}$ and $\hat{p}^{(pz)}$ are the pressure fields associated with the following three fundamental problems

(F P p x)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\hat{p}^{(px)}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}}^{(px)}=-\unicode[STIX]{x1D707}H(\hat{z}_{i}-\hat{z})\boldsymbol{e}_{x},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{(px)}=0;\end{eqnarray}$$

(F P p y)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\hat{p}^{(py)}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}}^{(py)}=-\unicode[STIX]{x1D707}H(\hat{z}_{i}-\hat{z})\boldsymbol{e}_{y},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{(py)}=0;\end{eqnarray}$$

(F P p z)$$\begin{eqnarray}-\unicode[STIX]{x1D735}\hat{p}^{(pz)}+\unicode[STIX]{x1D707}\unicode[STIX]{x0394}\hat{\boldsymbol{u}}^{(pz)}=-\unicode[STIX]{x1D707}H(\hat{z}_{i}-\hat{z})\boldsymbol{e}_{z},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\hat{\boldsymbol{u}}^{(pz)}=0.\end{eqnarray}$$

These problems describe the response to the pressure gradient forcing along the three coordinates and have been previously derived by Lācis & Bagheri (Reference Lācis and Bagheri2016) using formal multi-scale expansion. Keep in mind that in (FPpx)–(FPpz) the fields $\hat{\boldsymbol{u}}^{(px)}$, $\hat{\boldsymbol{u}}^{(py)}$ and $\hat{\boldsymbol{u}}^{(pz)}$ have units of m$^{2}$, similar to the permeability of a porous medium, see table 2.

Table 2.

Summary of pressure and velocity fields with the corresponding units.

From the five fundamental problems (FPsx)–(FPpz), we can determine the resistance vectors $\boldsymbol{f}^{(1)}$ and $\boldsymbol{f}^{(2)}$ for the pressure condition (3.3). We start with $\boldsymbol{f}^{(2)}$. It generates a pressure jump

(3.10)$$\begin{eqnarray}p^{-}-p=\boldsymbol{f}^{(2)}\boldsymbol{\cdot }\boldsymbol{u}_{t}=\tilde{\boldsymbol{f}}^{(2)}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D749}}{\unicode[STIX]{x1D707}},\end{eqnarray}$$

where $\tilde{\boldsymbol{f}}^{(2)}=\boldsymbol{f}^{(2)}\boldsymbol{\cdot }\unicode[STIX]{x1D647}$. The pressure field response (3.7) due to the shear is

(3.11)$$\begin{eqnarray}\hat{p}=\hat{p}^{(\unicode[STIX]{x1D70F}x)}\frac{\unicode[STIX]{x1D70F}_{x}}{\unicode[STIX]{x1D707}}+\hat{p}^{(\unicode[STIX]{x1D70F}y)}\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

Then, we need to relate the effective pressures in porous and free-fluid regions to the linear pressure responses $\hat{p}^{(\unicode[STIX]{x1D70F}x)}$ and $\hat{p}^{(\unicode[STIX]{x1D70F}y)}$ in the interface cell. For the velocity, a simple plane average at the interface (2.7) is sufficient. However, for the pressure condition a single pressure value will not provide the necessary information about the pressure jump. Therefore, we define the effective pressure in the interior and free fluid as

(3.12a,b)$$\begin{eqnarray}p^{-}=\frac{1}{V_{f}}\int _{\ell ^{2}}\int _{\hat{z}_{b}}^{\hat{z}_{b}+\ell }\hat{p}\,\text{d}V=\langle \hat{p}\rangle ^{-},\quad p=\frac{1}{V_{f}}\int _{\ell ^{2}}\int _{\hat{z}_{t}-\ell }^{\hat{z}_{t}}\hat{p}\,\text{d}V=\langle \hat{p}\rangle ^{+},\end{eqnarray}$$

with $V_{f}$ corresponding to the fluid volume in the integration region. To neglect any transition effects of the pressure field near the interface, these volume averages are taken at the bottom and at the top of the interface cell. In this way, the averaging operation is sufficiently far away from the interface to obtain a representative pressure value for the interior and the free fluid.

Now we insert the pressure field decomposition (3.11) into (3.12) and we take the difference between the interior pressure and the free-fluid pressure,

(3.13)$$\begin{eqnarray}p^{-}-p=\left(\langle \hat{p}^{(\unicode[STIX]{x1D70F}x)}\rangle ^{-}-\langle \hat{p}^{(\unicode[STIX]{x1D70F}x)}\rangle ^{+}\right)\frac{\unicode[STIX]{x1D70F}_{x}}{\unicode[STIX]{x1D707}}+\left(\langle \hat{p}^{(\unicode[STIX]{x1D70F}y)}\rangle ^{-}-\langle \hat{p}^{(\unicode[STIX]{x1D70F}y)}\rangle ^{+}\right)\frac{\unicode[STIX]{x1D70F}_{y}}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

By comparing the above to (3.10), we obtain

(3.14a,b)$$\begin{eqnarray}\tilde{f}_{x}^{(2)}=\langle \hat{p}^{(\unicode[STIX]{x1D70F}x)}\rangle ^{-}-\langle \hat{p}^{(\unicode[STIX]{x1D70F}x)}\rangle ^{+},\quad \tilde{f}_{y}^{(2)}=\langle \hat{p}^{(\unicode[STIX]{x1D70F}y)}\rangle ^{-}-\langle \hat{p}^{(\unicode[STIX]{x1D70F}y)}\rangle ^{+}.\end{eqnarray}$$

We emphasise that $\hat{p}^{(\unicode[STIX]{x1D70F}x)}$ and $\hat{p}^{(\unicode[STIX]{x1D70F}y)}$ are pressure fields in the interface cell generated due to shear stress forcing (figure 8c,d) and can be computed from the fundamental problems (FPsx), (FPsy). Finally, the resistance vector $\boldsymbol{f}^{(2)}$, appearing in the front of the slip velocity in (3.3), is obtained from

(3.15)$$\begin{eqnarray}\boldsymbol{f}^{(2)}=\tilde{\boldsymbol{f}}^{(2)}\boldsymbol{\cdot }\unicode[STIX]{x1D647}^{-1}.\end{eqnarray}$$

The procedure to get this friction coefficient is similar to the one reported by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) and Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013).

Figure 9.

A lid-driven cavity with a porous bed. (a) Shows the dimensions of the computational domain. The dashed vertical line indicates the streamwise position where the DNS and the effective models are compared. (b) Depicts an enlarged view of the porous materials showing the microscale geometry of the three test cases considered. The interface is located at a distance of $0.1\ell$ above the solid structure. (c) Shows the domain for the continuum description.

We turn our attention to the resistance coefficient $\boldsymbol{f}^{(1)}$. The pressure jump condition (3.3) due to the Darcy velocity is

(3.16)$$\begin{eqnarray}p^{-}-p=\boldsymbol{f}^{(1)}\boldsymbol{\cdot }\boldsymbol{u}^{-}=-\tilde{\boldsymbol{f}}^{(1)}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D735}p^{-}}{\unicode[STIX]{x1D707}},\end{eqnarray}$$

where $\tilde{\boldsymbol{f}}^{(1)}=-\boldsymbol{f}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D646}$. The pressure field response (3.7), corresponding to the pore pressure gradient forcing, is

(3.17)$$\begin{eqnarray}\hat{p}=-\hat{p}^{(px)}\frac{\unicode[STIX]{x2202}_{x}p^{-}}{\unicode[STIX]{x1D707}}-\hat{p}^{(py)}\frac{\unicode[STIX]{x2202}_{y}p^{-}}{\unicode[STIX]{x1D707}}-\hat{p}^{(pz)}\frac{\unicode[STIX]{x2202}_{z}p^{-}}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

Using (3.12), we can express the pressure jump as

(3.18)$$\begin{eqnarray}\displaystyle p-p^{-} & = & \displaystyle \langle \hat{p}\rangle ^{+}-\langle \hat{p}\rangle ^{-}=\left(\langle \hat{p}^{(px)}\rangle ^{-}-\langle \hat{p}^{(px)}\rangle ^{+}\right)\frac{\unicode[STIX]{x2202}_{x}p^{-}}{\unicode[STIX]{x1D707}}\nonumber\\ \displaystyle & & \displaystyle +\,\left(\langle \hat{p}^{(py)}\rangle ^{-}-\langle \hat{p}^{(py)}\rangle ^{+}\right)\frac{\unicode[STIX]{x2202}_{y}p^{-}}{\unicode[STIX]{x1D707}}+\left(\langle \hat{p}^{(pz)}\rangle ^{-}-\langle \hat{p}^{(pz)}\rangle ^{+}\right)\frac{\unicode[STIX]{x2202}_{z}p^{-}}{\unicode[STIX]{x1D707}}.\end{eqnarray}$$

The comparison of (3.16) and (3.18) yields

(3.19a-c)$$\begin{eqnarray}\tilde{f}_{x}^{(1)}=\langle \hat{p}^{(px)}\rangle ^{-}-\langle \hat{p}^{(px)}\rangle ^{+},\quad \tilde{f}_{y}^{(1)}=\langle \hat{p}^{(py)}\rangle ^{-}-\langle \hat{p}^{(py)}\rangle ^{+},\quad \tilde{f}_{z}^{(1)}=\langle \hat{p}^{(pz)}\rangle ^{-}-\langle \hat{p}^{(pz)}\rangle ^{+}.\end{eqnarray}$$

We recall that the pressure fields in the interface cell are generated by the pore pressure gradient forcing below the interface (figure 8e–g), and they are computed from fundamental problems (FPpx)–(FPpz). The final form of the friction coefficient is,

(3.20)$$\begin{eqnarray}\boldsymbol{f}^{(1)}=-\tilde{\boldsymbol{f}}^{(1)}\boldsymbol{\cdot }\unicode[STIX]{x1D646}^{-1}.\end{eqnarray}$$

This friction coefficient term, which to the best of authors’ knowledge is reported for the first time, is particularly important for capturing the correct pressure jump across the interface for layered problems, as we demonstrate in the next section.

3.3 Validation of the TR model for porous surfaces

We consider the same flow configuration as in § 2.3, but we replace the bottom textured wall with a porous surface, see figure 9(a). The porous medium consists of a periodic distribution of solid inclusions with a characteristic length scale $\ell$ (see $\ell \times \ell$ square in figure 9a). The width and the height of the cavity are $H$, while the depth of the porous material is $H/2$. The scale separation parameter is again $\unicode[STIX]{x1D716}=\ell /H=0.1$. The flow reaches the interior seepage velocity quickly (Lācis & Bagheri Reference Lācis and Bagheri2016; Lācis, Zampogna & Bagheri Reference Lācis, Zampogna and Bagheri2017); therefore, a porous material containing only five repeating structures in depth is sufficient for using Darcy’s equation in the interior.

To demonstrate the generality of the TR model, we consider three kinds of porous geometries, see figure 9(b). Configuration (i) has circular solid inclusions, which results in an isotropic porous medium. The anisotropic elliptic inclusions considered in configuration (ii) are the same as investigated by Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013). The last geometry (iii) has isotropic circular inclusions with the interface layer different from the interior. The porosity and geometrical details of the three configurations are listed in table 3.

We carry out fully resolved DNS. For each configuration, the mean over an ensemble of $50$ shifted porous beds is computed. The free flow is similar to the flow in the lid-driven cavity with the textured bottom (§ 2.3) except that there are mass fluxes in and out of the porous material. The averaged DNS will be compared to effective representations of the porous bed (figure 9c). Since the flow configuration is similar to the one reported in § 2.3, we will again neglect $\unicode[STIX]{x2202}_{x}u_{z}$ in the effective simulations. We have checked that $\unicode[STIX]{x2202}_{x}u_{z}$ is much smaller than $\unicode[STIX]{x2202}_{z}u_{x}$. Within the porous domain, we employ Darcy’s law, where the only unknown quantity is the pore pressure $p^{-}$. The interior permeability tensors ($\unicode[STIX]{x1D646}$) for all geometries are listed in the last column of table 3. They were obtained by solving a set of Stokes equations in a periodic unit cell in the bulk (Whitaker Reference Whitaker1998; Mei & Vernescu Reference Mei and Vernescu2010). A Neumann condition on pore pressure $\unicode[STIX]{x1D735}p^{-}\boldsymbol{\cdot }\boldsymbol{n}=0$ is enforced at solid boundaries of the cavity. This condition corresponds to zero fluid flux through the wall. Boundary conditions for the free fluid remain the same as in § 2.3. We place the interface at a distance $z_{i}=0.1\ell$ above the solid structures; the interface location $z_{i}=0.0$ was not accessible due to meshing issues. We compute the effective parameters appearing in boundary conditions (2.1), (3.1), (3.3) using the procedure explained in §§ 2.1, 2.2 and 3.2. The slip and the transpiration lengths, reported in table 4, have nearly the same value for the three configurations because the porous materials have similar porosity near the interface (table 3).

Table 3.

Geometrical properties of the porous media considered in this work (see graphical representation in figure 9b). The subscript $i$ corresponds to the interface and the subscript $b$ corresponds to the bulk. The porosity $\unicode[STIX]{x1D719}$ is defined as the ratio between the solid volume and the fluid volume. The last column shows the interior permeability tensor $\unicode[STIX]{x1D646}$.

Figure 10.

The pressure (a) and the transpiration velocity (b) profiles of the lid-driven cavity with layered isotropic porous bed (iii). (a) Shows the distribution of pressure along the vertical dashed line in figure 9(a). The grey shaded region corresponds to the porous material.

Table 4.

The slip length, the transpiration length and the resistance coefficients $\boldsymbol{f}^{(1)}$ and $\boldsymbol{f}^{(2)}$ for the porous medium geometries shown graphically in figure 9(b). The last three columns show the ratio between the model and the DNS results for the slip and the transpiration velocities.

To validate the velocity conditions (2.1), (3.1), we sample the ensemble-averaged DNS and the effective model at coordinates $(0.5H,z_{i})$ and $(0.25H,z_{i})$ for slip and transpiration velocities, respectively. The ratios between the model predictions and the DNS are given in last columns of table 4. It is clear that the slip velocity is predicted as accurately as for the textured surfaces. In the second to last column of table 4 we list the ratio between the Darcy transpiration velocity $u_{z}^{-}$ – sampled just below point $(0.25H,z_{i})$ – and the DNS result. The Darcy velocity alone is a rather inaccurate predictor and has a relative error of up to $37\,\%$. The agreement between the TR model (3.1) – that augments the Darcy contribution with the term containing the transpiration length – and the DNS is better as the relative error is smaller than $15\,\%$. If a smaller error is desired, one can use the full multi-scale expansion model to reduce the error to below $8\,\%$ (Sudhakar et al. Reference Sudhakar, Lācis, Pasche and Bagheri2019).

To validate the pressure condition (3.3), we analyse the ensemble-averaged DNS of configuration (iii). The pressure $\hat{p}$ along the vertical dashed lined in figure 9(a) is shown in figure 10(a) with a solid black line. The region corresponding to the porous domain is shadowed. We observe that the pressure field undergoes a sharp variation when transitioning from the free-fluid region to the porous medium, as shown by an inset in figure 10(a). The sharp variation indicates that there could be a pressure jump in the effective representation. Note that there are microscale oscillations in the pressure field, because the employed ensemble average filters out only the microscale variations in the $x$ direction. For comparison, the pressure obtained from the effective model is shown using dotted blue symbols in figure 10(a). We can observe that the agreement between model predictions and DNS is good both in free fluid and interior, while the sharp variation in the near vicinity of the interface is not modelled. This is a direct consequence of having an infinitely thin interface in the effective representation, which condenses all variations near the boundary to a single line. The same quantitative agreement we observed for the configurations (i) and (ii).

To show the importance of the resistance coefficients $\boldsymbol{f}^{(1)}$ and $\boldsymbol{f}^{(2)}$, we carried out two more effective simulations. In the first one – called ‘leading order’ – we set $\boldsymbol{f}^{(1)}=0$ in the pressure condition (3.3). This corresponds to pressure condition proposed by Marciniak-Czochra & Mikelić (Reference Marciniak-Czochra and Mikelić2012) Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013) and Carraro et al. (Reference Carraro, Marušić-Paloka and Mikelić2018). In the second one – called ‘pressure continuity’ – we impose $p=p^{-}$, which has been a common approach in the past (Ene & Sanchez-Palencia Reference Ene and Sanchez-Palencia1975; Levy & Sanchez-Palencia Reference Levy and Sanchez-Palencia1975; Hou et al. Reference Hou, Holmes, Lai and Mow1989; Lācis & Bagheri Reference Lācis and Bagheri2016). Results from the leading-order model and the pressure continuity model are reported in figure 10(a) using crosses and a dashed curve, respectively. We observe that both conditions result in a poor agreement between the model and the DNS if compared to the TR model.

In addition, an inaccurate pressure condition has consequences for the flow field, as illustrated in figure 10(b). There, the transpiration velocity at the interface is shown using the same symbols as in figure 10(a). We observe that the error in the pressure condition can lead to significantly different – and inaccurate – vertical velocity predictions. The coefficient $\boldsymbol{f}^{(1)}$ imposes a larger resistance for the wall-normal velocity, and thus it decreases the transpiration by precisely the correct amount.

5.1 One-to-one comparison between the TR model and the MSE

An expansion up to the order $O(\unicode[STIX]{x1D716})$ results in the following boundary conditions (A 32) for the free-fluid velocity at the interface,

(5.1)$$\begin{eqnarray}\boldsymbol{u}=\underbrace{\unicode[STIX]{x1D647}_{e}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D749}}{\unicode[STIX]{x1D707}}}_{O(U^{s})}-\underbrace{\unicode[STIX]{x1D646}_{e}\boldsymbol{\cdot }\frac{\unicode[STIX]{x1D735}p^{-}}{\unicode[STIX]{x1D707}}+\unicode[STIX]{x1D648}_{e}\boldsymbol{ : }\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D749}}{\unicode[STIX]{x1D707}}}_{O(\unicode[STIX]{x1D716}U^{s})}+O(\unicode[STIX]{x1D716}^{2}U^{s}),\end{eqnarray}$$

where $U^{s}$ is a characteristic magnitude of the slip velocity. The tensor $\unicode[STIX]{x1D648}_{e}$ is a third-rank tensor with $81$ elements, while the tensors $\unicode[STIX]{x1D647}_{e}$ and $\unicode[STIX]{x1D646}_{e}$ are both second-rank tensors with $9$ elements. The double dot operation between a third-rank tensor $\unicode[STIX]{x1D63C}$ and a second-rank tensor $\unicode[STIX]{x1D63D}$ is defined as $\unicode[STIX]{x1D63C}\boldsymbol{ : }\unicode[STIX]{x1D63D}=\unicode[STIX]{x1D608}_{ijk}\unicode[STIX]{x1D609}_{jk}$, where summation over repeating indices is implied.

The leading-order term $O(U^{s})$ of the velocity condition (5.1) is the slip term. The slip tensor $\unicode[STIX]{x1D647}$ in the TR model (2.1) corresponds to the upper left $2\times 2$ block of $\unicode[STIX]{x1D647}_{e}$. From mass conservation arguments, it can be shown that the last row and column of tensor $\unicode[STIX]{x1D647}_{e}$ appearing in (5.1) are zero. Thus there is no transpiration velocity at $O(U^{s})$.

There are two higher-order $O(\unicode[STIX]{x1D716}U^{s})$ terms in the velocity condition (5.1); a Darcian term related to the pore pressure gradient and a term related to the variation of the shear stress. In the TR model, the Darcy contribution to the tangential velocity components at the interface is neglected. In other words, the TR model has the Darcy contribution only for the wall-normal transpiration component (3.1). One may again show from mass conservation, that the last row of $\unicode[STIX]{x1D646}_{e}$ in equation (5.1) is equal to the last row of the interior permeability tensor $\unicode[STIX]{x1D646}$. Therefore, the Darcy term $u_{z}^{-}$ in the TR model (3.2) corresponds to $(\unicode[STIX]{x1D646}_{e}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p^{-}/\unicode[STIX]{x1D707})\cdot \hat{n}$. Finally, the term related to the variation of the shear stress in (5.1) is compared with the transpiration boundary conditions in the TR model (2.10). We observe that the TR model contains some of the next-order terms arising from the variation of the shear stress; the TR model contains only a second-rank tensor corresponding to the tangential shear stress variations in the tangential directions, while the full third-rank tensor in the MSE corresponds to all shear stress variations in all directions. The contribution from the variation of the shear stress to the tangential velocity components is neglected in the TR model.

The boundary condition for the pressure derived using MSE (A 33) is

(5.2)$$\begin{eqnarray}p^{-}-p=\underbrace{\boldsymbol{b}\boldsymbol{\cdot }\unicode[STIX]{x1D749}}_{O(\unicode[STIX]{x0394}P)}-\underbrace{\boldsymbol{a}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p^{-}+\unicode[STIX]{x1D63E}\boldsymbol{ : }\unicode[STIX]{x1D735}\unicode[STIX]{x1D749}}_{O(\unicode[STIX]{x1D716}\unicode[STIX]{x0394}P)}+O(\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x0394}P),\end{eqnarray}$$

where $\unicode[STIX]{x0394}P$ is the characteristic magnitude of a pressure drop in the system. Here, $\boldsymbol{b}$ and $\boldsymbol{a}$ are vectors and $\unicode[STIX]{x1D63E}$ is a second-rank tensor. The leading term $O(\unicode[STIX]{x0394}P)$ of expression (5.2) induces a pressure jump proportional to the shear stress. It can be shown through mass conservation and force balance that the last element in $\boldsymbol{b}$ (corresponding to shear stress $2\unicode[STIX]{x1D707}\unicode[STIX]{x2202}_{z}u_{z}$) is always equal to $-1$. Therefore this term can be transferred to the left-hand side and grouped together with ‘$-p$’ to yield the total free-fluid stress, as appearing in the TR model (3.3). Furthermore, we assert that the first two elements of the vector $\boldsymbol{b}$ in expression (5.2) correspond to the friction factor $\boldsymbol{f}^{(2)}$ in (3.3), which can be confirmed by replacing the shear stress in expression (5.2) with $\unicode[STIX]{x1D707}\unicode[STIX]{x1D647}^{-1}\boldsymbol{u}_{t}$. Thus there is a full overlap of the leading-order pressure jump terms between the TR and the MSE. Regarding the last two higher-order $O(\unicode[STIX]{x1D716}\unicode[STIX]{x0394}P)$ terms in the pressure condition (5.2): the vector $\boldsymbol{a}$ corresponds to the friction factor $\boldsymbol{f}^{(1)}$, which is confirmed by replacing the pore pressure in equation (5.2) with $\unicode[STIX]{x1D707}\unicode[STIX]{x1D646}^{-1}\boldsymbol{u}^{-}$, while the term corresponding to the variations of shear stress is completely neglected for the pressure condition in the TR model.

Table 6.

Summary of the TR model boundary condition terms at the orders at which corresponding terms emerge from the formal multi-scale expansion.

In table 6 we group the different terms of the TR model according to the order at which they emerge in the multi-scale expansion. The slip velocity in the TR model contains only the leading-order term with a reduced shear stress vector, while the transpiration condition and the pressure condition contains all the leading-order contributions and some of the next-order corrections.

5.2 Accuracy of the TR model compared to the MSE

We compare the flow and the pressure in the lid-driven cavity with a textured and a porous surface computed from fully resolved ensemble-averaged DNS with three different effective models; (i) the zeroth-order model, containing only leading-order terms from the MSE condition (5.1)–(5.2), (ii) the first-order model, containing all terms from the MSE condition (5.1)–(5.2) and (iii) the TR model.

Figure 12.

The transpiration velocity along the interface between the flow in the lid-driven cavity and the rough surface (a). Pressure distribution along the vertical slice for the cavity flow over the layered porous surface (b). Interface location $z_{i}=0.1\ell$.

The transpiration velocity along the interface – at the same textured wall as discussed in § 2.3 – is shown in figure 12(a). The pressure distribution along the vertical slice – of the same layered geometry as discussed in § 3.3 – is shown in figure 12(b). The interface in both cases is located at $z_{i}=0.1\ell$. As expected, the TR model has a clear improvement over the zeroth-order MSE model. It is observed that the first-order MSE model and the TR model nearly overlap for both the transpiration and the pressure fields. We can thus conclude that the TR model provides nearly as good approximation as the MSE model, but with significantly reduced complexity.

As a final remark, we note that the TR model is not mathematically (or asymptotically) fully consistent. If mathematical rigour is sought, all the next-order terms for the slip velocity, transpiration velocity and pressure condition should be taken into account. However, as we have demonstrated for the turbulent channel flow, the transpiration velocity from the TR model induces a large shift of the mean velocity profile changing the sign of the roughness function. Thus this example shows that higher-order terms can be physically important.