2.1. Governing equations and domain decomposition

Let us consider the turbulent flow of a viscous, incompressible, Newtonian fluid in a channel delimited from one side (at $\hat y = 2H$ ) by a smooth, impermeable wall and from the other side (at $\hat y \leq 0$ ) by a permeable substrate constructed with spanwise-elongated ( $\hat z$ -aligned) solid inclusions, regularly arranged with given periodicity $\ell$ in the streamwise and wall-normal directions ( $\hat x$ and $\hat y$ , respectively), see figure 1. The velocity components ( $\hat {u}_1=\hat {u}$ , $\hat {u}_2=\hat {v}$ , $\hat {u}_3=\hat {w}$ ) and the pressure $\hat {p}$ are the dependent variables, to be evaluated over space ( $\hat {x}_1= \hat {x}$ , $\hat {x}_2= \hat {y}$ , $\hat {x}_3= \hat {z}$ ) and time $\hat t$ . The conservation equations governing the flow can be expressed as follows:

(2.1) \begin{align} \frac {\partial \hat {u}_{i}}{\partial \hat {x}_{i}}=0, \quad \rho \left (\frac {\partial \hat {u}_{i}}{\partial \hat {t}}+\hat {u}_{j} \frac {\partial \hat {u}_{i}}{\partial \hat {x}_{j}}\right )=-\frac {\partial \hat {p}}{\partial \hat {x}_{i}}+\mu \frac {\partial ^{2} \hat {u}_{i}}{\partial \hat {x}_{j}^{2}}, \end{align}

with $\rho$ and $\mu$ the fluid density and dynamic viscosity, respectively.

We identify two characteristic length scales: a microscopic one, $\tilde \ell$ , characterizing the porous bed, and a macroscopic one, $\mathcal {L}$ , related to the large-scale motion in the channel. Provided that the two length scales are well-separated, i.e. $\tilde \ell \ll \mathcal {L}$ , it is possible to manipulate the microscale problem by means of an asymptotic analysis in terms of a small parameter $\epsilon =\tilde \ell /\mathcal {L} \ll 1$ . As illustrated in figure 1, the flow domain is decomposed into three distinct subdomains: a channel-flow region away from the porous–free-fluid interface (superscript ‘ $\mathcal {C}$ ’), an interfacial region (superscript ‘ $\mathcal {I}$ ’) and a region within the porous layer away from boundaries, governed by Darcy’s law (superscript ‘ $\mathcal {P}$ ’). Correspondingly, the following three sets of normalized variables are proposed:

(2.2a) \begin{align} X_{i} & = \frac {\hat {x}_i}{\mathcal {L}}, \quad U_{i}^{\mathcal {C}}=\frac {\hat {u}}{\mathcal {U}}, \quad P^{\mathcal {C}}=\frac {\hat {p}}{\rho \mathcal {U}^{2}}, \end{align}

(2.2b) \begin{align} x_{i} & = \frac {\hat {x}_i}{\tilde \ell }, \quad U_{i}^{\mathcal {I}}=\frac {\hat {u}}{\epsilon \mathcal {U}}, \quad P^{\mathcal {I}}=\frac {\hat {p}}{\mu \mathcal {U}/\mathcal {L}}, \end{align}

(2.2c) \begin{align} x_{i} & = \frac {\hat {x}_i}{\tilde \ell }, \quad U_{i}^{\mathcal {P}}=\frac {\hat {u}}{\epsilon ^2 \mathcal {U}}, \quad P^{\mathcal {P}}=\frac {\hat {p}}{\mu \mathcal {U}/\mathcal {L}}, \\[6pt] \nonumber \end{align}

where $\mathcal {U}$ is a suitable macroscopic velocity scale; a discussion on the proper selection of scales is provided later. Based on the normalization above, the governing (2.1) can be recast into the following dimensionless forms in the $\bullet ^{\mathcal {C}}$ , $\bullet ^{\mathcal {I}}$ and $\bullet ^{\mathcal {P}}$ regions, respectively:

(2.3a) \begin{align} \frac {\partial U_{i}^{\mathcal {C}}}{\partial X_{i}}=0, \quad \frac {\partial U_{i}^{\mathcal {C}}}{\partial t}+U_{j}^{\mathcal {C}} \frac {\partial U_{i}^{\mathcal {C}}}{\partial X_{j}}=-\frac {\partial P^{\mathcal {C}}}{\partial X_{i}}+\frac {1}{Re} \frac {\partial ^{2} U_{i}^{\mathcal {C}}}{\partial X_{j}^{2}} , \end{align}

(2.3b) \begin{align} \frac {\partial {U}^{\mathcal {I}}_{i}}{\partial x_{i}}=0, \quad \epsilon ^{2} {Re}\left (\frac {\partial {U}_{i}^{\mathcal {I}}}{\partial t}+{U}_{j}^{\mathcal {I}} \frac {\partial {U}_{i}^{\mathcal {I}}}{\partial x_{j}}\right )=-\frac {\partial P^{\mathcal {I}}}{\partial x_{i}}+\frac {\partial ^{2} {U}_{i}^{\mathcal {I}}}{\partial x_{j}^{2}} , \end{align}

(2.3c) \begin{align} \epsilon \frac {\partial U_{i}^{\mathcal {P}}}{\partial x_{i}}=0, \quad \epsilon ^{4} Re \, U_{j}^{\mathcal {P}} \frac {\partial U_{i}^{\mathcal {P}}}{\partial x_{j}}=-\frac {\partial P^{\mathcal {P}}}{\partial x_{i}}+\epsilon \frac {\partial ^{2} U_{i}^{\mathcal {P}}}{\partial x_{j}^{2}} , \\[6pt] \nonumber \end{align}

with $Re= {\rho \,\mathcal {U} \mathcal {L}}/{\mu }$ . Note that the time scale is the same in the interface and free-fluid region ( $t = \hat {t} \, \mathcal {U}/\mathcal {L}$ ) and that in the bulk of the porous domain the motion is assumed steady. In the intermediate and porous regions, the dependent variables are function of both the fast (microscopic) and the slow (macroscopic) coordinates ( $x_{i}$ , $X_{i}$ , respectively), whilst in the channel-flow region, the dependent variables vary spatially with the macroscopic coordinates, $X_{i}$ , only. A fictitious dividing surface between the channel-flow region and the interfacial layer is defined at $\hat x_2= \hat {y}_\infty$ , and continuity of the velocity and the traction vectors is applied there. With $y_\infty =\hat {y}_\infty /\tilde \ell$ and $\mathcal {Y}_\infty = \hat {y}_\infty /\mathcal {L}=\epsilon y_\infty$ the microscopic and the macroscopic vertical coordinates of this interface, respectively, the matching conditions can be written as follows:

(2.4a) \begin{align} \lim _{x_2 \rightarrow y_\infty }\,U^{\mathcal {I}}_{i} & = \frac {1}{\epsilon } \lim _{X_2 \rightarrow \mathcal {Y}_\infty }\,U^{\mathcal {C}}_{i}, \end{align}

(2.4b) \begin{align} \displaystyle {\lim _{x_2 \rightarrow y_\infty }\,\left (-P^{\mathcal {I}}\delta _{i2}+\frac {\partial V^{\mathcal {I}}}{\partial x_i}+\frac {\partial U^{\mathcal {I}}_i}{\partial y}\right )} & = {\lim _{X_2 \rightarrow \mathcal {Y}_\infty }\,\left (-Re\,P^{\mathcal {C}}\delta _{i2}+\frac {\partial V^{\mathcal {C}}}{\partial X_i}+\frac {\partial U^{\mathcal {C}}_i}{\partial Y}\right )}, \\[6pt] \nonumber \end{align}

with $\delta _{ij}$ the Kronecker index. For the conditions above to be valid, $y_\infty$ must be sufficiently large such that the $\bullet ^{\mathcal {I}}$ variables become uniform in $x$ and $z$ at the virtual interface.

2.2. Asymptotic analysis of the microscale problem

The dependent variables in the interfacial and the porous subdomains are expanded in terms of $\epsilon$ as

(2.5) \begin{align} U_i^{\mathcal {I}}= U_{i}^{\mathcal {I}(0)}+\epsilon U_{i}^{\mathcal {I}(1)}+ \epsilon ^2 U_{i}^{\mathcal {I}(2)}+\ldots .,\quad P^{\mathcal {I}}= P^{\mathcal {I}(0)}+\epsilon P^{\mathcal {I}(1)}+ \epsilon ^2 P^{\mathcal {I}(2)}+\ldots ., \end{align}

(2.6) \begin{align} U_i^{\mathcal {P}}= U_{i}^{\mathcal {P}(0)}+\epsilon U_{i}^{\mathcal {P}(1)}+ \epsilon ^2 U_{i}^{\mathcal {P}(2)}+\ldots .,\quad P^{\mathcal {P}}= P^{\mathcal {P}(0)}+\epsilon P^{\mathcal {P}(1)}+ \epsilon ^2 P^{\mathcal {P}(2)}+\ldots .. \\[6pt] \nonumber \end{align}

Furthermore, the gradients are recast using the chain rule $\displaystyle \ ({\partial }/{\partial x_i}) \rightarrow \displaystyle \ ({\partial }/{\partial x_i})+\epsilon \,\ ({\partial }/{\partial X_i})$ . The asymptotic expressions are substituted into the equations governing the flow in the microscopic regions, for the microscale problems to be reconstructed at different orders of $\epsilon$ . It has been shown (Naqvi & Bottaro Reference Naqvi and Bottaro2021) that the resulting systems of equations for the interfacial and the porous regions can be combined by defining a composite description of the asymptotic expansions, that is

(2.7a) \begin{align} u_{i}= u_{i}^{(0)}+\epsilon u_{i}^{(1)}+ \mathcal{O}(\epsilon ^{2}), \quad p= p^{(0)}+\epsilon p^{(1)}+ \mathcal{O}(\epsilon ^{2}), \end{align}

with

(2.7b) \begin{align} u_{i}^{(0)}=\left \{\begin{array}{l} U_{i}^{\mathcal {I}(0)}, \quad y \in \mathcal {I} \\ \epsilon U_{i}^{\mathcal {P}(0)}, \quad y \in \mathcal {P} \end{array}\right ., \quad u_{i}^{(1)}=\left \{\begin{array}{l} U_{i}^{\mathcal {I}(1)}, \quad y \in \mathcal {I} \\ \epsilon U_{i}^{\mathcal {P}(1)}, \quad y \in \mathcal {P} \end{array}\right ., \end{align}

and

(2.7c) \begin{align} p^{(0)}=\left \{\begin{array}{l} P^{\mathcal {I}(0)}, \quad y \in \mathcal {I} \\ P^{\mathcal {P}(0)} + \epsilon P^{\mathcal {P}(1)}, \quad y \in \mathcal {P} \end{array}\right ., \quad p^{(1)}=\left \{\begin{array}{l} P^{\mathcal {I}(1)}, \quad y \in \mathcal {I} \\ \epsilon P^{\mathcal {P}(2)}, \quad y \in \mathcal {P} \end{array}\right .. \end{align}

The following composite system, valid over the whole region below the dividing interface (i.e. $x_2\lt y_\infty$ ), is thus obtained:

(2.8) \begin{align} \left \{\begin{array}{l} \partial _{i} u_{i}=-\epsilon \partial _{i}^{\prime } u_{i}^{(0)}+ \mathcal{O}\left (\epsilon ^{2}\right ) \\ -\partial _{i} p+ \partial _{j}^{2} u_{i}= \mathcal {R} \,u_{j}\partial _{j} u_{i}+ \epsilon \left [\partial _{i}^{\prime } p^{(0)}-2 \partial _{j} \partial _j^{\prime } u_{i}^{(0)}+ \mathcal {R} \, u_{j}^{(0)} \partial _{j}^{\prime } u_{i}^{(0)}\right ]+ \mathcal{O}\left (\epsilon ^{2}\right ) \end{array}\right . \end{align}

with $\mathcal {R}=\epsilon ^2 Re$ a microscopic Reynolds number and with derivatives indicated by $\displaystyle \partial _i=\ ({\partial }/{\partial x_i})$ and $\displaystyle \partial ^{\prime }_i=\ ({\partial }/{\partial X_i}).$

In order to treat the problem above, we first simplify it by linearizing the convective terms applying an Oseen approximation. In particular, a constant value is assigned to the streamwise velocity component, $u_1$ , near the interface, chosen as the surface-averaged slip velocity $\displaystyle u_{slip}= \ ({\hat u_{slip}}/{\epsilon \mathcal {U})}$ (with $\hat u_{slip}$ the dimensional slip velocity at the plane $\hat y=0$ ), i.e. $u_j^{(0)} \simeq (u_{slip},0,0)$ (other choices are clearly possible. For example, Bottaro (Reference Bottaro2019) tested the friction velocity as advective speed; results shown in the following support the present choice of the slip velocity). Thus, the advection term in (2.8) simplifies as $\mathcal {R}\,{ u_{slip}} \partial _1 u_i$ , with

(2.9) \begin{align} \mathcal {R}\,u_{slip}= \frac {\rho \, \hat {u}_{slip}\,\tilde \ell }{\mu }= Re_{slip}. \end{align}

The quantity $Re_{slip}$ is a slip-velocity Reynolds number, based on the microscopic length scale $\tilde \ell$ ; as we will see later, its value is not necessarily small. The composite system (2.8) is now approximated as

(2.10) \begin{align} \left \{\begin{array}{l} \partial _{i} u_{i}=-\epsilon \partial _{i}^{\prime } u_{i}^{(0)}+ \mathcal{O}\left (\epsilon ^{2}\right ) \\ -\partial _{i} p+ \partial _{j}^{2} u_{i}= Re_{slip} \, \partial _1 u_i + \epsilon \left [\partial _{i}^{\prime } p^{(0)}-2 \partial _{j} \partial _j^{\prime } u_{i}^{(0)}+ Re_{slip} \, \partial ^{\prime }_1 u_i^{(0)}\right ]+ \mathcal{O}\left (\epsilon ^{2}\right ) \end{array}\right . .\end{align}

The leading-order problem reads

(2.11) \begin{align} \mathcal{O}(1):\left \{\begin{array}{l} \partial _{i} u_{i}^{(0)}=0, \\ -\partial _{i} p^{(0)}+\partial _{j}^{2} u_{i}{ }^{(0)}=Re_{slip} \, \partial _{1} u_{i}^{(0)}, \\ \left (-p^{(0)} \delta _{i_{2}}+\partial _{2} u_{i}^{(0)}+\partial _{i} u_{2}^{(0)}\right )_{x_2=y_{\infty }}=S_{i2}^{\mathcal C}, \end{array}\right . \end{align}

with $S_{i2}^{\mathcal C}$ the macroscopic traction vector evaluated at $X_2={\mathcal {Y}_\infty }$ , i.e.

(2.12) \begin{align} S_{i2}^{\mathcal C}= \boldsymbol {\sigma ^{\mathcal C}} \cdot \boldsymbol {e}_2|_{X_2={\mathcal {Y}_\infty }} =\left .\left (\frac {\partial U^{\mathcal C}}{\partial Y}+\frac {\partial V^{\mathcal C}}{\partial X},-Re P^{\mathcal C}+2 \frac {\partial V^{\mathcal C}}{\partial Y},\frac {\partial W^{\mathcal C}}{\partial Y}+\frac {\partial V^{\mathcal C}}{\partial Z} \right )\right |_{X_2 = {\mathcal {Y}}_\infty }, \end{align}

where $\boldsymbol {\sigma }^{\mathcal C}$ is the stress tensor. From now on, the outer dependent variables are written without the superscript $\bullet ^{\mathcal C}$ .

At next order, we have

(2.13) \begin{align} \mathcal{O}(\epsilon ):\left \{\begin{array}{l} \partial _{i} u_{i}^{(1)}=-\partial _{i}^{\prime } u_{i}^{(0)}, \\[3pt] -\partial _{i} p^{(1)}+\partial _{j}^{2} u_{i}^{(1)}=Re_{slip}\left (\partial _{1} u_{i}^{(1)}+\partial ^{\prime }_{1} u_{i}^{(0)}\right )+\partial _{i}^{\prime } p^{(0)}-2 \partial _{j} \partial _{j}^{\prime } u_{i}^{(0)} , \\[5pt] \left (-p^{(1)} \delta _{i_{2}}+\partial _{2} u_{i}^{(1)}+\partial _{i} u_{2}^{(1)}\right )_{x_2=y_{\infty }}= -\left (\partial ^{\prime }_{2} u_{i}^{(0)}+\partial ^{\prime }_{i} u_{2}^{(0)}\right )_{x_2=y_{\infty }}. \end{array}\right . \end{align}

The linearity of (2.11) and (2.13) permits us to assume generic solutions of the problems. For the leading-order problem, we express the dependent variables as

(2.14) \begin{align} \left \{\begin{array}{l} u_{i}^{(0)}=u_{i j}^\dagger S_{j2}, \\[5pt] p^{(0)}=p^\dagger _{j} S_{j2}, \end{array}\right . \end{align}

with the closure variables, $u_{i j}^\dagger$ and $p^\dagger _{j}$ , dependent on only the microscopic coordinates, $x_i$ . Decoupled ad hoc auxiliary systems arise from plugging the generic solutions into (2.11); they are

(2.15) \begin{align} \left \{\begin{array}{l} \partial _{i} {u}^\dagger _{i j}=0, \\[3pt] -\partial _{i} {p}_{j}^\dagger +\partial _{l}^{2} {u}_{ij}^\dagger =Re_{slip}\, \partial _{1} {u}^\dagger _{ij},\\[3pt] \left .\left (-{p}_{j}^\dagger \delta _{i 2}+\partial _{2} {u}_{i j}^\dagger +\partial _{i} {u}_{2j}^\dagger \right )\right | _{x_2= y_{\infty }}=\delta _{ij}, \end{array}\right . \end{align}

where the microscopic problems correspond to $j=1,2,3$ . For the problem forced by $S_{22}$ (i.e. with $j=2$ ), the analytical solution ${u}_{i2}^\dagger =0, \, {p}_{2}^\dagger =-1$ is easily retrieved. At $O(\epsilon )$ the following generic forms hold:

(2.16) \begin{align} \left \{\begin{array}{l} u_{i}^{(1)}=u_{ijk}^{\ddagger } \,\partial ^{\prime }_{k} S_{j2}, \\[5pt] p^{(1)}=p_{jk}^{\ddagger } \,\partial _{k}^{\prime } S_{j2}, \end{array}\right . \end{align}

leading to

(2.17) \begin{align} \left \{\begin{array}{l} \partial _{i} u_{ijk}^{\ddagger }=-{u}_{kj}^\dagger , \\ Re_{slip}\left (\partial _{1} u_{ijk}^{\ddagger }+{u}^\dagger _{ij} \delta _{k1}\right )=-\partial _{i} p_{jk}^{\ddagger }-{p}^\dagger _{j} \delta _{ki}+\partial _{l}^{2} u_{ijk}^{\ddagger }+2 \partial _{k} {u}^\dagger _{ij},\\ \left .\left (-{p}^\ddagger _{jk} \delta _{i2}+\partial _{2} u^\ddagger _{ijk}+\partial _{i} u^\ddagger _{2jk}\right )\right |_{x_2= y_{\infty }}=-\left .\left ({u}_{ij}^\dagger \delta _{k2}+ {u}_{2j}^\dagger \delta _{ik}\right )\right |_{x_2=y_{\infty }}; \end{array}\right . \end{align}

these are nine decoupled systems, i.e. corresponding to $j,k=1,2,3$ . The closure problems (2.15) and (2.17) are to be solved in a representative unit cell of the microscopic region, subject to periodicity of all the dependent variables along $x$ and $z$ and to the boundary conditions $u_{ij}^\dagger = 0$ and $u_{ijk}^\ddagger =0$ on the solid grains, arising from the no-slip condition. Further, the microscopic unit cell is delimited from the bottom (theoretically at $y \to -\infty$ ) by the bulk of the porous domain, where dependent variables are cyclic of period 1 also along $y$ ; from a numerical perspective, results do not change provided the domain is at least two rows deep.

2.3. Formal expressions of the effective boundary conditions

Numerical solutions are sought for systems (2.15) and (2.17), with focus on the values of the fields at $x_2=y_\infty$ since $ {u}_{ij}^\dagger |_{y_{\infty }}$ and $ u^\ddagger _{ijk} |_{y_{\infty }}$ are eventually the coefficients needed to close the macroscopic effective boundary conditions for the velocity; these conditions result from matching the velocity vector at the fictitious interface between the channel-flow and the interfacial regions, as per (2.4a ). The upscaled conditions, second-order accurate in terms of $\epsilon$ , are

(2.18) \begin{align} \left .U_i\right |_{\mathcal {Y}_{\infty }}=\epsilon \left (\left .u_i^{(0)}\right |_{y_{\infty }}+\epsilon \left .u_i^{(1)}\right |_{y_{\infty }} \right ) + \mathcal{O}(\epsilon ^3)= \epsilon \left .{u}^\dagger _{ij}\right |_{y_{\infty }} S_{j2} +\epsilon ^2 \left .u^\ddagger _{ijk}\right |_{y_{\infty }}\, \frac {\partial S_{j2}}{\partial X_k} + \mathcal{O}(\epsilon ^3). \end{align}

The numerical procedure to solve the closure problems is similar to that followed by Naqvi & Bottaro (Reference Naqvi and Bottaro2021) and Ahmed et al. (Reference Ahmed, Naqvi, Buda and Bottaro2022b ) for porous media of either isotropic (such as spherical grains) or transversely isotropic microstructures in the $\hat {x}-\hat {z}$ plane (such as spanwise- or streamwise-elongated elements). We focus on the same parameters that do not vanish at the matching interface found in these references,

(2.19) \begin{align} \begin{array}{l} \left .{u}_{11}^\dagger \right |_{y_{\infty }} = y_{\infty } + \lambda _x, \quad \left .{u}_{33}^\dagger \right |_{y_{\infty }} = y_{\infty } + \lambda _z,\\ - \left .{u}_{211}^{\ddagger }\right |_{y_{\infty }}= \left .{u}_{121}^{\ddagger }\right |_{y_{\infty }} = 0.5\,y_{\infty }^2 + \lambda _x\,y_{\infty } + \mathcal {K}^{itf}_{xy}, \\ -\left .{u}_{233}^{\ddagger }\right |_{y_{\infty }}= \left .{u}_{323}^{\ddagger }\right |_{y_{\infty }} = 0.5\,y_{\infty }^2 + \lambda _z\,y_{\infty } + \mathcal {K}^{itf}_{zy}, \\ \left .{u}_{222}^{\ddagger }\right |_{y_{\infty }} = \mathcal {K}_{yy}, \end{array} \end{align}

with $\lambda _x$ and $\lambda _z$ the dimensionless Navier’s slip coefficients in the streamwise and the spanwise directions, respectively, $\mathcal {K}^{itf}_{xy}$ and $\mathcal {K}^{itf}_{zy}$ the interface permeability coefficients and $\mathcal {K}_{yy}$ an intrinsic permeability component. The novel contribution here is the incorporation of the effect of near-interface inertia on the microscale flow behaviour, which renders the aforementioned parameters sensitive to the value of $Re_{slip}$ .

Once relations (2.19) are plugged into (2.18) macroscopic matching conditions at the interface $Y_\infty = \epsilon y_\infty$ between the intermediate and the outer region are obtained. These conditions can then be transferred to $Y=0$ by a second-order Taylor expansion, to eventually yield the following effective boundary conditions:

(2.20a) \begin{align} \left .U\right |_{Y=0} =\left .\epsilon \lambda _{x} S_{12}\right |_{Y=0}+\left .\epsilon ^{2} \mathcal {K}_{x y}^{i t f} \frac {\partial S_{22}}{\partial X}\right |_{Y=0}+ \mathcal{O} \big(\epsilon ^{3}\big), \end{align}

(2.20b) \begin{align} \left .V\right |_{Y=0}=-\left .\epsilon ^{2} \mathcal {K}_{x y}^{i t f} \frac {\partial S_{12}}{\partial X}\right |_{Y=0}-\left .\epsilon ^{2} \mathcal {K}_{z y}^{i t f} \frac {\partial S_{32}}{\partial Z}\right |_{Y=0}+\left .\epsilon ^{2} \mathcal {K}_{yy} \frac {\partial S_{22}}{\partial Y}\right |_{Y=0}+ \mathcal{O}\big (\epsilon ^{3}\big ), \end{align}

(2.20c) \begin{align} \left .W\right |_{Y=0}=\left .\epsilon \lambda _{z} S_{32}\right |_{Y=0}+\left .\epsilon ^{2} \mathcal {K}_{z y}^{i t f} \frac {\partial S_{22}}{\partial Z}\right |_{Y=0}+ \mathcal{O} \big (\epsilon ^{3}\big). \end{align}

At this point, something should be said about the scales ( $\tilde \ell , \mathcal L, \mathcal U$ ) used to normalize the preceding equations, and on whether the principle of separation of scales is satisfied. If the magnitude of the macroscopic pressure gradient driving the flow in the channel is $\displaystyle \mathcal {M}=\left |{\Delta \hat p}/{\hat {L}_x}\right |$ , one may derive a stress $ \tau _{\mathcal {M}}= \mathcal {M}\, H= ({\tau _{\mathcal {B}}+\tau _{\mathcal {T}}})/{2}$ , with $\tau _{\mathcal {B}}$ and $\tau _{\mathcal {T}}$ the total shear stresses at $Y=0$ (bottom) and $Y=2$ (top), respectively. The corresponding shear velocity $ u_{\tau (\mathcal {M})}= \sqrt {{\tau _{\mathcal {M}}}/{\rho }}$ is chosen here as the macroscopic velocity scale, i.e. $\mathcal {U}=u_{\tau (\mathcal {M})}$ . This is an appropriate characterization of the velocity of near-wall eddies, since it is known that the root-mean-square (r.m.s.) of the fluctuating speed scales with the friction velocity. As far as the macroscopic length scale is concerned, it has been proposed first by Luchini (Reference Luchini and Desideri1996) that the important boundary condition for turbulence is that experienced by quasistreamwise vortices. The relevant length scale should thus be the vortex diameter which is around $20$ viscous units, i.e. $\mathcal L \sim \alpha \, \displaystyle \ ({\nu }{u_{\tau /(\mathcal M)}})$ , with $\alpha$ a constant close to 20. As far as the microlength scale is concerned, we observe that in the asymptotic analysis by Saffman (Reference Saffman1971) for the case of the flow over an isotropic porous substrate of permeability $K$ , the choice $\tilde \ell = \sqrt {K}$ was made. Recently, Hao & García-Mayoral (Reference Hao and García-Mayoral2025) have argued that the effect of deep porous substrates (those considered in the present paper can also be characterized as ‘deep’) on the near-wall flow is essentially governed by the intrinsic permeability of the medium, which, in all cases considered here, is a small fraction of $\ell ^2$ (cf. also Mei & Vernescu (Reference Mei and Vernescu2010); Zampogna et al. (Reference Zampogna, Magnaudet and Bottaro2019a ); Lācis et al. (Reference Lācis, Sudhakar, Pasche and Bagheri2020); Naqvi & Bottaro (Reference Naqvi and Bottaro2021); Sudhakar et al. (Reference Sudhakar, Lācis, Pasche and Bagheri2021)), with $\ell$ the periodicity of the grains. For regularly roughened walls, when a substrate permeability cannot be defined, the proper length scale should be a measure of the slip length of the texture, also a small fraction of the periodicity, $\ell$ . At this point, we anticipate that $\ell ^+$ is approximately 20 times larger than $\sqrt {\mathcal K_{yy}^+}$ (this will be shown later for all porous substrates considered.) Thus, it seems appropriate to say that $\tilde \ell \sim \ell /\alpha$ , with the same value of the constant $\alpha$ as in the definition of the macrolength $\mathcal L$ . We are now able to estimate the small parameter $\epsilon$ of the expansion; it is found that $\epsilon \sim \ell ^+/\alpha ^2$ ranges from around 0.025 up to 0.1 for cylindrical inclusions having $\ell ^+$ values between $10$ and $40$ . For the ‘modified’ grains described later, $\epsilon$ would be even smaller, considering that such inclusions tend to block the flow and the permeability is much lower than in the previous case.

On the basis of the arguments presented, microscopic and macroscopic length scales are sufficiently well separated, for all the cases treated in this paper. The microscopic length scale in our analysis is the displacement of the origin of the near-wall vortex, caused by the presence of either a rough or a porous substrate; the macroscopic scale is the diameter of the vortex itself. One referee of this work objected vigorously to our choice of microscopic length scale, arguing that the only proper microscale is the pattern periodicity. If, as they objected, this were indeed the case, then $\tilde \ell \sim \ell$ and the parameter of the expansion would become $\epsilon \sim \ell ^+/\alpha$ , which exceeds $1$ for $\ell ^+\gt 20$ . Beyond $\ell ^+ \approx 20$ , they argued, separation of scale, and the expansion proposed, would be untenable. In the end, we believe that only a posteriori verifications against feature-resolving results can inform on the domain of validity of the upscaling procedure adopted; this crucial point will be addressed in §§ 3.1 and 4.

In dimensional form, the effective boundary conditions are the same that have been found before (Naqvi & Bottaro Reference Naqvi and Bottaro2021) and read

(2.21a) \begin{align} \left .\left .\hat {u}\right |_{0} \approx \hat {\lambda }_x\left (\frac {\partial \hat {u}}{\partial \hat {y}}+\frac {\partial \hat {v}}{\partial \hat {x}}\right )\right |_{0}+\left .\frac {\hat {\mathcal {K}}_{xy}^{itf}}{\mu } \frac {\partial }{\partial \hat {x}}\left (-\hat {p}+2 \mu \frac {\partial \hat {v}}{\partial \hat {y}}\right )\right |_{0}, \end{align}

(2.21b) \begin{align} \left .\left .\hat {v}\right |_{0} \approx \frac {\hat {\mathcal {K}}_{yy}}{\mu } \frac {\partial }{\partial \hat {y}}\left (-\hat {p}+2 \mu \frac {\partial \hat {v}}{\partial \hat {y}}\right )\right |_{0}-\left .\hat {\mathcal {K}}^{i t f}_{xy} \frac {\partial }{\partial \hat {x}}\left (\frac {\partial \hat {u}}{\partial \hat {y}}+\frac {\partial \hat {v}}{\partial \hat {x}}\right )\right |_{0} -\left .\hat {\mathcal {K}}^{i t f}_{zy} \frac {\partial }{\partial \hat {z}}\left (\frac {\partial \hat {w}}{\partial \hat {y}}+\frac {\partial \hat {v}}{\partial \hat {z}}\right )\right |_{0}, \end{align}

(2.21c) \begin{align} \left .\left .\hat {w}\right |_{0} \approx \hat {\lambda }_z\left (\frac {\partial \hat {w}}{\partial \hat {y}}+\frac {\partial \hat {v}}{\partial \hat {z}}\right )\right |_{0}+\left .\frac {\hat {\mathcal {K}}_{zy}^{itf}}{\mu } \frac {\partial }{\partial \hat {z}}\left (-\hat {p}+2 \mu \frac {\partial \hat {v}}{\partial \hat {y}}\right )\right |_{0}. \end{align}

Further considerations, simplifications and implementation-related details concerning the condition (2.21b ) are given in appendix A. After having established the effective conditions which hold at $\hat y = 0$ , we can render them adimensional in the most convenient way. Thus, we now choose to scale the governing (2.1) for the free-fluid region, as well as the corresponding interface conditions (2.21a )–(2.21c ), by the use of geometric scales (cf. figure 1); this corresponds to setting $\mathcal L = H$ in equations (2.2a – 2.2c ). By the same token, the microscopic problems in the unit cells is rescaled with the periodicity of the pattern and this amounts to setting $\tilde \ell = \ell$ in equations (2.2a – 2.2c) so that, eventually, $\epsilon$ is defined by the ratio $\ell$ over $H$ , like in the laminar case (Naqvi & Bottaro Reference Naqvi and Bottaro2021; Ahmed & Bottaro Reference Ahmed and Bottaro2024). Having rescaled the problem for computational convenience, the dimensional model coefficients now read

(2.22) \begin{align} \hat {\lambda }_{x,\,z}=\lambda _{x,\,z}\,l, \quad \hat {\mathcal {K}}_{xy,\,zy}^{itf}=\mathcal {K}_{xy,\,zy}^{itf}\,l^2, \quad \hat {\mathcal {K}}_{yy}=\mathcal {K}_{yy}\,l^2. \end{align}

We want to emphasize that these coefficients are not empirical, but arise from the solution of auxiliary systems of equations solved in the $\hat x$ - or $\hat z$ -periodic elementary cell of figure. 1. The two terms, $\hat {\mathcal {K}}_{xy}^{itf}$ and $\hat {\mathcal {K}}_{zy}^{itf}$ , are interface permeabilities since, by analogy to Darcy’s law in the bulk of the porous domain, they multiply the streamwise and spanwise gradients of the pressure in the expressions of $\hat u|_0$ and $\hat w|_0$ .They differ from the corresponding intrinsic permeability components which come from the solution of Stokes problems in a triply periodic unit cell taken in the bulk of the porous region and, as such, have little in common with the flow around the porous–free-fluid interface (Bottaro Reference Bottaro2019).

2.4. Evaluation of the macroscopic coefficients for selected geometries

Typical geometries of the inclusions used to construct the porous media under study are illustrated in figure 2; they are aligned in either the streamwise direction (substrates $LC$ and $LM$ ) or the spanwise direction (substrates $TC$ and $TM$ ), all satisfying a porosity $\theta =0.5$ , where $\theta$ is defined by considering a cubic unit cell within the porous region and evaluating the ratio of the volume occupied by the fluid to the total volume of the cell.

Figure 2. The problems under study. The computational domain is displayed in (a), with the dimensions indicated in the macroscopic coordinates (normalized by half the channel height). (b) The bulk unit cell of the different porous media considered are drawn in microscopic dimensionless coordinates. All media have porosity $\theta = 0.5$ .

A simple method, similar to that followed by Ahmed et al. (Reference Ahmed, Naqvi, Buda and Bottaro2022b ), is used to numerically evaluate the macroscopic coefficients in the effective boundary conditions. First, the systems governing the microscale fields $u_{11}^\dagger$ and $u_{33}^\dagger$ are solved on a microscopic domain with a sufficiently large value of $y_\infty$ (like, for instance, the domain sketched in figure 3 a). In this work, the solution of the closure problems is conducted using finite-volume discretization, as by the implementation of Simcenter STAR-CCM + software; in general, the microscopic domain is discretized into polygonal/polyhedral cells with sufficient mesh refinement in close vicinity of the porous–free-fluid interface such that grid-independent results for the closure fields are eventually obtained. Second, the Navier’s slip coefficients ( $\lambda _x,\lambda _z$ ) are estimated by averaging ${u}_{11}^\dagger$ and ${u}_{33}^\dagger$ , respectively, over the plane $y=0$ . The numerical values of the interface permeability coefficients can be computed via the following volume integrals:

(2.23) \begin{align} \begin{array}{l} \mathcal {K}^{itf}_{xy}=\int _{\mathcal {V}_{0}} {u}_{11}^\dagger \,\mathrm {d}V, \\[3pt] \mathcal {K}^{itf}_{zy}=\int _{\mathcal {V}_{0}} {u}_{33}^\dagger \,\mathrm {d}V, \end{array} \end{align}

where $\mathcal {V}_{0}$ denotes the whole fluid’s volume in the elementary cell below the interface at $y=0$ . This renders the dimensionless Navier-slip and the interface permeability coefficients dependent, in general, on the geometry of the inclusions and the slip-velocity Reynolds number, $Re_{slip}$ , which appears in the microscopic auxiliary systems (note that, after rescaling, $Re_{slip}$ is now defined with $\ell$ as length scale).

On the other hand, the medium permeability $\mathcal {K}_{yy}$ is intrinsic to the geometry of the porous region, where the velocity level is much smaller than $u_{slip}$ and the inertial effects are thus negligible; $\mathcal {K}_{yy}$ can be estimated by solving a Stokes system on a triply periodic cell of the porous domain, imposing unit forcing along $y$ , and evaluating the superficial average of the corresponding microscopic field over that cell (Mei & Vernescu Reference Mei and Vernescu2010).

Figure 3. Contours of the microscopic variables $u_{11}^\dagger$ , $u_{21}^\dagger$ and $u_{33}^\dagger$ at (i) $Re_{slip} = 0$ and (ii) $Re_{slip} = 30$ , shown over an $x-y$ plane for the case of transverse cylinders of porosity $\theta =0.5$ . Close-ups of the contours near the fluid–porous interface are presented, while the typical domain considered in the simulations is shown in (a). Slip and permeability coefficients are independent of the value of $y_\infty$ , provided it is larger than 2.

Transverse ( $\hat z$ -elongated) and longitudinal ( $\hat x$ -elongated) inclusions allow for further simplification of the microscopic, auxiliary problems, by setting either $\partial /\partial x_3$ or $\partial /\partial x_1$ to zero, respectively, yielding two-dimensional systems of equations. For the case of spanwise-elongated inclusions, we get the following two systems of interest at leading order:

(2.24) \begin{align} \left \{\begin{array}{l} \partial _{1} {u}^\dagger _{11}+\partial _{2} {u}^\dagger _{21}=0, \\[3pt] -\partial _{1} {p}_{1}^\dagger +\partial _{1}^{2} {u}_{11}^\dagger +\partial _{2}^{2} {u}_{11}^\dagger =Re_{slip}\, \partial _{1} {u}^\dagger _{11},\\[3pt] -\partial _{2} {p}_{1}^\dagger +\partial _{1}^{2} {u}_{21}^\dagger +\partial _{2}^{2} {u}_{21}^\dagger =Re_{slip}\, \partial _{1} {u}^\dagger _{21}, \\[3pt] \left .\left (\partial _{2} {u}_{11}^\dagger +\partial _{1} {u}_{21}^\dagger \right )\right | _{x_2= y_{\infty }}=1, \\[3pt] \left .\left (-{p}_{1}^\dagger +2\partial _{2} {u}_{2 1}^\dagger \right )\right | _{x_2= y_{\infty }}=0, \end{array}\right . \end{align}

and

(2.25) \begin{align} \left \{\begin{array}{l} \partial _{1}^{2} {u}_{33}^\dagger +\partial _{2}^{2} {u}_{33}^\dagger =Re_{slip}\, \partial _{1} {u}^\dagger _{33}, \\[3pt] \left .\left (\partial _{2} {u}_{33}^\dagger \right )\right | _{x_2= y_{\infty }}=1. \end{array}\right . \end{align}

Figure 4. Behaviours of the homogenization model parameters. Panel (a) displays results of the closure problems for the Navier-slip and interface permeability coefficients as functions of $Re_{slip}$ for the porous substrates TC (solid lines) and TM (dashed lines). In panel (b), the linear relation (2.28) between $\lambda _x$ and $Re_{slip}$ is plotted (black lines) for four values of $\epsilon$ , fixing $Re_{\tau (\mathcal M)}=193$ , in order to evaluate $Re_{slip}$ at the intersection points.

The numerical solutions of the previous systems under Stokes conditions and at $ Re_{slip}=30$ are shown in figure 3 for the case of transverse cylindrical inclusions, while the dependence of the macroscopic coefficients on $Re_{slip}$ is displayed in figure 4(a) for the substrates $TC$ and $TM$ . A preliminary estimation of the value of the slip velocity can be obtained from the first-order term in the effective boundary condition of $U$ , equation (2.20a ), which may be recast in terms of the wall distance in viscous units ( $Y^+ = Y Re_{\tau (\mathcal {M})}$ ) and the mean velocity, already normalized by $u_{\tau (\mathcal {M})}$ and hence from now on indicated as ${\overline {U}}^+$ , as follows:

(2.26) \begin{align} \left .{\overline {U}}^+\right |_{Y=0}\approx \left .\lambda _x^+ \frac {\partial {\overline {U}}^+}{\partial Y^+}\right |_{Y=0}, \end{align}

where $\displaystyle \lambda _x^+= ({\rho u_{\tau (\mathcal {M})} \hat {\lambda }_x}/{\mu })=\epsilon Re_{\tau (\mathcal {M})} \lambda _x$ . Provided the roughness maintains a sufficiently small amplitude so that the velocity gradient $\displaystyle \left . \ ({\partial {\overline {U}}^+}/{\partial Y^+})\right |_{Y=0}$ at the virtual wall remains close to $1$ , (2.26) simplifies to a Dirichlet boundary condition, i.e.

(2.27) \begin{align} \left .{\overline {U}}^+\right |_{Y=0}=\frac {\hat {u}_{slip}}{u_{\tau (\mathcal {M})}}\approx \epsilon Re_{\tau (\mathcal {M})} \lambda _x, \end{align}

which means that the slip-velocity Reynolds number can be written as

(2.28) \begin{align} Re_{slip}=\frac {\rho \hat {u}_{slip} \ell }{\mu }= \frac {\rho u_{\tau (\mathcal {M})} \ell }{\mu } \,\epsilon Re_{\tau (\mathcal {M})} \lambda _x= \epsilon ^2 Re_{\tau (\mathcal {M})}^2 \,\lambda _x. \end{align}

With $\displaystyle \lambda _x= \ ({Re_{slip}}/{\epsilon ^2 Re_{\tau (\mathcal {M})}^2})$ , a linear relation between $\lambda _x$ and $Re_{slip}$ can be drawn for different values of $\epsilon$ , at the fixed value of the friction Reynolds number, $Re_{\tau (\mathcal {M})}=193$ ; cf. figure 4(b). The value of $Re_{slip}$ at the intersection point is evaluated as shown in figure 4(b), yielding as an immediate consequence all the macroscopic coefficients at this value; table 1 reports all coefficients for $\epsilon$ ranging from $0.05$ to $0.2$ . The reader is referred here to the work by Fairhall et al. (Reference Fairhall, Abderrahaman-Elena and García-Mayoral2019) where useful findings regarding the possible deviations of the slip lengths from the viscous predictions are presented for surface textures different from those considered here.

Table 1. Values of the macroscopic coefficients for the 16 porous substrates considered in the present study. For all patterns, the porosity is $\theta = 0.5$ and $Re_{\tau (\mathcal M)}=193$ , while $\epsilon = \ell /H$ is varied from $0.05$ (subscript $5$ ) to $0.2$ (subscript $20$ ).

Streamwise-elongated inclusions (substrates $LC$ and $LM$ ) represent a special case since inertial effects at the microscale level disappear as a consequence of setting $\partial /\partial x_1$ to $0$ in the auxiliary systems (Luchini et al. Reference Luchini, Manzo and Pozzi1991); as such, the macroscopic coefficients are independent of $Re_{slip}$ . The coefficients for these substrates are directly available by revisiting the results for $TC$ and $TM$ , at $Re_{slip}=0$ , and simply switching the streamwise and spanwise coordinates.

3.1. Validation of the model

The applicability of the upscaling approach followed is assessed here by considering the turbulent flow ( $Re_{\tau (\mathcal {M})}=193$ ) in a channel delimited from the bottom ( $Y\leq 0$ ) by the substrate $TC_{20}$ (transverse cylinders, $\epsilon =0.2$ ), and validating sample results of the homogenized simulation, based on the effective boundary conditions (2.20a – 2.20c ) with the macroscopic coefficients given in table 1, against a classical fine-grained simulation. The mesh requirements, and thus the numerical cost, of the latter are much higher since it needs to resolve the seepage flow in the bulk of the porous domain and to account for the interactions occurring across the interfacial region, where significant ejection and sweep events take place (cf. figure 5). Quantitatively, the number of finite-volume cells in the fully resolving DNS ( $N_{cells}\approx 7.3 \times 10^{6}$ ) is four times that in the homogenized DNS ( $N_{cells}\approx 1.8 \times 10^{6}$ ). Another point to be taken into account with respect to the full DNS is the technical complexity associated with the mesh generation, especially in the interfacial region. An unstructured grid was used in the porous substrate and in the lower-half of the channel; the cells are polygonal in section (on the $X-Y$ plane) and are extruded in the spanwise direction with a uniform spacing of $\approx 6$ viscous units. For the top row of cylinders (the closest to the porous–free-fluid interface), the mesh is refined such that the first cell centre is at a distance of around $0.3$ viscous units from the cylindrical grain, measured in the direction normal to the boundary.

Figure 5. Full feature-resolving simulation of the coupled flow problem including the flow through and the turbulent flow over the porous substrate $TC_{20}$ at $Re_{\tau (\mathcal {M})}=193$ : profiles of the $X$ – $Z$ -averaged mean velocity across the free-fluid region and closely below the fluid–porous interface are plotted. Instantaneous distributions (examples) of the interface-normal velocity component, $V^+$ , captured during ‘suction’ and ‘blowing’ events are also displayed.

Figure 6. Turbulent channel flow ( $Re_{\tau (\mathcal {M})}=193$ ) over the porous substrate $TC_{20}$ : predictions of the homogenized model, indicated by green lines with filled symbols, for (a,b) the mean velocity profile across the channel and for the near-interface distributions of (c) the r.m.s. of the turbulent fluctuations in the three velocity components, (d) the turbulence intensities and (e,f) the Reynolds/viscous shear stresses are validated against results of the full simulation (red lines). The dashed black profiles refer to the corresponding smooth, impermeable channel case.

The results in the free-fluid region are presented and compared (homogenization-based versus fine-grained) in figures 6 and 7, in terms of the following dimensionless parameters: the mean velocity, $\overline {U}^+$ ; the r.m.s. of the fluctuations in the velocity components, $(U_{rms},\,V_{rms},\,W_{rms})=(\overline {U^{\prime}U^{\prime}}^{1/2},\, \overline {V^{\prime}V^{\prime}}^{1/2},\, \overline {W^{\prime}W^{\prime}}^{1/2})$ where the turbulent fluctuations are defined as $U^{\prime}_i={U^+_i}-\overline {U}^+_i$ ; the intensity of the fluctuations, $\displaystyle (I_U,\,I_V,\,I_W)= (({U_{rms}}/{\overline {U}^+}),\,\ ({V_{rms}}/{\overline {U}^+}),\,\ ({W_{rms}}/{\overline {U}^+}) )$ ; the Reynolds shear stress, $\tau _{XY}^R=-\overline {U^{\prime}V^{\prime}}$ ; the viscous shear stress, $\displaystyle \tau _{XY}^V=\ ({1}/{Re_{\tau (\mathcal {M})}}) \ ({\partial \overline {U}^+}/{\partial Y})$ ; and the production rate of the turbulent kinetic energy, $\displaystyle P_T= \ ({-1}/{Re_{\tau (\mathcal {M})}}) \overline {U^{\prime}_i U^{\prime}_j}\ ({\partial \overline {U}_i^+}/{\partial X_j})$ . While figure 6 focuses on the validation of the present model with the effective boundary conditions of the three velocity components imposed at $Y=0$ , figure 7 shows, in addition, the corresponding macroscopic results when the interface-normal velocity component is suppressed (i.e. $V|_{Y=0}=0$ ) and only the in-plane slip velocities are applied; this is important since it highlights the need of accounting for transpiration at the virtual boundary (Gómez-de-Segura et al. Reference Gómez-de-Segura, Sharma and García-Mayoral2018a ; Bottaro Reference Bottaro2019; Lācis et al. Reference Lācis, Sudhakar, Pasche and Bagheri2020).

From inspection of figure 6, it is clear that the model captures well the trends of the mean velocity and the turbulence statistics displayed. The velocity profile can be analysed in terms of the slip velocity, $U_{slip}^+= \overline {U}^{+}\big |_{Y=0}$ ; the shift in the intercept of the logarithmic velocity profile, $\Delta U^+$ (taking the smooth channel case as a reference for the measurement and averaging the shift over the region $30 \lesssim Y^+ \lesssim 120$ (Ahmed et al. Reference Ahmed, Naqvi, Buda and Bottaro2022b )); the percentage change in the bulk (channel-averaged) velocity through the free-fluid region, $\Delta U_{ch}^+ \, \%$ (taking the bulk velocity in a fully smooth channel, $U_{ch}^+ \approx 15.69$ , as a reference); and the corresponding percentage change in skin-friction coefficient, $\Delta C_f \, \%$ (taking the smooth-channel value, $C_f= 2/(U_{ch}^+)^2 \approx 0.00813$ , as a reference). The analysis performed here shows that, for the turbulent flow over the perturbed boundaries considered, the log-law is still valid (over $30 \lesssim Y^+ \lesssim 120$ ), yet it is shifted (relative to that for a smooth wall) by $\Delta U^+$ such that the logarithmic profile reads

(3.1) \begin{align} \overline {U}^+ = \frac {1}{\kappa }\,{\textrm {ln}}\big(Y^+ \big) + B + \Delta U^+, \end{align}

Figure 7. Distribution of the mean velocity (a) and behaviours of quantities of interest related to turbulence statistics (b–f) over the porous substrate $TC_{20}$ : predictions of the homogenized simulation when the effective boundary conditions of the three velocity components are imposed (green lines with filled circles) or when transpiration is neglected (blue lines) are validated against results of the fine-grained simulation (red lines), while the dashed profiles are related to the smooth, impermeable channel case.

where $\kappa$ is the von Kármán constant and $B$ is the intercept of the logarithmic profile for the flow over a corresponding smooth wall. Based on (3.1), if $\Delta U^+ \lt 0$ (respectively, $\Delta U^+ \gt 0$ ), the logarithmic profile is shifted downwards (upwards), and in general the skin-friction drag increases (decreases); this is consistent with the definition of the roughness function, $\Delta U^+$ adopted by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019),Ibrahim et al. (Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021) and Khorasani et al. (Reference Khorasani, Lācis, Pasche, Rosti and Bagheri2022, Reference Khorasani, Luhar and Bagheri2024), which differs in sign from that originally introduced by Hama (Reference Hama1954) and Clauser (Reference Clauser1954). According to (3.1), $\Delta U^+ \lt 0$ is generally accompanied by $\Delta U_{ch}^+ \, \% \lt 0$ and $\Delta C_f \, \% \gt 0$ . The full feature-resolving simulation for the case chosen for validation ( $TC_{20}$ ) yields $\displaystyle (U_{slip}^+, \Delta U^+, \Delta U_{ch}^+ \, \%, \Delta C_f \, \%) \approx (1.37, -2.76, -12.2\, \%, +29.6\, \%)$ , while the values obtained from the model are, respectively, $(1.44, -2.33,-10.3\, \% , + 24.4\, \%)$ ; a decrease in flow rate and, therefore, an increase in skin-friction coefficient is realized in both simulations. Moreover, the model predictions for the r.m.s. fluctuations of the velocity components at the fictitious interface ( $Y=0$ ) match well the results of the full simulation and deviate significantly from zero. In general, the accuracy of the macroscopic model is reasonable taking into account that the value of $\epsilon = 0.2$ ( $\ell ^+ \approx 40$ ) related to the porous substrate chosen for validation ( $TC_{20}$ ) is rather large, meaning that microscopic and macroscopic length scales do not differ widely. On the other hand, it is obvious from figure 7 that the comparison with the fine-grained simulation is not satisfactory when the transpiration-free model is applied, where the mechanism of drag increase is idle, $\Delta U^+$ is close to $0$ , and the trend of the turbulence statistics next to the fictitious boundary is similar to that of a smooth, impermeable channel (cf. Ibrahim et al. (Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021)). A third, important, scenario is presented later in § 4, where the transpiration velocity boundary condition is imposed while ${\mathcal {K}}_{yy}$ is set to zero, for the porous substrate to be modelled as a rough, impermeable wall. The discussion there centres around the evaluation of ${\mathcal {K}}_{yy}$ for a porous bed which is bounded from the bottom, and on how the accuracy of the model is affected by neglecting the medium permeability of a deep, yet finite, substrate such as the one considered here for validation (figure 5).

Figure 8. From (i) to (iii), instantaneous distributions of $U^{\prime}$ , $V^{\prime}$ and $W^{\prime}$ at the porous–free-fluid interface ( $Y = 0$ ) for case $TC_{20}$ . The fully resolved results (a) are compared with the homogenized ones (b).

Finally, figure 8 displays a comparison between the results of the full texture-resolving simulation and the homogenized one, concerning the fluctuating patterns of the three velocity components at the porous–free-fluid interface. This figure is added following one referee’s advice, with the purpose of providing the readers with the information needed to assess on their own how well the upscaled boundary conditions mimic the effect of the porous substrate on the turbulence.

3.2. Case studies: results and discussion

Numerical simulations were run for the channel flow over the different porous substrates (TC, LC, TM, LM), with four values of $\epsilon$ tested for each ( $\epsilon = 0.05, 0.1, 0.15, 0.2$ ). The macroscopic model already validated was employed to study the 16 problems under consideration, with the Oseen-based upscaled coefficients contributing to the effective boundary conditions available in table 1. The main results are presented and discussed below.

3.2.1. Mean velocity and skin-friction drag

In table 2, the pitch distance and the macroscopic coefficients for each porous pattern are expressed in wall units based on the velocity scale $u_{\tau (\mathcal {M})}$ ; they are defined by

(3.2a) \begin{align} \ell ^+ = \frac {\rho \,u_{\tau (\mathcal {M})}\,\ell }{\mu } = \epsilon Re_{\tau (\mathcal {M})}, \nonumber \end{align}

(3.2b) \begin{align} \lambda _x^+ = \frac {\rho \,u_{\tau (\mathcal {M})}\,\hat \lambda _x}{\mu } = \epsilon Re_{\tau (\mathcal {M})} \lambda _x, \quad \lambda _z^+ = \epsilon Re_{\tau (\mathcal {M})} \lambda _z, \end{align}

(3.2c) \begin{align} \mathcal {K}_{xy}^{itf,+} = \epsilon ^2 Re_{\tau (\mathcal {M})}^2 {\mathcal {K}}_{xy}^{itf}, \quad \mathcal {K}_{zy}^{itf,+} = \epsilon ^2 Re_{\tau (\mathcal {M})}^2 {\mathcal {K}}_{zy}^{itf}, \quad \mathcal {K}_{yy}^{+} = \epsilon ^2 Re_{\tau (\mathcal {M})}^2 {\mathcal {K}}_{yy}. \end{align}

Values of the major quantities related to the behaviour of mean velocity through the free-fluid region are also listed in the table (refer to the definitions in § 3.1). The most significant finding is that reduction of the skin-friction drag coefficient (negative values of $\Delta C_f\, \%$ , associated with positive $\Delta U^+$ and $\Delta U_{ch}^+\, \%$ ) is attainable only by the porous substrates formed by longitudinal inclusions (LC and LM), those characterized by streamwise-preferential slip lengths and interface permeabilities ( $\lambda _x^+ \gt \lambda _z^+$ , $\mathcal {K}_{xy}^{itf,+}\gt \mathcal {K}_{zy}^{itf,+}$ ). Such a favourable influence (up to $5\, \%$ reduction in $C_f$ ) takes place exclusively at relatively small values of $\ell ^+$ , a behaviour similar to that found by Gómez-de-Segura & García-Mayoral (Reference Gómez-de-Segura and García-Mayoral2019) for this kind of permeable boundaries and analogous to that exhibited by riblets (Bechert & Bartenwerfer Reference Bechert and Bartenwerfer1989; Garcia-Mayoral & Jiménez Reference Garcia-Mayoral and Jiménez2011a ; Endrikat et al. Reference Endrikat, Modesti, García-Mayoral, Hutchins and Chung2021a ,Reference Endrikat, Modesti, MacDonald, García-Mayoral, Hutchins and Chung b ; Wong et al. Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024). On the other hand, permeable beds consisting of transverse grains yield only drag increase, and this becomes more pronounced with $\ell ^+$ . For comparison purposes, the results obtained by normalizing results with the shear velocity of the bottom surface, $u_{\tau (\mathcal {B})}$ , are given in appendix B.

Table 2. Values of the macroscopic coefficients characterizing the different configurations considered for the porous substrate, estimated in wall units with $Re_{\tau (\mathcal {M})}=193$ and $\epsilon$ varied from $0.05$ (subscript $5$ ) to $0.2$ (subscript $20$ ). Major results are presented, with the normalization based on $u_{\tau (\mathcal {M})}$ . Monitoring the progress of the mean bulk velocity $U_{ch}^+$ during $10$ additional units of time, $\Delta U_{ch}^+\,\%$ is found to differ by $\pm 0.2\,\%$ at the most (and $\pm 0.07\,\%$ on average) from the final values reported in the table.

The behaviour of the sample quantities reported on the right-hand side of table 2 for the four substrate configurations are graphically presented as function of the streamwise Navier-slip lengths in figure 9. It is important to highlight the following features with reference to the trends of the figure.

Figure 9. Dependence of (a) the slip velocity $U_{slip}^+$ on $\lambda _x^+$ and of (b) the shift of the logarithmic profile intercept $\Delta U^+$ , (c) the percentage change in the bulk mean velocity $\Delta U_{ch}^+ \, \%$ and (d) the percentage change in the skin-friction coefficient $\Delta C_{f} \, \%$ on $\Delta \lambda ^+ =\lambda _x^+ - \lambda _z^+$ , for turbulent channel flows ( $Re_{\tau (\mathcal {M})}=193$ ) over the four types of permeable beds under study; cf. table 2. Simple linear relations fitting the behaviour of $U_{slip}^+$ with $\lambda _x^+$ and the performance of the other quantities at small values of $\Delta \lambda ^+$ are presented.

1. (i) As discussed in § 2.4, the first-order term in the effective boundary condition of the streamwise velocity yields a slip velocity at the permeable interface $\displaystyle U_{slip}^+\approx \lambda _x^+ \ ({\partial \overline {U}^+}/{\partial Y^+}) |_{Y=0} \approx \lambda _x^+$ . Figure 9(a) shows that this linear dependence fits well with the results of the simulations, for the roughness amplitudes considered. Besides the omission of the higher-order term, the small percentage errors (up to $\approx 11\, \%$ in absolute value) may be attributed to the deviation of $\displaystyle \ ({\partial \overline {U}^+}/{\partial Y^+}) |_{Y=0}$ from $1$ because the parameters are expressed in wall units based on $u_{\tau (\mathcal {M})}$ , and not the permeable-interface shear velocity $u_{\tau (\mathcal {B})}$ . However, one can write $\displaystyle \ ({\partial \overline {U}^+}/{\partial Y^+}) |_{Y=0} = ({\partial \overline {U}^{+ (\mathcal {B})}}/{\partial Y^{+ (\mathcal {B})}}) |_{Y=0} \, [\ ({u_{\tau (\mathcal {B})}}/{u_{\tau (\mathcal {M})}}) ]^2$ , where $\displaystyle \ ({\partial \overline {U}^{+ (\mathcal {B})}}/{\partial Y^{+ (\mathcal {B})}} ) |_{Y=0} \approx 1$ (provided that the Reynolds stress at $Y=0$ is much smaller than the viscous stress), which results in the modified relation $\displaystyle U_{slip}^+ \approx [({u_{\tau (\mathcal {B})}}/{u_{\tau (\mathcal {M})}}) ]^2 \, \lambda _x^+$ . This expression enhances the predictions of the slip velocity (maximum error below $4\, \%$ ), yet it cannot be employed a priori, since the values of the shear-velocity ratio are available only after numerical simulations have been conducted (cf. appendix B).

2. (ii) It has been found convenient to express the roughness function as the difference between the shifts of the virtual origins of mean and turbulent flows, i.e. $\displaystyle \Delta U^+ \approx \ell _U^+- \ell _{Turb.}^+$ (Ibrahim et al. Reference Ibrahim, Gómez-de-Segura, Chung and García-Mayoral2021). For small protrusion heights, Luchini et al. (Reference Luchini, Manzo and Pozzi1991) have shown that $\Delta U^+$ takes the form $\Delta U^+ = \lambda _x^+ - \lambda _z^+= \Delta \lambda ^+$ ; in the present settings, this assumption holds only up to $|\Delta \lambda ^+|\lesssim 0.25$ ; cf. figure 9(b). Also Gómez-de-Segura et al. (Reference Gómez-de-Segura, Sharma, García-Mayoral, Moin and Urzay2018b ) plotted the roughness function against the difference between the displacements of the virtual origins, and they did it for a variety of complex surfaces (anisotropic porous substrates, superhydrophobic surfaces, riblets and canopies), highlighting a behaviour (figure 5 in their paper) qualitatively similar to that displayed in figure 9(b).

3. (iii) As far as the trends of $\displaystyle \Delta C_f\, \%= \ ({C_f - C_{f,\,smooth}}/{C_{f,\,smooth}}) \times 100\, \%$ are concerned, the classical linearized relation $\displaystyle \Delta C_f\, \% = (({-\Delta U^+}/ {(2C_{f,\,smooth})^{-0.5}+(2\kappa )^{-1}} ) \times 100\, \%$ is expected to align well with the results for small changes in $C_f$ (Luchini Reference Luchini and Desideri1996; Bechert et al. Reference Bechert, Bruse, Hage, Van der Hoeven and Hoppe1997). With $\Delta U^+ = \Delta \lambda ^+$ , $C_{f,\,smooth}= 0.00813$ and the von Kármán constant $\kappa = 0.4$ , the linear dependence $\displaystyle \Delta C_f\, \% \approx -0.1 \,\Delta \lambda ^+ \times 100\, \%$ is valid provided that $|\Delta \lambda ^+|$ remains sufficiently small, as confirmed in figure 9(d). Under the same condition, it can be shown that $\Delta U_{ch}^+\, \% \approx -0.5 \Delta C_f\, \% \approx +0.05\,\Delta \lambda ^+ \times 100\, \%$ , which fits well the results in figure 9(c).

4. (iv) It is notable that, at any fixed value of $\Delta \lambda ^+$ , the porous substrates LM and TM outperform the configurations LC and TC in terms of either maximizing the drag reduction or minimizing the drag increase. One possible justification is that the permeable beds constructed with modified cylinders (LM and TM) exhibit much smaller values of the medium permeability $\mathcal {K}_{yy}^+$ compared with those designed based on flat cylinders (LC and TC), as can be realized from table 2. This favourable feature enhances $\Delta U^+$ by attenuating the transpiration velocity at the fictitious interface ( $Y=0$ ), an effect which can be perceived as a mitigation of the blowing and suction events. The influence of transpiration on $\Delta U^+$ will be discussed in further detail in § 3.3.

Figure 10. Dependence of major quantities characterizing the turbulent channel flow over the porous substrates under study on the pitch distance of the inclusions measured in wall units, $\ell ^+=\epsilon Re_{\tau (\mathcal {M})}$ . Results of the homogenization-based DNS are plotted with filled circles. The lines in panel (a) represent the simple relation $U_{slip}^+ = \ell ^+ \lambda _x$ , while those in the other panels are simple fitting curves.

In figure 10, the different results are plotted against the pitch distance, $\ell ^+$ . With regard to the slip velocity, assuming the simple linear relation $\displaystyle U_{slip}^+ \approx \lambda _x^+ = \ell ^+ \lambda _x$ and recalling the trends of $\lambda _x$ from table 1, one can expect that $U_{slip}^+$ changes linearly with $\ell ^+$ for the porous beds LC and LM since $\lambda _x$ is independent of $\ell ^+$ for these streamwise-elongated patterns, in contrast to the spanwise-elongated patterns TC and TM for which the coefficient $\lambda _x$ decreases with the increase of $\ell ^+$ on account of near-interface advection. These expectations agree with the behaviours displayed in figure 10(a). For small $\ell ^+$ values, the quantities $\Delta U^+$ , $\Delta U_{ch}^+\, \%$ and $\Delta C_f\, \%$ are directly proportional to $\Delta \lambda ^+= \ell ^+ \Delta \lambda$ (Luchini Reference Luchini and Desideri1996), where $\Delta \lambda$ is equal to $\lambda _x-\lambda _z$ . Table 1 implies that $\Delta \lambda \big |_{LM}\gt \Delta \lambda \big |_{LC}\gt \Delta \lambda \big |_{TC}\gt \Delta \lambda \big |_{TM}$ with the first two positive and the last two negative. For a small value of $\ell ^+$ , one should therefore expect $\Delta U^+ \big |_{LM}\gt \Delta U^+ \big |_{LC}\gt \Delta U^+ \big |_{TC}\gt \Delta U^+ \big |_{TM}$ (and likewise for $\Delta U_{ch}^+\, \%$ and $-\Delta C_f\, \%$ ) with the substrate LM yielding the maximum drag reduction and TM resulting in the maximum drag increase; cf. figure 10(b–d). Departing from the viscous regime, it is found that the drag reduction attainable by LM and LC peaks at some value of $\ell ^+$ between $10$ and $20$ . The performance of these porous substrates then degrades, yet drag reduction is still achievable until a threshold within $20 \lesssim \ell ^+ \lesssim 30$ is reached, beyond which drag increase takes place. Gómez-de-Segura et al. (Reference Gómez-de-Segura, Sharma and García-Mayoral2018a ) studied highly connected porous media with streamwise-preferential permeability and attributed the aforementioned behaviour to the formation of drag-increasing spanwise-coherent rollers associated with a Kelvin–Helmholtz-like instability whose initiation is governed by the intrinsic permeability component $\mathcal {K}_{yy}^+$ of the medium. Finally, we should not forget that all results plotted in figures 9 and 10 are obtained via homogenization-based DNS. The accuracy of the trends displayed is thus dependent on the accuracy of the upscaling approach for each of the cases considered (in this respect, one may want to go back to the validation conducted considering the pattern $TC_{20}$ in § 3.1).

3.2.2. The mechanism of drag increase/reduction

The influence of porous substrates on the near-interface turbulence is considered next. For sufficiently small values of $\ell ^+$ , the wall texture alters the structure of turbulence merely by shifting down its virtual origin by a distance $\ell _{Turb.}^+ \approx \lambda _z^+$ , whereas the effect is much more complicated beyond the viscous regime, especially with the increase in transpiration velocity. It is therefore useful to present and discuss some turbulence statistics of interest for the channel flow over selected porous beds of relatively large grain spacings/sizes ( $LM_{10}$ , longitudinal modified inclusions, $\ell ^+ \approx 20$ , the maximum drag reduction reported; $LM_{20}$ , longitudinal modified inclusions, $\ell ^+ \approx 40$ , drag increase; $TM_{10}$ , transverse modified inclusions, $\ell ^+ \approx 20$ , drag increase; $TM_{20}$ , transverse modified inclusions, $\ell ^+ \approx 40$ , the maximum drag increase reported). The velocity profiles are plotted in figure 11(a). The turbulence-characterizing quantities plotted in figure 11(b) are chosen since, as shown by Ahmed et al. (Reference Ahmed, Naqvi, Buda and Bottaro2022b ), their behaviours near the porous–free-fluid interface can be linked to the favourable/adverse effects of the permeable boundaries on friction drag. With focus on the peak values of $V_{rms}$ , $W_{rms}$ , $\tau _{XY}^R$ and $P_T$ , and the distributions of $I_W$ , it can be realized that the drag-reducing substrate ( $LM_{10}$ ) yields results comparable to those in the reference case of turbulence over a smooth, impermeable wall; this applies also to the other drag-reducing patterns not considered in the figure, i.e. $LM_{5}$ , $LC_{5}$ , $LC_{10}$ . Conversely, the drag-increasing ones result in intensified levels of these quantities. For instance, with $TM_{20}$ , the peak values $V_{rms}$ , $\tau _{XY}^R$ and $P_T$ are larger than the values in a smooth channel by approximately $16 \, \%$ , $24 \, \%$ and $50 \, \%$ , respectively. The values of the quantities at the fictitious interface, $Y=0$ , are of particular interest in the present work and their correlations with $\Delta U^+$ are explored in § 3.3; it is evident from the figure that significant values of $V_{rms}$ are obtained at the plane $Y=0$ , in particular when $\ell ^+$ is sufficiently large, an important effect (Jiménez et al. Reference Jiménez, Uhlmann, Pinelli and Kawahara2001; Orlandi et al. Reference Orlandi, Leonardi, Tuzi and Antonia2003, Reference Orlandi, Leonardi and Antonia2006; Orlandi & Leonardi Reference Orlandi and Leonardi2006, Reference Orlandi and Leonardi2008) which would obviously be absent if transpiration were unaccounted for in the formulation of the model.

Figure 11. Predictions of the homogenization-based model for (a) the mean velocity profiles and (b–f) sample statistics for the channel flow ( $Re_{\tau (\mathcal {M})} = 193$ ) over four different porous substrates.

Figure 12. Quadrant analysis of the Reynolds shear stress, $\tau _{xy}^R$ , for turbulent channel flows ( $Re_{\tau (\mathcal {M})} = 193$ ) over two different porous substrates ( $LM_{10}$ and $TM_{20}$ ). Instantaneous values of ( $U^{\prime}$ , $V^{\prime}$ ) throughout the planes at $Y = 0.005$ and $Y=0.416$ (evaluated at all grid points) are shown in panels (a) to (d), while contributions to $\tau _{xy}^R$ from each quadrant are plotted in the bottom frames against $Y^+ = Y Re_{\tau (\mathcal {M})}$ up to the centreline of the channel.

The quadrant analysis in figure 12 reveals details of the generation of the Reynolds stress, $\tau _{XY}^R$ , from the turbulent events taking place in the flow near the substrates $LM_{10}$ and $TM_{20}$ . In figure 12, the instantaneous distributions of $(U^{\prime},V^{\prime})$ are displayed over the plane at $Y^+ \approx 1$ , directly adjacent to the substrate–channel interface, and the plane at $Y^+ \approx 80$ , well above the substrate. The phenomena can be classified into negative-production events (first and third quadrants, with $-U^{\prime}V^{\prime}\lt 0$ ) and positive-production ones (second and fourth quadrants, with $-U^{\prime}V^{\prime}\gt 0$ ); refer to, for instance, Wallace et al. (Reference Wallace, Eckelmann and Brodkey1972). Eventually, the Reynolds stress generated from the sum of the positive contributions from the ejection (second quadrant, bursting of low-speed fluid) and the sweep (fourth quadrant, inrush of high-speed fluid) events at any $Y^+$ level is generally larger than that arising from the sum of the contributions of the other two quadrants. The production of turbulence is dominated by the sweep event in the close vicinity of the boundary (cf. figure 12 a,b), while ejection is dominant away from the wall (cf. figure 12 c,d). For a better understanding, the contributions from the four quadrants to the Reynolds shear stress at a given time instant, evaluated over different $X-Z$ planes up to $Y^+ \approx 80$ , are plotted in figure 12(e,f); they are obtained by integrating the values of $-U^{\prime}V^{\prime}$ related to each of the quadrants, separately, over the area occupied by the specific event and using the overall area of the $X-Z$ plane ( $=2\pi \times \pi$ ) as a weight. It is notable that ejection becomes dominant beyond a threshold within $Y^+= 12$ – $15$ . All the findings above agree qualitatively with the results by Kim et al. (Reference Kim, Moin and Moser1987) in a channel delimited by smooth, impermeable walls. From a quantitative perspective, the production of turbulence via both ejection and sweep is clearly intensified for case $TM_{20}$ (the porous substrate of maximum drag increase) compared with the levels with $LM_{10}$ (the substrate of largest drag reduction), at all the values of $Y^+$ considered.

In addition to the material presented in this section, it is beneficial to provide, via visualizations, some qualitative insights into the effects of the surface texture on the coherent structures (e.g. the pattern of streaks) and the turbulent events occurring in the inner region of the boundary layer; this is available in the Supplemental movie published online alongside this article. It might also be of interest, for future research focussed on the flow physics, to explore how different substrate topologies affect the spectral density of the Reynolds stress and the premultiplied spectra of the velocity components next to the interface. This would probably need also a more extensive comparison between texture-resolving and modelled simulations, considering a large variety of porous microstructures. For these reasons, in the present contribution we prefer to address attention to the relation between the roughness function and the upscaled coefficients of the model.

3.3. In pursuit of a correlation for $\Delta U^+$ over porous/textured walls

Figure 13. Dependence of $\Delta U^+$ (a–c) and the related quantities $\mathcal {D}$ (d–f) and $\mathcal {F}$ (g–i) on turbulence-characterizing parameters of interest measured at the fictitious interface (at $Y=0$ ). The filled symbols indicate results of the homogenized simulations for turbulent flow over the four substrate configurations under study (cf. figure 2), with $\ell ^+$ varied for each pattern as described in table 2, while the fitting relations (3.3–3.8) are plotted with solid lines in (d–i).

We proceed from the earlier discussion on figure 11, concerning how the near-wall distributions of some turbulence-characterizing parameters can control the mechanism of drag reduction/increase over the permeable boundaries, to explore the correlation between the roughness function $\Delta U^+$ and the fictitious-interface values of quantities of particular interest: $\displaystyle \tilde V_{rms}=V_{rms}\big |_{Y=0}$ , $\displaystyle \tilde W_{rms}=W_{rms}\big |_{Y=0}$ and $\displaystyle \tilde \tau _{XY}^R= \tau _{XY}^{R}\big |_{Y=0}$ . Figure 13(a–c) reveals that the dependence of $\Delta U^+$ on $\displaystyle \tilde V_{rms}$ cannot be described by a universal function valid for all permeable boundaries; the same can be said for $\displaystyle \tilde W_{rms}$ and $\displaystyle \tilde \tau _{XY}^R$ . Conversely, each configuration yields a unique relationship, and even the general trends differ when porous substrates of streamwise-preferential permeability (LC and LM, non-monotonic behaviour) are compared with those consisting of spanwise-elongated grains (TC and TM, strictly monotonic decrease). To explain this, let us assume conditions corresponding to a small value of $\tilde V_{rms}$ fixed for the four patterns (e.g. $\tilde V_{rms}= 0.01$ ) and analyse the resulting $\Delta U^+$ . While a fixed value of $\tilde V_{rms}$ may imply that, for all the boundaries, the virtual origin of turbulence has the same shift from the $Y=0$ plane (i.e. constant $\ell _{Turb.}^+$ ), the position of the virtual origin of the mean flow, $\ell _U^+\approx \lambda _x^+$ , can significantly differ according to the value of the streamwise Navier-slip length, $\lambda _x^+$ , for each wall, and, consequently, different values of the roughness function $\Delta U^+ = \ell _U^+ - \ell _{Turb.}^+$ are obtained. In the search of a function displaying a universal behaviour, we follow two separate paths.

The first path relies on analysing the mean velocity profile, $\overline {U}^+(Y^+)$ , over each of the permeable substrates to monitor the upward shifts of the velocity at matched $Y^+$ values, taking the profile over a smooth, impermeable wall as a reference (for instance, cf. figure 11 a). Such a velocity shift is, by definition, equal to $U_{slip}^+$ at $Y=0$ and to $\Delta U^+$ in the logarithmic region. Whether $\Delta U^+$ is positive or negative, it is $U_{slip}^+ \gt \Delta U^+$ for all textured boundaries (for $\ell ^+ = 0$ , the smooth, impermeable wall is retrieved, and the limit $U_{slip}^+ = \Delta U^+=0$ is reached). The function $\mathcal {D}=U_{slip}^+ - \Delta U^+ \geq 0$ can therefore be defined to indicate the depression in the velocity shift when moving from the wall to the logarithmic region; it is plotted against the turbulence parameters in figure 13(d–f).

In the second path, we proceed from the fact that the approximation $\Delta U^+ = \Delta \lambda ^+$ holds only for small surface roughness, while a further reduction in the value of the roughness function occurs with the increase of $\ell ^+$ , i.e. $\Delta U^+ = \Delta \lambda ^+ - \mathcal {F}$ , with the newly defined function $\mathcal {F} \geq 0$ . The behaviour of $\mathcal {F}$ is shown in figure 13(g–i).

Both functions $\mathcal {D}$ and $\mathcal {F}$ increase monotonically with each of $\displaystyle \tilde V_{rms}$ , $\displaystyle \tilde W_{rms}$ and $\displaystyle \tilde \tau _{XY}^R$ , and it can be realized from figure 13(d–i) that, even from a quantitative point of view, general trends emerge. Eventually, the following fitting relationships can be proposed (together with their accuracy levels):

(3.3) \begin{align} \mathcal {D}= 11.5 \times \left [\tilde V_{rms}\right ]^{0.6}, \quad \quad NRMS_{error}\approx 11\, \%, \end{align}

(3.4) \begin{align} \mathcal {D}= 6 \times \tilde W_{rms}, \quad \quad NRMS_{error} \approx 13\, \%, \end{align}

(3.5) \begin{align} \mathcal {D}= 12.5 \times \left [\tilde \tau _{XY}^R\right ]^{0.4}, \quad \quad NRMS_{error} \approx 11\, \%, \end{align}

(3.6) \begin{align} \mathcal {F}= 11.5 \times \tilde V_{rms}, \quad \quad NRMS_{error} \approx 18\, \%, \end{align}

(3.7) \begin{align} \mathcal {F}= 1.8 \times \left [\tilde W_{rms}\right ]^{2} + 1.6 \times \tilde W_{rms}, \quad \quad NRMS_{error}\approx 37\, \%, \end{align}

(3.8) \begin{align} \mathcal {F}= 11 \times \left [\tilde \tau _{XY}^R\right ]^{0.6}, \quad \quad NRMS_{error}\approx 29\, \%, \end{align}

where the normalized r.m.s., $NRMS_{error}$ , is evaluated by dividing the conventional $RMS_{error}$ by the mean value of either $\mathcal {D}$ or $\mathcal {F}$ . The ranges of validity of the relations proposed are

(3.9) \begin{align} 0 \leq \tilde V_{rms} \lesssim 0.25, \quad 0 \leq \tilde W_{rms} \lesssim 0.85, \quad 0 \leq \tilde \tau _{XY}^R \lesssim 0.085. \end{align}

In the remainder of this section, we aim to demonstrate (i) that transpiration strongly controls the depression in the velocity shift over a wide range of textured boundaries and (ii) that correlating $\tilde V_{rms}$ to the macroscopic coefficients of the homogenization model permits the use of (3.3) and (3.6) for an a priori estimate of the roughness-function-related parameters $\mathcal {D}$ and $\mathcal {F}$ .

Figure 14. Values of the parameter $\mathcal {D}$ plotted against the r.m.s. of the turbulent fluctuations in the wall-normal velocity at the plane $Y=0$ . In panel (a), results from the literature for channels roughened with streamwise-elongated, spanwise-elongated or three-dimensional elements are shown: blank square, Cheng & Castro (Reference Cheng and Castro2002); red circles, Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003); purple triangles, Orlandi & Leonardi (Reference Orlandi and Leonardi2006); green squares, Burattini et al. (Reference Burattini, Leonardi, Orlandi and Antonia2008); grey diamonds, Orlandi & Leonardi (Reference Orlandi and Leonardi2008); blank circles, Hao & García-Mayoral (Reference Hao and García-Mayoral2025). In panel (b), the results of Hao & García-Mayoral (Reference Hao and García-Mayoral2025) for symmetric channels bounded by either deep (red diamonds) or shallow (grey squares) porous substrates are plotted, together with the values of the present homogenization-based simulations (light-blue triangles). Solid lines refer to correlation (3.3), while the linear relationship by Orlandi & Leonardi (Reference Orlandi and Leonardi2008) is plotted with dashed lines.

Orlandi et al. (Reference Orlandi, Leonardi, Tuzi and Antonia2003) demonstrated that the principal characteristics of the flow over a rough surface are closely related to the presence of wall-normal velocity distribution at the interface between the protrusions and the overlying turbulent boundary layer. A more formal description of this dependence has been proposed by Orlandi et al. (Reference Orlandi, Leonardi and Antonia2006) and Orlandi & Leonardi (Reference Orlandi and Leonardi2006) who found good correlation between the quantity $\mathcal {D}= U_{slip}^+ - \Delta U^+$ and the r.m.s. fluctuations of the wall-normal velocity at the plane passing through the crests of the roughness elements. Later, Orlandi & Leonardi (Reference Orlandi and Leonardi2008) explored the relationship between $\mathcal {D}$ and $\tilde V_{rms}$ for walls with different textures, by collecting and plotting many results from the literature (Cheng & Castro Reference Cheng and Castro2002; Leonardi et al. Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003; Orlandi & Leonardi Reference Orlandi and Leonardi2006; Burattini et al. Reference Burattini, Leonardi, Orlandi and Antonia2008; Flores & Jiménez Reference Flores and Jiménez2006) together with new ones related to the flow over surfaces roughened with longitudinal/transverse bars or various three-dimensional patterns. They concluded their study proposing the correlation $\displaystyle \mathcal {D}=\ ({B}/{\kappa }) \, \tilde V_{rms}$ , with $B$ and $\kappa$ as by (3.1). Most of the data considered by Orlandi & Leonardi (Reference Orlandi and Leonardi2008) in addition to the recent results by Hao & García-Mayoral (Reference Hao and García-Mayoral2025) are presented in figure 14(a); the strong correlation between $\mathcal {D}$ and $\tilde V_{rms}$ is evident, and the linear relationship by Orlandi & Leonardi (Reference Orlandi and Leonardi2008), plotted with $B=5.5$ and $\kappa =0.4$ is found to perform well, where $NRMS_{error}$ is below $12\, \%$ . Interestingly, good correlation between $\mathcal {D}$ and $\tilde V_{rms}$ can also be realized in figure 14(b) for the turbulent flow over permeable boundaries, based on the values reported by Hao & García-Mayoral (Reference Hao and García-Mayoral2025) plotted next to the results of the present macroscopic DNS; the $NRMS_{error}$ for Orlandi–Leonardi relationship is approximately $23\, \%$ . With regard to the present correlation (3.3), the deviations are comparable to those reported above, with $NRMS_{error}\approx 14\, \%$ for the rough walls and $\approx 18\, \%$ for the porous boundaries, even for values of $\tilde V_{rms}$ much larger than the validity limit (3.9) of our simulations. Figure 14 thus confirms that $\tilde V_{rms}$ is a key parameter which controls the roughness function in the turbulent flow over rough/porous boundaries and that (3.3) performs well even for quite large values of $\tilde V_{rms}$ . The major difficulty in putting (3.3), or Orlandi–Leonardi correlation, to practical use is that $\tilde V_{rms}$ is not available until a full simulation of the turbulent flow above a textured wall is conducted.

The crux of the matter is thus the search of a simplified expression for $\tilde V_{rms}$ , as function of the macroscopic coefficients which permit to describe the near wall flow. After some efforts, we have found that the parameter $\Psi$ defined as

(3.10) \begin{align} \Psi = \left (\frac {\mathcal {K}_{xy}^{itf, +}}{\lambda _x^+} + \frac {\mathcal {K}_{zy}^{itf, +}}{\lambda _z^+}+ \sqrt {\mathcal {K}_{yy}^+}\right ) \, \left (\frac {\lambda _z^+}{\lambda _x^+}\right )^{0.25} \end{align}

is well correlated to $\tilde V_{rms}$ , as shown in figure 15. It is worth highlighting that (3.10) is based on the coefficients present in the boundary condition for the transpiration velocity (2.21b ): the parameters $\displaystyle \ ({\mathcal {K}_{xy}^{itf, +}}/{\lambda _x^+})$ and $\displaystyle ({\mathcal {K}_{zy}^{itf, +}}/{\lambda _z^+})$ appear when the streamwise/spanwise Navier-slip conditions $\displaystyle (\ ({\partial \hat {u}}/{\partial \hat {y}})+\ ({\partial \hat {v}}/{\partial \hat {x}}) ) |_{0}=\ ({\hat {u} |_{0}}/{\hat \lambda _x})$ and $\displaystyle ( ({\partial \hat {w}}/{\partial \hat {y}})+\ ({\partial \hat {v}}/{\partial \hat {z}} ) ) |_{0}=\ ({\hat {w} |_{0}}/{\hat \lambda _z})$ are substituted into the second and the third terms on the right-hand side of (2.21b ) and the equation is recast in wall units, while $\displaystyle \sqrt {\mathcal {K}_{yy}^+}$ quantifies the role of the intrinsic permeability for porous boundaries. The presence of $\displaystyle ({\lambda _z^+}/{\lambda _x^+})$ in (3.10) permits us to differentiate walls with spanwise-preferential slip ( $\lambda _z^+\gt \lambda _x^+$ ) from those exhibiting preferential streamwise slip ( $\lambda _x^+\gt \lambda _z^+$ ), and implies that, for the same values of $\displaystyle ({\mathcal {K}_{xy}^{itf, +}}/{\lambda _x^+})$ , $\displaystyle \ ({\mathcal {K}_{zy}^{itf, +}}/{\lambda _z^+})$ and $\displaystyle \sqrt {\mathcal {K}_{yy}^+}$ , relatively stronger transpiration is associated with the former wall patterns (e.g. substrates with transverse inclusions). Based on the data plotted in figure 15, we can propose the fitting equation

(3.11) \begin{align} \tilde V_{rms} = 0.00075 \, \Psi ^3+ 0.002 \, \Psi ^2, \end{align}

for which the $NRMS_{error}$ is less than $10\, \%$ . Substituting (3.11) into (3.3) and (3.6), we finally obtain the following expressions for the roughness-function-related quantities:

(3.12) \begin{align} \mathcal {D}= U_{slip}^+ - \Delta U^+= 11.5 \times \big (0.00075 \,\Psi ^3+ 0.002 \,\Psi ^2\big )^{0.6}, \end{align}

(3.13) \begin{align} \mathcal {F}= \Delta \lambda ^+ - \Delta U^+= 11.5 \times \big (0.00075 \,\Psi ^3+ 0.002 \,\Psi ^2\big ), \\[6pt] \nonumber \end{align}

valid up to $\Psi \approx 6$ . These expressions are plotted in figure 16 together with the results obtained from the homogenization-based DNS conducted for the porous patterns TC, LC, TM and LM. Estimates of $\mathcal {D}$ and $\mathcal {F}$ for the turbulent flow over a perturbed wall of given microstructure and given value of $\ell ^+= \epsilon Re_{\tau (\mathcal {M})}$ are thus available provided (i) $\ell ^+$ is lower than approximately 40, and (ii) outer layer similarity is maintained.

Figure 15. The r.m.s. of turbulent fluctuations in the transpiration velocity at $Y=0$ , plotted against the compound macroscopic parameter $\Psi$ for the different porous patterns considered (same symbols as in figure 13). The solid line represents a third-order polynomial fitting.

As a side remark, we observe that the relations obtained in this section have been generated by fitting data that pertain to the turbulent flow in a channel with asymmetric boundaries, and all the quantities (mean streamwise velocity, turbulence statistics and macroscopic coefficients) have been normalized with the macroscopic-pressure-gradient-based shear velocity, $u_{\tau (\mathcal {M})}$ . Since different choices appear in the literature, we provide in appendix Bkey quantities scaled with the total stress at the bottom wall. It is also shown that (3.12) and (3.13) remain reasonably accurate, independently of the choice of $u_\tau$ .

Figure 16. The roughness-function-related quantities $\mathcal {D}$ and $\mathcal {F}$ , plotted against the parameter $\Psi$ for the different porous patterns considered (same symbols as in figure 13, filled for $\mathcal {D}$ and empty for $\mathcal {F}$ ). Correlations (3.12) and (3.13) are plotted with solid lines.

3.4. Can we make a priori predictions?

It is useful to assess the accuracy of (3.12) and (3.13) for the turbulent flow over perturbed boundaries different from the porous ones based on which these correlations have been generated. In particular, we choose to check the generality of the relations above by validating them against existing numerical/experimental results for the motion over rough, impermeable walls ( $\mathcal {K}_{yy}=0$ ) with either two- or three-dimensional wall corrugations. The results in figure 17 are related to the turbulent flow $ (\displaystyle Re_{\tau (\mathcal {M})}=182.70 )$ in a symmetric channel delimited by walls roughened with in-line patterns of cubical protrusions having side length $e$ and pitch $\ell =2e$ . The Oseen-based upscaled coefficients, sensitive to the level of near-interface advection and hence to the value of $\ell ^+$ , are evaluated for $\ell ^+ = (0, 12, 23.9, 35.9, 47.8)$ and are plotted in figure 17(a). The corresponding values of $\Psi$ are $(0, 2.84, 4.39, 5.78, 7.29)$ , and are used to predict the behaviour of the quantity $\mathcal {D}$ in figure 17(b). These predictions are compared with the numerical results of the feature-resolving simulations by Hao & García-Mayoral (Reference Hao and García-Mayoral2025), and good agreement is observed. To highlight the need of incorporating near-wall advection into the homogenization model, the calculations have been repeated by setting $Re_{slip}$ equal to $0$ in (2.24) and (2.25), and significant errors in the predictions of $\mathcal {D}$ are found when $\ell ^+$ exceeds 10. In fact, the Stokes-based results coincide with the Oseen-based ones for $\ell ^+ \lesssim 10$ , a threshold similar to that reported by Ahmed & Bottaro (Reference Ahmed and Bottaro2024) for laminar channel flows. Since $\Delta U^+ = U_{slip}^+ - \mathcal {D}$ , the calculation of the roughness function based on the values obtained for $\mathcal {D}$ requires knowledge of the slip velocity, $U_{slip}^+$ . It may be tempting to use the approximation $\displaystyle U_{slip}^+= \lambda _x^+ ({\partial \overline {U}^+}/{\partial Y^+}) \big |_{Y=0} \approx \lambda _x^+$ , but care must be exerted, since significant errors appear in the predicted $U_{slip}^+$ as $\ell ^+$ becomes large, as a result of the large Reynolds stress generated at the channel virtual boundary in $Y=0$ (Hao & García-Mayoral Reference Hao and García-Mayoral2025). It is the approximation $\displaystyle ({\partial \overline {U}^+}/{\partial Y^+}) \big |_{Y=0} =1$ which eventually breaks down. In fact, based on the values of $U_{slip}^+$ and $\lambda _x^+$ reported by Hao & García-Mayoral (Reference Hao and García-Mayoral2025), we observe that the absolute deviations between the two quantities, for $\ell ^+ = (12, 23.9, 35.9, 47.8)$ , are, respectively, $(1\, \%, 4\, \%, 28\, \%, 53\, \%)$ . If we substitute $U_{slip}^+ = \lambda _x^+$ into (3.12), with the Oseen-based coefficients, we obtain $\Delta U^+ \approx (-0.73, -1.48, -2.49,-4.01)$ , progressively deviating from the values computed by Hao & García-Mayoral (Reference Hao and García-Mayoral2025): $\Delta U^+ \approx (-0.50, -1.69, -3.77, -4.91)$ .

Figure 17. Turbulent flow ( $Re_\tau \approx 180$ ) in a symmetric channel whose top/bottom boundaries are roughened with cubes (in-line arrangement) of size-to-pitch ratio $e/\ell = 0.5$ , with the spacing in wall units, $\ell ^+=\epsilon Re_{\tau (\mathcal {M})}$ , varied up to $50$ . Values of the macroscopic coefficients are plotted against $\ell ^+$ in panel (a). In panel (b), the behaviour of the parameter $\mathcal {D}$ based on (3.12) is shown (blue curve), and is validated against the results by Hao & García-Mayoral (Reference Hao and García-Mayoral2025) obtained from full simulations (squares). The black dashed curve refers to the predictions of (3.12) when $\Psi$ is evaluated with the Stokes-based upscaled coefficients, neglecting near-wall inertia; they are $\lambda _x= \lambda _z \approx 0.0653$ and $\mathcal {K}_{xy}^{itf} \approx \mathcal {K}_{zy}^{itf}= 0.0083$ .

Figure 18. Behaviour of $\Delta U^+$ with the increase in $\ell ^+$ , for the turbulent flow over surfaces with different shapes of riblets. The proposed correlations ( $\mathcal {D}$ -based (3.12), black solid lines; $\mathcal {F}$ -based (3.13), blue solid lines) are validated against relevant DNS/experimental results from the literature (red symbols). The literature results plotted are by (a–e) Wong et al. (Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024) and (f) Bechert et al. (Reference Bechert, Bruse, Hage, Van der Hoeven and Hoppe1997); the latter were reported originally in terms of $ ({\Delta C_f}/{C_{f,\, smooth}})$ and the corresponding values of $\Delta U^+$ are obtained here employing the relation $\Delta U^+= -\ ({\Delta C_f}/{C_{f,\,smooth}}) \, [(2C_{f,\,smooth})^{-0.5}+1.25 ]$ . In all panels, the thick black dashed lines represent the simple linear dependence $\Delta U^+={\lambda _x^+}-{\lambda _z^+}= (\lambda _x - \lambda _z)\, \ell ^+$ , while the grey dashed lines (wide dashes) show the predictions for $\Delta U^+$ given by Wong et al. (Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024) based on the so-called ’viscous vortex model’.

The next case examined is that of riblets. Rather than explicitly using the expression recalled earlier, $\Delta U^+ = \ell ^+_U - \ell ^+_{Turb.}$ , which is not predictive unless turbulent simulations are conducted for any shape of the riblets (or a model allowing for a priori predictions of $\ell ^+_{Turb.}$ is formulated, e.g. the ‘viscous vortex model’ by Wong et al. (Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024)), we employ (3.13). Clearly, when $\Psi$ is vanishingly small the relationship (3.13) yields the classical viscous approximation $\Delta U^+ = \Delta \lambda ^+$ (Luchini Reference Luchini and Desideri1996; Garcia-Mayoral & Jiménez Reference Garcia-Mayoral and Jiménez2011a ). Conversely, the behaviour of the roughness function can deviate significantly from this linear equation as $\ell ^+$ increases and transpiration becomes more pronounced. Different ribleted surfaces are shown in figure 18. For each geometry, the macroscopic coefficients are calculated (appendix C) and expressed in wall units by applying (3.2b ) and (3.2c ) for different spacings $\ell ^+= \epsilon Re_{\tau (\mathcal {M})}$ within the range considered ( $0 \leq \ell ^+ \leq 36$ ); the values of $\Psi$ are accordingly between $0$ and $7.4$ . The predictions in the form $\Delta U^+$ versus $\ell ^+$ , plotted with blue solid lines, are validated against the DNS results by Wong et al. (Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024) and the experimental findings by Bechert et al. (Reference Bechert, Bruse, Hage, Van der Hoeven and Hoppe1997); a reasonably good agreement can be ascertained from the figure, including the deviation from the linear dependence departing from the viscous regime, the performance degradation when the pitch distance exceeds a threshold between $15$ and $20$ , and eventually the drag increase for large riblets’ periodicity. For larger values of $\ell ^+$ , not considered in the figure, predictions of the correlation are questionable since the resulting values of $\Psi$ are significantly beyond the applicability range of the present correlation. For instance, Gatti et al. (Reference Gatti, von Deyn, Forooghi and Frohnapfel2020) studied the turbulent flow over trapezoidal riblets (similar to those in figure 18 d but with angle of $53.5^\circ$ and height equal to $0.476\,\ell$ ) and found that $\Delta U^+$ tends to become almost constant for $\ell ^+$ larger than approximately $60$ , a behaviour which cannot be captured by (3.13). Another point to be mentioned is that, for the ribleted surfaces examined in figure 18, the $\mathcal {D}$ -based predictive relationship (3.12), with $U_{slip}^+$ set to $\lambda _x^+$ , yields results for $\Delta U^+$ (plotted with orange solid lines) which are generally of lower accuracy than those obtained from (3.13). In addition, the viscous-vortex-model-based predictions of the roughness function, provided by Wong et al. (Reference Wong, Camobreco, García-Mayoral, Hutchins and Chung2024), are added to the figure; they generally exhibit reasonable accuracy up to the optimal $\ell ^+$ value for each ribleted surface (i.e. that corresponding to the maximum attainable drag reduction).

It is worth concluding this section with some notes of caution as to the applicability of our predictive correlations.

1. (i) The expressions (3.12) and (3.13) are formulated based on the results of 16 homogenization-based DNS, considering only four configurations of the porous substrate; all cases relate to transversely isotropic patterns and are characterized by a unique value of porosity, $\theta = 0.5$ . More work is certainly needed before the generality of these predictive relationships can be fully confirmed. This is also brought to attention in view of the possible deviations between the upscaling-based numerical results (employed to obtain all the trends/relationships in this paper) and the real values of the quantities of interest, attainable for example via high-fidelity feature-resolving DNS.

2. (ii) There are configurations for which the two proposed correlations may not be tenable. An example is the case of a superhydrophobic wall, with a flat and undeformable liquid–gas interface. In this case $\tilde V_{rms}=0$ and (3.12) reduces to $\displaystyle \Delta U^+= U^+_{slip}$ (incidentally, also the expression $\displaystyle \mathcal D = \ ({B}/{\kappa }) \, \tilde V_{rms}$ of Orlandi & Leonardi (Reference Orlandi and Leonardi2008) yields the same result). Since $U_{slip}^+ \approx \lambda _x^+$ the roughness function would then have a value larger than that of the conventional viscous approximation, $\Delta U^+ = \lambda _x^+ - \lambda _z^+$ (which is retrieved by the $\mathcal {F}$ -based relationship (3.13)). For the case of superhydrophobic ribs it has been shown by Luchini (Reference Luchini2015) that results for $\Delta U^+$ obtained with either no-slip/no-shear boundary conditions or with a homogenized condition collapse well with the linear, viscous approximation until $\ell ^+ \approx 30$ .

3. (iii) The good agreement between the predictions of (3.13) and the reference results for $\Delta U^+$ in figure 18 does not imply that the correlation captures, for example, the initiation of a Kelvin–Helmholtz instability past some threshold value of $\ell ^+$ . We believe that our correlations represent an improvement over linear, viscous results (dashed lines in figure 18) to predict the roughness function, but we would not want to push this as far as stating that they capture the physics at large values of $\ell ^+$ . From a mathematical perspective, (3.13) appears to reasonably quantify drag reduction up to and beyond its maximum attainable value, for each ribleted surface considered.