Turbulent drag reduction by anisotropic permeable substrates $-$ analysis and direct numerical simulations

TL;DR

This study uses direct numerical simulations to analyze how anisotropic permeable substrates can reduce turbulent skin-friction, revealing a linear regime of drag reduction and a critical permeability beyond which drag increases due to flow instabilities.

Contribution

It provides the first detailed DNS analysis confirming theoretical predictions on anisotropic permeable substrates and identifies the critical permeability threshold for drag performance.

Findings

Drag reduction up to 25% at low permeabilities.

Degradation of drag reduction occurs beyond a critical permeability.

Spanwise-coherent structures linked to Kelvin-Helmholtz instability emerge at high permeabilities.

Abstract

We explore the ability of anisotropic permeable substrates to reduce turbulent skin-friction, studying the influence that these substrates have on the overlying turbulence. For this, we perform DNSs of channel flows bounded by permeable substrates. The results confirm theoretical predictions, and the resulting drag curves are similar to those of riblets. For small permeabilities, the drag reduction is proportional to the difference between the streamwise and spanwise permeabilities. This linear regime breaks down for a critical value of the wall-normal permeability, beyond which the performance begins to degrade. We observe that the degradation is associated with the appearance of spanwise-coherent structures, attributed to a Kelvin-Helmholtz-like instability of the mean flow. This feature is common to a variety of obstructed flows, and linear stability analysis can be used to predict it. For large permeabilities, these structures become prevalent in the flow, outweighing the drag-reducing effect of slip and eventually leading to an increase of drag. For the substrate configurations considered, the largest drag reduction observed is $\approx 20-25\%$ at a friction Reynolds number $\delta^+ = 180$.

Tables

Table 1. Table 1: Characteristics of the simulations for VANS approach (BB_E80) and Brinkman’s (BB_Br). Porosity, ϵ italic-ϵ \epsilon ; viscosity ratio ν ~ / ν ~ 𝜈 𝜈 \tilde{\nu}/\nu ; friction Reynolds number, δ + = u τ ​ δ / ν superscript 𝛿 subscript 𝑢 𝜏 𝛿 𝜈 \delta^{+}=u_{\tau}\delta/\nu ; permeability, K + = K ​ u τ 2 / ν 2 superscript 𝐾 𝐾 superscript subscript 𝑢 𝜏 2 superscript 𝜈 2 K^{+}=Ku_{\tau}^{2}/\nu^{2} ; location of zero total stress, δ w subscript 𝛿 𝑤 \delta_{w} ; friction coefficient, c f = 2 ​ ( u τ / U b ) 2 subscript 𝑐 𝑓 2 superscript subscript 𝑢 𝜏 subscript 𝑈 𝑏 2 c_{f}=2(u_{\tau}/U_{b})^{2} ; friction coefficient of the corresponding smooth channel, c f 0 subscript 𝑐 subscript 𝑓 0 c_{f_{0}} ; and change of drag defined as Δ ​ D = ( c f − c f 0 ) / c f 0 Δ 𝐷 subscript 𝑐 𝑓 subscript 𝑐 subscript 𝑓 0 subscript 𝑐 subscript 𝑓 0 \Delta D=(c_{f}-c_{f_{0}})/c_{f_{0}} . Viscous units are defined with u τ subscript 𝑢 𝜏 u_{\tau} measured at the interface plane y = 0 𝑦 0 y=0 .

Cases	$ϵ$	$\tilde{\nu }/\nu$	${\delta }^{+}$	$K^{+}$	${\delta }_w/\delta$	${\delta }_i/\delta$	$c_f(\times {10}^{-3})$	$c_{f_0}(\times {10}^{-3})$	$\Delta D$
BB_E80	0.8	-	203	1.14	1.12	0.04	10.41	8.19	27.15
BB_Br	-	1.0	204	1.19	1.11	-	10.34	8.07	28.06

Table 2. Table 2: DNS parameters. K + ​ x limit-from 𝐾 𝑥 \sqrt{K\raisebox{4.52083pt}{\scriptsize$+$}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsize$x$}} , K + ​ y limit-from 𝐾 𝑦 \sqrt{K\raisebox{4.52083pt}{\scriptsize$+$}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsize$y$}} and K + ​ z limit-from 𝐾 𝑧 \sqrt{K\raisebox{4.52083pt}{\scriptsize$+$}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsize$z$}} are the streamwise, wall-normal and spanwise permeability lengths, h + superscript ℎ h^{+} is the thickness of the substrate, Δ ​ U + Δ superscript 𝑈 \Delta U^{+} is the shift of the velocity profile in the logarithmic region, and D ​ R 180 𝐷 subscript 𝑅 180 DR_{180} and D ​ R 5000 𝐷 subscript 𝑅 5000 DR_{5000} are the values of drag reduction for δ + = 180 superscript 𝛿 180 \delta^{+}=180 and δ + = 5000 superscript 𝛿 5000 \delta^{+}=5000 , respectively, obtained using expression ( 5 ). The values D ​ R 5000 𝐷 subscript 𝑅 5000 DR_{5000} have been calculated using the smooth-channel centreline velocity from Lee & Moser ( 2015 ) . The first three substrate configurations A, B and C have thickness h = 100 ​ K y ℎ 100 subscript 𝐾 𝑦 h=100\sqrt{K_{y}} and different anisotropy ratios ϕ x ​ y subscript italic-ϕ 𝑥 𝑦 \phi_{xy} . The last three substrate configurations, C ′ , C ′′ and C ′′′ , have ϕ x ​ y ≈ 11.4 subscript italic-ϕ 𝑥 𝑦 11.4 \phi_{xy}\approx 11.4 , same as substrate C, but different thickness h / K x ℎ subscript 𝐾 𝑥 h/\sqrt{K_{x}} .

	Cases	$\sqrt{K_x^{+}}$	$\sqrt{K_y^{+}}$	$\sqrt{K_z^{+}}$	$h^{+}$	$U_b/U_{b_{sm}}$	$\Delta U^{+}$	$D{R}_{180}$	$D{R}_{5000}$
Smooth		0	0	0	0	1.0	-	-	-
${\varphi }_{xy}=\sqrt{\frac{K_x}{K_y}}\approx 3.6$	A1	0.71	0.20	0.20	19.5	1.037	0.51	5.64	3.93
	A2	1.00	0.28	0.28	28.1	1.045	0.68	7.26	5.08
	A3	1.42	0.39	0.39	38.8	1.052	0.80	8.44	5.92
	A4	1.74	0.48	0.48	48.1	1.041	0.54	6.10	4.25
	A5	2.45	0.68	0.68	68.1	0.963	-0.68	-7.38	-4.99
	A6	3.61	1.00	1.00	100.2	0.819	-3.02	-42.31	-26.58
	A7	5.50	1.52	1.52	152.7	0.616	-6.59	-143.84	-76.46
	A8	10.97	3.04	3.04	304.2	0.381	-11.03	-546.15	-194.20
${\varphi }_{xy}=\sqrt{\frac{K_x}{K_y}}\approx 5.5$	B1	1.00	0.18	0.18	18.0	1.053	0.84	8.63	6.06
	B2	1.79	0.32	0.32	32.1	1.085	1.29	12.71	9.01
	B3	2.12	0.39	0.39	39.0	1.086	1.31	12.93	9.17
	B4	2.45	0.45	0.45	45.0	1.070	1.01	10.22	7.20
	B5	3.61	0.66	0.66	65.7	0.979	-0.46	-5.24	-3.56
	B6	5.48	1.00	1.00	100.0	0.792	-3.66	-56.35	-34.47
	B7	10.89	1.99	1.99	198.4	0.517	-8.66	-261.34	-120.00
${\varphi }_{xy}=\sqrt{\frac{K_x}{K_y}}\approx 11.4$	C1	1.00	0.09	0.09	9.0	1.062	0.98	9.89	6.96
	C2	1.73	0.15	0.15	14.0	1.106	1.67	16.01	11.45
	C3	2.45	0.21	0.21	22.0	1.145	2.24	20.63	14.93
	C4	3.6	0.32	0.32	32.0	1.178	2.84	25.10	18.38
	C5	4.48	0.39	0.39	39.1	1.183	2.87	25.34	18.56
	C6	5.47	0.48	0.48	47.9	1.152	2.34	21.38	15.50
	C7	10.89	0.96	0.96	95.6	0.898	-2.21	-29.35	-18.92
$\frac{h}{\sqrt{K_x}}=1.5$	C′1	2.45	0.21	0.21	3.67	1.130	2.00	18.74	13.49
	C′2	3.61	0.32	0.32	5.40	1.171	2.70	24.12	17.62
	C′3	5.49	0.48	0.48	8.23	1.156	2.40	21.87	15.88
	C′4	10.84	0.95	0.95	16.26	0.962	-0.90	-10.84	-7.27
$\frac{h}{\sqrt{K_x}}=1.0$	C′′1	3.61	0.32	0.32	3.61	1.154	2.42	22.02	15.99
	C′′2	5.48	0.48	0.48	5.51	1.163	2.53	22.86	16.64
	C′′3	7.01	0.62	0.62	7.01	1.127	1.90	17.93	12.88
	C′′4	9.03	0.79	0.79	9.03	1.066	0.84	8.62	6.05
	C′′5	10.85	0.95	0.95	11.03	1.001	-0.12	-1.32	-0.91
$\frac{h}{\sqrt{K_x}}=0.5$	C′′′1	2.45	0.21	0.21	1.22	1.063	0.93	9.46	6.65
	C′′′2	3.62	0.32	0.32	1.86	1.091	1.36	13.35	9.48
	C′′′3	5.47	0.48	0.48	2.74	1.133	2.04	19.11	13.77
	C′′′4	7.01	0.62	0.62	3.50	1.153	2.39	21.81	15.83
	C′′′5	9.03	0.79	0.79	4.52	1.129	1.95	18.34	13.19
	C′′′6	10.83	0.95	0.95	5.42	1.092	1.30	12.88	9.13

Equations

$$\nabla p=-\nu {\mathbf{K}}^{-1}\mathbf{u}+\tilde{\nu }{\nabla }^2\mathbf{u}.$$

$$\nu {\left[\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right]}_{y=0^{+}}$$

$$\nu {\left[\frac{\partial w}{\partial y}+\frac{\partial v}{\partial z}\right]}_{y=0^{+}}$$

$${\left[-p+2\nu \frac{\partial v}{\partial y}\right]}_{y=0^{+}}$$

$${\hat{u}|}_{y=0^{+}}$$

$${\hat{w}|}_{y=0^{+}}$$

$${\hat{v}|}_{y=0^{+}}$$

$$c_f=2\frac{{\tau }_w}{U_{\delta }^2}=2\frac{1}{U_{\delta }^{+2}},$$

$$DR=-\frac{c_f-c_{f0}}{c_{f0}}=1-\frac{1}{{\left(1+\Delta U^{+}/U_{\delta 0}^{+}\right)}^2}.$$

$$\Delta U^{+}\approx {\ell }_U^{+}-{\ell }_T^{+},$$

$${u^{+}|}_{y^{+}=0}={\ell }_x^{+}{{\frac{\partial u}{\partial y}}^{+}|}_{y^{+}=0},$$

$${w^{+}|}_{y^{+}=0}={\ell }_z^{+}{{\frac{\partial w}{\partial y}}^{+}|}_{y^{+}=0},$$

$${\ell }_x^{+}=\xi \sqrt{K_x^{+}}\tanh \left(\frac{h^{+}}{\sqrt{K_x^{+}}}\right),$$

$${\ell }_z^{+}=\xi \sqrt{K_z^{+}}\tanh \left(\frac{h^{+}}{\sqrt{K_z^{+}}}\right),$$

$$\Delta U^{+}\approx \xi \left(\sqrt{K_x^{+}}-\sqrt{K_z^{+}}\right).$$

$$K_{Br}^{+}=K_y^{+}\tanh \left(\frac{\sqrt{2K_x^{+}}}{y_c^{+}}\right){\tanh }^2\left(\frac{h^{+}}{\sqrt{12K_y^{+}}}\right).$$

$$K_{Br}^{+}\approx K_y^{+}.$$

$$\Delta U^{+}=\left({\varphi }_{xy}-1\right)\sqrt{K_y^{+}}.$$

$${\tau }_w=-P_x\left(\delta +h\right),$$

$$c_f=2\frac{{\tau }_w}{U_b^2}=2\frac{{\tau }_w}{{\left(U_{b0}+\Delta U_b\right)}^2}=c_{f0}\frac{1+h/\delta }{{\left(1+\Delta U_b/U_{b0}\right)}^2}=c_{f0}\frac{{\left(1+h/\delta \right)}^3}{{\left(1+\Delta q/q_0\right)}^2},$$

$$\Delta q\approx 2\delta \left(U_{Br}+U_{Darcy}\right)+2q_{substrate},$$

$$c_f\approx c_{f0}\frac{1+3h/\delta }{1+\frac{2\sqrt{K_x}}{U_{b0}}{\frac{dU}{dy}|}_{y=0}}.$$

$$\text{bnabla}\text{bcdot}\mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\text{bcdot}\text{bnabla}\mathbf{u}=-\nabla p+\frac{1}{Re}{\nabla }^2\mathbf{u},$$

$$\begin{split}\left[I-\Delta t\frac{\beta_{k}}{\mathrm{Re}\mathrm{L}}\right]\boldsymbol{u}^{k}=\boldsymbol{u}^{k-1}+\Delta t\Big{[}&\frac{\alpha_{k}}{\mathrm{Re}}\mathrm{L}(\boldsymbol{u}^{k-1})-\gamma_{k}\mathrm{N}(\boldsymbol{u}^{k-1})\\ &-\zeta_{k}\mathrm{N}(\boldsymbol{u}^{k-2})-\left(\alpha_{k}+\beta_{k}\right)\mathrm{G}(p^{k})\Big{]}\end{split}$$

$$\left[\begin{array}[]{cc}\mathrm{A}&\mathrm{G}\\ \displaystyle\mathrm{D}&0\\ \end{array}\right]\left(\begin{array}[]{c}\mathbf{u}^{k}\\ \displaystyle p^{k}\\ \end{array}\right)=\left(\begin{array}[]{c}r^{k-1}\\ \displaystyle 0\\ \end{array}\right),$$

$$\left[\begin{array}[]{cc}\mathrm{A}&0\\ \displaystyle\mathrm{D}&-\mathrm{D}\mathrm{A^{-1}}\mathrm{G}\\ \end{array}\right]\left[\begin{array}[]{cc}\mathrm{I}&\mathrm{A^{-1}}\mathrm{G}\\ \displaystyle 0&\mathrm{I}\\ \end{array}\right]\left(\begin{array}[]{c}\mathbf{u}^{k}\\ \displaystyle p^{k}\\ \end{array}\right)=\left(\begin{array}[]{c}r^{k-1}\\ \displaystyle 0\\ \end{array}\right)$$

$$\mathrm{A}{\mathbf{u}}^{\ast }=r^{k-1},$$

$$\mathrm{D}{\mathrm{A}}^{-1}\mathrm{G}{p}^{n+1}=\mathrm{D}{\mathbf{u}}^{\ast },$$

$${\mathbf{u}}^k={\mathbf{u}}^{\ast }-{\mathrm{A}}^{-1}\mathrm{G}{p}^k,$$

$$\frac{u_c/t_c}{\nu u_c/K}\sim \frac{K}{\nu t_c}=\frac{K^{+}}{t_c^{+}}.$$

$$\Delta U_{max}^{+}\approx 0.8\times 0.38\times \left({\varphi }_{xy}-1\right).$$

$$-{\overline{u^{\prime }{v}^{\prime }}}^{+}+\frac{d{U}^{+}}{d{y}^{+}}=\frac{{\delta }^{\prime +}-y^{\prime +}}{{\delta }^{\prime +}},$$

$${\int }_{{\ell }_{T,P}^{+}}^{H^{+}}-{\overline{u^{\prime }{v}^{\prime }}}^{+}d{y}^{\prime +}+U^{+}(H^{+})-U^{+}({\ell }_{T,P}^{+})=H^{+}-{\ell }_{T,P}^{+}-\frac{H^{+2}-{\ell }_{T,P}^{+2}}{2{\delta }^{\prime +}}.$$

$$\begin{array}{rl}\Delta U^{+}& =U_P^{+}(H^{+})-U_S^{+}(H^{+})\\ & =\underset{\begin{array}{c}{U}_{slip}^{+}\end{array}}{\underbrace{U_P^{+}({\ell }_{T,P}^{+})}}-U_S^{+}({\ell }_{T,P}^{+})\underset{\begin{array}{c}{\mathcal{T}}_{uv}\end{array}}{\underbrace{-{\int }_{{\ell }_{T,P}^{+}}^{H^{+}}\left[\left(-{\overline{u^{\prime }{v}^{\prime }}}_P^{+}\right)-\left(-{\overline{u^{\prime }{v}^{\prime }}}_S^{+}\right)\right]d{y}^{\prime +},}}\end{array}$$

Full text

1 Introduction

The high skin friction experienced in turbulent flows represents a problem for several engineering applications, such as pipelines and transportation vehicles. The need is therefore to develop new technologies that reduce turbulent drag, preferably passive, since in contrast with active technologies, these do not require an energy input and have generally lower manufacturing costs. In this paper we present the potential of anisotropic permeable substrates, a passive technology, to reduce turbulent skin friction, as has recently been proposed by Abderrahaman-Elena & García-Mayoral (2017).

Most of the literature in turbulent flows over permeable substrates has focused on isotropic materials, observing a substantial increase in drag with respect to a smooth wall (Breugem et al., 2006; Rosti et al., 2015; Kuwata & Suga, 2016). This increase has often been attributed to the onset of large spanwise-coherent structures, which increase the momentum transfer and thus the Reynolds stresses near the wall. Here, we study the effect of anisotropy and provide physical insight into the behaviour of anisotropic permeable substrates in turbulent flows for drag-reducing purposes, when the permeability is preferential in the streamwise direction. Recent studies have also covered anisotropic substrates, albeit not considering the case of streamwise-preferential permeability (Kuwata & Suga, 2017; Suga et al., 2018).

Previous studies have shown that streamwise-preferential complex surfaces can reduce drag in turbulent flows (Bechert et al., 1997; Luchini et al., 1991; Jiménez, 1994; Gómez-de-Segura et al., 2018b). This is indeed the case for some of the most common passive technologies for drag reduction, such as riblets or superhydrophobic surfaces. Recently, Abderrahaman-Elena & García-Mayoral (2017) suggested that the drag reduction ability of anisotropic permeable substrates is based on the same mechanism. The general idea is that complex surfaces can reduce drag if they offer more resistance to the cross flow than to the streamwise mean flow. When the surface texture is vanishingly small compared to the near-wall turbulent structures, the effect of complex surfaces can be reduced to an apparent slip in the tangential directions. Luchini et al. (1991), Luchini (1996) and Jiménez (1994) showed that the change in drag is proportional to the difference between the streamwise and spanwise slips. Hahn et al. (2002) observed this behaviour also in turbulent flows over substrates permeable in the streamwise and spanwise directions only. They observed that the streamwise slip is beneficial for drag reduction, while the spanwise slip has an opposite effect. Their substrates, however, were ideal, in the sense that they were impermeable in the wall-normal direction. Hence, the work by Hahn et al. (2002) is closely connected to studies where only tangential slips are allowed, while the surface remains impermeable, such as those carried out by Min & Kim (2004) or Busse & Sandham (2012) in the context of superhydrophobic surfaces. Recently, Rosti et al. (2015) have studied permeable substrates with very low wall-normal permeability, which would also fall under this category. The analysis by Gómez-de-Segura et al. (2018a) shows that the deleterious effect of the spanwise slip saturates if this is not accompanied by a corresponding wall-normal transpiration. Therefore, surfaces with isotropic slip can also reduce drag, although suboptimally.

The linear theory of Luchini et al. (1991) and Jiménez (1994) is valid only as long as the texture lengthscales are small compared to the characteristic lengthscales of near-wall turbulence. As the texture size increases, additional deleterious effects set in, breaking down the drag-reducing performance and eventually leading to an increase of drag. The mechanisms behind these deleterious effects vary from one technology to another. In riblets, for instance, the degradation of performance is due to the appearance of spanwise-coherent rollers, which arise from a Kelvin-Helmholtz instability (García-Mayoral & Jiménez, 2011). These structures are in fact a common feature to a variety of obstructed flows (Ghisalberti, 2009).

Several studies on permeable substrates have also reported the existence of such structures (Breugem et al., 2006; Kuwata & Suga, 2016; Zampogna & Bottaro, 2016; Suga et al., 2017). In these studies, the large increase of the Reynolds stresses compared to that over a smooth wall and the subsequent increase in drag was associated to the presence of Kelvin-Helmholtz rollers. Abderrahaman-Elena & García-Mayoral (2017) suggested the formation of these rollers as a possible drag-degrading mechanism for anisotropic permeable substrates. They proposed a model to bound the maximum achievable drag reduction based on the onset of the Kelvin-Helmholtz-like instability. Gómez-de-Segura et al. (2018b) extended the analysis and identified the wall-normal permeability as the governing parameter in this instability. This result agrees with the work performed by Jiménez et al. (2001), who observed the formation of Kelvin-Helmholtz rollers over substrates which were permeable in the wall-normal direction only, and inferred that the relaxation of the impermeability condition at the wall was sufficient to elicit the rollers.

Several drag-reducing surfaces show a linear regime, where the drag reduction increases linearly with a certain characteristic length of the texture, followed by a saturation and an eventual increase of drag (García-Mayoral & Jiménez, 2011). Although the same has not been shown for anisotropic permeable substrates, the similarities between the drag reduction curves of riblets and those of seal fur by Itoh et al. (2006) suggest a similar behaviour (Abderrahaman-Elena & García-Mayoral, 2017). The effect of the seal fur studied by Itoh et al. (2006) would be to some extend that of an anisotropic permeable material, since it is a layer of hairs preferentially aligned in the streamwise direction.

In the current work, we investigate the drag reduction ability of anisotropic permeable substrates. The aim of this work is to understand how the overlying turbulent flow is modified by the presence of such substrates and build predictive models to estimate their drag-reducing behaviour. For that, we perform a series of DNSs of channel flows bounded by permeable substrates, which are selected using the information obtained from a linear stability theory and the linearised theory of Luchini et al. (1991) and Jiménez (1994) for drag reduction.

The present paper is organised as follows. In §2 we discuss several models to characterise the flow within the permeable substrates and present the analytic solution to the model subsequently used, Brinkman’s model. How streamwise-preferential permeable substrates can reduce drag is explained in §3, where we also discuss the theoretical models derived by Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b). The former provides estimates for the expected drag reduction in the linear regime, while the latter bounds the achievable drag reduction based on linear stability theory. These models allow us to select particular permeable substrates for the subsequent DNS study. Details for the DNS setup are presented in §4. In §5, we present the DNS results for the permeable substrates selected and assess the validity of the theoretical models. Drag reduction curves for different anisotropic permeable substrates are also included, allowing to define design guidelines for optimal substrate configurations. Finally, conclusions are summarised in §6.

2 Flow within the permeable substrate

Following Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b), we focus on permeable materials where the pores are much smaller than any near-wall turbulent lengthscale. We therefore opt for a macroscopic, homogenised approach to model the flow within the permeable medium, due to the high resolution required otherwise to explicitly solve the flow within the pores. The permeable medium is modelled as homogeneous, by defining a local, instantaneous average solution of the flow within the fluid-solid matrix.

A classical approach to characterise the homogenised flow within a permeable medium is Darcy’s equation (Darcy, 1856). This is the simplest model amongst the continuum approaches, and results from a volume average of the Stokes equation over many pores/particulate obstacles. Note that under the assumption of vanishingly small pore size, such averages could still be conducted in small volumes compared to the scales of the overlying flow. Darcy’s equation is a balance between the pressure gradient across the permeable medium and the viscous drag caused by the pressure of the solid matrix. More sophisticated continuum approaches used in the literature include homogenisation techniques (Zampogna & Bottaro, 2016; Lācis & Bagheri, 2017) or the Volume Averaged Navier-Stokes equations (VANS) (Whitaker, 1996; Ochoa-Tapia & Whitaker, 1995b, a). Several authors have recently used the latter to study flows over permeable substrates (Breugem et al., 2006; Tilton & Cortelezzi, 2008; Rosti et al., 2015).

The volume average, implicit in Darcy’s equation, accounts for the viscous stresses caused by velocity gradients over lengths smaller than the averaging one. This effectively filters out diffusive effects acting over larger lengthscales. If the latter are relevant, they can be accounted for by including a macroscopic diffusive term, yielding Brinkman’s equation (Brinkman, 1947),

[TABLE]

The first two terms in equation (1) constitute Darcy’s equation, and the last term, $\tilde{\nu}\nabla^{2}\mathbf{u}$, is the Brinkman term, with $\mathbf{u}$ the velocity vector, p the kinematic pressure, and $\nu$ and $\tilde{\nu}$ the molecular viscosity of the fluid and the effective macrocopic viscosity, respectively. The homogenised flow within the permeable substrate and the different lengthscales accounted for by the various terms in equation (1) are illustrated in figure 1. Panel (c) portrays the flow between the obstacles, which results in Darcy’s equation when averaged, while panel (b) portrays the large scale diffusion missed by the volume averaging and captured by the Brinkman term. Brinkman’s model is suitable for substrates made up of open matrices of obstacles, where fluid regions are significantly interconnected and diffusion can act efficiently over large scales. But, it does not represent correctly substrates made up of microducts essentially isolated from each other, where diffusion cannot act over scales larger than the pores (Lévy, 1983; Auriault, 2009). Abderrahaman-Elena & García-Mayoral (2017) used this distinction to characterise substrates as ‘highly-connected’ or ‘poorly-connected’, and argued that the former offered better properties for drag reduction. The two types of materials are illustrated in figure 2.

In poorly-connected substrates, Darcy’s equation provides a reasonable model for the flow within (Lévy, 1983; Auriault, 2009), but it cannot capture the interfacial layer that forms immediately below the substrate-fluid interface, where the velocity transitions from Darcy’s velocity deep inside the substrate to a certain slip velocity at the interface plane. For that, the ‘jump condition’ proposed by Beavers & Joseph (1967) is generally used, which imposes a slip velocity proportional to the external shear at the substrate-fluid interface. The constant of proportionality, $\alpha_{BJ}$, accounts for the structure of the permeable material and is determined empirically.

In highly-connected substrates, in contrast, the Brinkman model allows to capture the interfacial region, under certain assumptions. This equation is also a volume averaging model, so it implicitly assumes that any small volume within the substrate contains a large number of obstacles. However, as the averaging volume approaches the interface with the free flow, this assumption would eventually cease to hold. The specialised literature shows no general agreement regarding the treatment of the substrate-fluid interface (Lācis & Bagheri, 2017; Zampogna & Bottaro, 2016; Ochoa-Tapia & Whitaker, 1995a; Le Bars & Worster, 2006). Some studies impose jump conditions, as discussed previously, although these can be of different types, such as a jump in velocity (Beavers & Joseph, 1967), a jump in shear stress but not in velocity (Ochoa-Tapia & Whitaker, 1995a), or continuity of both velocity and shear stress (Vafai & Kim, 1990; Le Bars & Worster, 2006; Battiato, 2012, 2014). Previous studies have shown an analogy between Brinkman’s model and Beavers and Joseph’s ‘jump condition’ at the substrate-fluid interface (Taylor, 1971; Neale & Nader, 1974; Abderrahaman-Elena & García-Mayoral, 2017). Other studies, in contrast, define an adaptation region of certain thickness where the permeability transitions smoothly from its value within the substrate to infinity in the free flow. This is the case of Breugem et al. (2006), where they use the more general VANS approach with an adaptation region of thickness $\delta_{i}$. For substrates where the inertial terms are negligible, this approach would be analogous to using Brinkman’s model and ‘blurring’ the solution with a moving average of thickness $\delta_{i}$.

The analysis of Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b) suggested that highly-connected materials would yield greater drag reduction. Furthermore, for the small values of permeabilities considered in this study, the flow within the substrate would be dominantly viscous. In this scenario, the Brinkman model provides a simple but reasonable approximation. We therefore follow the above works and use Brinkman’s equation to model the flow within the substrate. For simplicity, we assume that pores are infinitely small, so the continuum hypothesis would hold for any vanishingly small volume, and Brinkman’s equation remains valid near the interface (Vafai & Kim, 1990). For larger permeabilities, the inertial terms might also become important and they need to be considered by including an additional Forchheimer term (Forchheimer, 1901; Joseph et al., 1982; Whitaker, 1996).

2.1 Analytic solution of Brinkman’s equation

In the present work we consider channels of height $2\delta$ delimited by two identical anisotropic permeable substrates of thickness $h$, as sketched in figure 1. The substrate-channel interfaces are located at $y=0$ and $y=2\delta$, and the substrates are bounded by impermeable walls at $y=-h$ and $y=2\delta+h$. Throughout the paper we will refer to the free-flow region between $y=0$ and $y=2\delta$ as ‘channel’ and to the permeable region below $y=0$ (or above $y=2\delta$) as ‘substrate’. The flow within the permeable substrates is modelled using equation (1), where the fluid density $\rho$ is assumed to be unity for convenience. The simplicity of Brinkman’s equation allows to solve it analytically, and the particularised solution at the substrate-channel interface can be implemented as boundary condition for the DNS of the channel, fully coupling the flow in both regions. The procedure to solve Brinkman’s equation is detailed in Appendix A. Here only the problem formulation and its solution are presented.

As discussed above, poorly-connected substrates have negligible macroscale viscous effects, which in equation (1) can be interpreted as having $\tilde{\nu}=0$, recovering Darcy’s equation. Highly-connected media, in turn, would asymptotically tend to have macroscale diffusion as efficient as a free flow, so $\tilde{\nu}\approx\nu$ (Tam, 1969; Lévy, 1983; Neale & Nader, 1974; Abderrahaman-Elena & García-Mayoral, 2017). Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b) suggested that such materials would have a better potential for drag reduction, as it will be discussed in §3. Here we follow them and assume $\tilde{\nu}=\nu$. The permeable substrates are characterised then by their thickness, $h$, and their permeabilities $K_{x}$, $K_{y}$ and $K_{z}$ in the streamwise, $x$, wall-normal, $y$, and spanwise, $z$, directions, respectively, which are considered to be the principal directions of the permeability tensor $\mathbf{K}$ in equation (1). The tensor has dimensions of length squared, and is a measure of the ability of the fluid to flow through a permeable medium. When $\mathbf{K}\rightarrow\infty$ the medium offers no resistance to the flow, and when $\mathbf{K}=0$ an impermeable medium is recovered.

Let us consider the lower substrate between $y=-h$ and $y=0$. To solve equation (1), we impose no slip and impermeability at $y=-h$, and continuity of the tangential and normal stresses at the substrate-channel interface, i.e. at $y=0$. The solution within the substrate is coupled to the flow within the channel by imposing the continuity of the three velocity components. The resulting boundary conditions at $y=0$ are then

[TABLE]

where $y=0^{+}$ and $y=0^{-}$ correspond to the channel and the substrate sides of the interface, respectively. Under the above assumptions, the boundary conditions (2) can be further simplified. The continuity of tangential stresses becomes that of $\partial u/\partial y$ and $\partial w/\partial y$, and the continuity of normal stresses that of $p$. Equation (1) is then solved by taking Fourier transforms in the tangential directions $(x,z)$. Following the derivations presented in Appendix A, the analytic solution particularised at $y=0$ provides the following expressions for the velocities,

[TABLE]

where the hat denotes variables in Fourier space. The coefficients $\mathcal{C}_{ij}$ are complex and depend on the structure of the permeable substrate through $K_{x}$, $K_{y}$, $K_{z}$ and $h$, as well as on the overlying flow through the streamwise and spanwise wavenumbers, $\alpha_{x}$ and $\alpha_{z}$, or the corresponding wavelengths, $\lambda_{x}=2\pi/\alpha_{x}$ and $\lambda_{z}=2\pi/\alpha_{z}$. The same procedure can be used to obtain a symmetric solution for the upper substrate, and the resulting expressions for the interface at $y=2\delta$ can be found in Appendix A. The effect of the permeable substrates on the channel flow is introduced through equations (3) and the corresponding equations at $y=2\delta$, which serve as boundary conditions.

To illustrate how the coefficients in equation (3) vary with the wavelengths, figure 3 shows maps of $\mathcal{C}_{uu}^{+}$, $\mathcal{C}_{ww}^{+}$ and $\mathcal{C}_{vp}^{+}$, which have zero imaginary part, as a function of $\lambda_{x}^{+}$ and $\lambda_{z}^{+}$ for a particular substrate. The superscript ‘$+$’ denotes viscous units, where magnitudes are normalised using the kinematic viscosity, $\nu$, and the friction velocity at the substrate-channel interface, $u_{\tau}=\sqrt{\tau_{w}}$. Note that the total stress at that location, $\tau_{w}$, accounts for both the viscous and the Reynolds stresses. $\mathcal{C}_{uu}^{+}$ and $\mathcal{C}_{ww}^{+}$ relate the streamwise and spanwise velocities with their corresponding wall-normal gradients, respectively, and are connected to the slip boundary conditions typically used in slip-only simulations (Hahn et al., 2002; Min & Kim, 2004; Busse & Sandham, 2012). $\mathcal{C}_{vp}^{+}$ represents an impedance relating the wall-normal velocity and the pressure (Jiménez et al., 2001). The slip coefficients $\mathcal{C}_{uu}^{+}$ and $\mathcal{C}_{ww}^{+}$ are purely real, so the tangential velocity is in phase with the tangential shear. The transpiration coefficient $\mathcal{C}_{vp}^{+}$ is also real but negative, so the wall-normal velocity is in anti-phase with the pressure. For the mean flow, i.e $\alpha_{x}^{+}=0$ and $\alpha_{z}^{+}=0$ (or alternatively $\lambda_{x}^{+}\rightarrow\infty$ and $\lambda_{z}^{+}\rightarrow\infty$), out of the 9 coefficients from equation (3) only $\mathcal{C}_{uu}^{+}$ and $\mathcal{C}_{ww}^{+}$ are non-zero and their value decreases as the wavenumbers increase, as shown in figure 3. In contrast, the transpiration coefficient $\mathcal{C}_{vp}^{+}$ is zero for the mean flow and becomes increasingly negative as the wavenumbers increase, since short wavelengths penetrate more easily through the substrate. In Gómez-de-Segura et al. (2018b), we conducted preliminary DNSs of channel flows with permeable substrates where only these three coefficients from equation (3) were included. The DNSs presented here in §5 show that the other coefficients modulate the results, and this modulation can become significant as the permeability increases.

3 Theoretical models

In this section, we present the theoretical models introduced by Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b) to estimate the drag reduction that permeable substrates can achieve. We also discuss the effect on internal and external flows and how they relate.

3.1 Drag reduction from surface manipulations

The friction coefficient, $c_{f}$, can be defined as

[TABLE]

where the density is assumed to be unity. The choice on the reference velocity $U_{\delta}$ depends on the type of flow studied. In external flows, the free stream velocity is typically used, while in internal flows the bulk velocity is more common. The substrates studied here would mainly be aimed at external flow applications, for instance as coatings in vehicle surfaces. The simulations, however, have been conducted in channels for simplicity. In this framework, García-Mayoral & Jiménez (2011) argued that choosing the centreline velocity as the reference for $c_{f}$ permitted a closer comparison with external-flow friction coefficients.

In the case of small surface textures, their effect is confined to the near-wall region. According to the classical theory of wall turbulence, sufficiently far away from the wall, the only effect of any surface manipulation is to modify the intercept of the logarithmic law, while the Kármán constant and the wake function remain unaltered (Clauser, 1956). The centreline velocity is then $U_{\delta}^{+}=U_{\delta 0}^{+}+\Delta U^{+}$, where the subscript ‘[math]’ indicates values for a reference smooth channel and $\Delta U^{+}$ is the shift of the logarithmic velocity profile with respect to the smooth wall. The drag reduction ($DR$) can then be expressed in terms of $\Delta U^{+}$,

[TABLE]

If $\Delta U^{+}>0$, the logarithmic region is shifted upwards and drag is reduced. Conversely, if $\Delta U^{+}<0$, the logarithmic region is shifted downwards and drag is increased. Note that $DR$ depends on the friction Reynolds number, $\delta^{+}$, through $U_{\delta 0}^{+}$, while $\Delta U^{+}$ does not. The latter therefore provides a more universal measure, as it can be extrapolated to higher $\delta^{+}$ (García-Mayoral & Jiménez, 2011; Spalart & McLean, 2011; Gatti & Quadrio, 2016; García-Mayoral et al., 2019). The change of $DR$ with $\Delta U^{+}$ given by equation (5) and its dependence with $\delta^{+}$ are depicted in figure 4. This figure shows a decrease of $DR$ with the Reynolds number, due to larger values of $U_{\delta 0}^{+}$. This can be expected to lead to discrepancies in $DR$ between simulations and experiments at low Reynolds numbers, and industrial applications at high Reynolds numbers. To circumvent this, in the present paper we quantify drag reduction in terms of $\Delta U^{+}$.

3.2 Drag reduction from virtual origins

Drag reduction from non-smooth, passive surfaces has recently been reviewed in García-Mayoral et al. (2019) as a virtual-origin effect, where the reduction of drag is essentially caused by an offset between the positions of the virtual, equivalent smooth walls perceived by the mean flow and the overlying turbulent flow. For vanishingly small surface textures, Luchini et al. (1991) proposed that $\Delta U^{+}$ produced by any complex surface is given by

[TABLE]

where $\ell_{U}^{+}$ refers to the virtual origin experienced by the mean flow, defined as the depth below a reference plane where the mean flow would perceive a non-slipping wall; and $\ell_{T}^{+}$ refers to the virtual origin experienced by turbulence. Luchini (1996) suggested that the latter could be identified as the origin experienced by the quasi-streamwise vortices. These virtual origins are measured from a reference plane, often taken at the top plane of the surface geometry, for instance at the riblet tips (Luchini et al., 1991) or at the substrate-fluid interface plane for permeable substrates (Abderrahaman-Elena & García-Mayoral, 2017). This is where we set $y=0$. The virtual origins perceived by the mean flow and the vortices are therefore at $y^{+}=-\ell_{U}^{+}$ and $y^{+}=-\ell_{T}^{+}$, respectively. As discussed below, these are directly connected to the concepts of ‘slip lengths’ and ‘protrusion heights’ typically used in the literature.

If the surface texture is small, the overlying flow does not perceive the detail of the texture, but a homogenised effect, quantified by these virtual origins. Complex surfaces therefore change drag by causing a relative y-displacement of turbulence with respect to the mean flow, as depicted in figure 5, but turbulence remains otherwise smooth wall-like (Luchini et al., 1991; Jiménez, 1994; Luchini, 1996; García-Mayoral et al., 2019). If $\ell_{T}^{+}<\ell_{U}^{+}$, quasi-streamwise vortices are, compared to a smooth wall, shifted farther away from the origin of the mean flow. As a result, the local momentum flux close to the surface decreases, thereby reducing the shear and the skin friction. Conversely, if $\ell_{U}^{+}<\ell_{T}^{+}$, the vortices perceive a deeper origin than the mean flow and friction drag increases.

Luchini et al. (1991) and Luchini (1996) proposed that the virtual origin of turbulence is given by that of the spanwise velocity. Given that for small surface textures the velocity profile near the surface is linear, the concept of virtual origins can be represented by Robin boundary conditions at the reference plane $y^{+}=0$,

[TABLE]

where the Robin coefficients $\ell_{x}^{+}$ and $\ell_{z}^{+}$, typically referred to as the streamwise and spanwise slip lengths, are roughly equal to the depths of the virtual origins, i.e. $\ell_{U}^{+}\approx\ell_{x}^{+}$ and $\ell_{T}^{+}\approx\ell_{z}^{+}$. In addition, the mean streamwise shear is $dU^{+}/dy^{+}|_{y^{+}=0}\approx 1$ and the slip length, $\ell_{x}^{+}$, is interchangeable with the slip velocity, $U_{slip}^{+}$.

Boundary conditions of the form of equation (7) are generally used in slip-only simulations, such as in Min & Kim (2004) or Busse & Sandham (2012). In these simulations, however, the effect of $\ell_{z}^{+}$ on $\Delta U^{+}$ saturates, that is the linear expression $\Delta U^{+}\approx\ell_{x}^{+}-\ell_{z}^{+}$ is valid only for $\ell_{z}^{+}\lesssim 1$ and increasing $\ell_{z}^{+}$ beyond $\approx 4$ has only a negligible effect on $\Delta U^{+}$. Gómez-de-Segura et al. (2018a) noted that this saturation effect is a result of the impermeability condition imposed at the interface, $v=0$, as the imposed impermeability impedes the displacement of the quasi-streamwise vortices further towards the interface. This effect would be present in the drag-reducing simulations of Hahn et al. (2002) and Rosti et al. (2018), which considered zero or very low values of wall-normal permeabilities, so that $v\approx 0$ at the interface. This would not be the case for the permeable substrates in general, or for those studied in this paper in particular. For the DNSs presented in §5, we consider equal wall-normal and spanwise permeabilities, $K_{y}^{+}=K_{z}^{+}$. The slip in the spanwise direction is then always accompanied by a corresponding wall-normal transpiration, and the virtual origin perceived by turbulence is roughly given by $\ell_{T}^{+}\approx\ell_{z}^{+}$ with no saturation. For a more general case where $K_{z}^{+}\neq K_{y}^{+}$, however, the virtual origin of turbulence would deviate from $\ell_{z}^{+}$ (Gómez-de-Segura et al., 2018a).

3.3 Virtual origins for anisotropic permeable substrates

Abderrahaman-Elena & García-Mayoral (2017) derived the streamwise and spanwise slip lengths, as well as $\Delta U^{+}$, for a permeable substrate. The authors calculated $\ell_{x}^{+}$ and $\ell_{z}^{+}$ by solving the flow within the permeable medium in response to an overlying shear, obtaining a solution of the form of equation (7), a procedure that has also been followed for riblets or superhydrophobic textures (Luchini et al., 1991; Ybert et al., 2007). Abderrahaman-Elena & García-Mayoral (2017) solved equation (1) for $u$ and $w$ under homogeneous shear, for which the pressure terms zero out. This is actually the solution for mode zero, i.e. $\alpha_{x}=0$ and $\alpha_{z}=0$, in Appendix A. Obtaining the relationships between the velocities and their corresponding shears at the interface, $\ell_{x}^{+}$ and $\ell_{z}^{+}$ are

[TABLE]

where $\xi$ is the ratio between the molecular and effective viscosities of the permeable substrate, and would be $\xi\approx 1$ for highly-connected substrates with $\tilde{\nu}\approx\nu$. Note that $\ell_{x}^{+}$ and $\ell_{z}^{+}$ in equation (8) are the coefficients $\mathcal{C}_{uu}^{+}$ and $\mathcal{C}_{ww}^{+}$ for mode $(0,0)$ in equation (3). For poorly-connected substrates, Abderrahaman-Elena & García-Mayoral (2017) obtained the same solution using Darcy’s equation with Beavers & Joseph’s jump conditions at the interface, in which case $\xi$ would be the inverse of Beavers & Joseph’s constant of proportionality, $\alpha_{BJ}$.

Abderrahaman-Elena & García-Mayoral (2017) concluded that the highest performance for a given anisotropic material would be achieved for sufficiently deep substrates, where h^{+}\gtrsim\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}},\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}. In this case, both hyperbolic tangents in equation (8) tend to unity and the slip lengths become \ell_{x}^{+}\approx\xi\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} and \ell_{z}^{+}\approx\xi\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}. Introducing these results into equation (6), $\Delta U^{+}$ becomes

[TABLE]

The microstructure of the substrate, represented by $\xi$, has therefore an important effect on the drag-reducing performance of the substrate. The optimum configuration would be obtained for highly-connected materials with $\xi\approx 1$ (i.e. $\tilde{\nu}\approx\nu$), which supports our previous assumption in §2.1. Furthermore, to maximise drag reduction, we seek highly anisotropic materials, maximising the streamwise permeability, $K_{x}^{+}$, while minimising the spanwise one, $K_{z}^{+}$. Note that equation (9) considers deep substrates, so that the flow near the substrate-channel interface does not perceive the bottom no-slipping wall. This assumption eliminates the substrate thickness, $h^{+}$, from the parameter space under consideration.

The linear theory that results in equation (9) is valid only if the texture lengthscales are small compared to the characteristic lengthscales of near-wall turbulence, so that the near-wall cycle is not altered. For a given permeable material (i.e. with fixed permeability values $K_{x}$, $K_{z}$ and $K_{y}$), the permeabilities $K_{x}^{+}$ and $K_{z}^{+}$ in viscous units would increase as the friction Reynolds number increases, thereby increasing $\Delta U^{+}$. As $K_{x}^{+}$ and $K_{z}^{+}$ increase, equation (9) would eventually stop holding, as other mechanisms set in, degrading the drag-reducing performance.

3.4 Onset of Kelvin-Helmholtz rollers

Equation (9) does not explicitly include the wall-normal permeability, or transpiration in general. However, most complex surfaces that produce slip produce also a non-zero wall-normal velocity at the reference plane, such as permeable substrates (Breugem & Boersma, 2005), riblets (García-Mayoral & Jiménez, 2011) or superhydrophobic surfaces (Seo et al., 2018), and this effect induces generally a degradation in drag. Abderrahaman-Elena & García-Mayoral (2017) and Gómez-de-Segura et al. (2018b) argued that the onset of the Kelvin-Helmholtz instability discussed in the introduction would disrupt the linear regime of equation (9), and could therefore be used to establish an a priori limit for its range of validity.

Kelvin-Helmholtz rollers are ubiquitous over permeable substrates and are known to increase drag (Breugem et al., 2006; Kuwata & Suga, 2016; Suga et al., 2017). Abderrahaman-Elena & García-Mayoral (2017) developed a model based on Darcy’s equation to characterise the onset of such structures, which is well-suited for poorly-connected permeable media. Gómez-de-Segura et al. (2018b) extended their analysis for highly-connected permeable substrates, which have a greater potential for drag reduction, and showed that the latter exhibit a different behaviour for the onset of the instability.

In this section, we summarise the procedure and results from Gómez-de-Segura et al. (2018b). The procedure is based on a linear stability analysis on the mean turbulent profile to capture the onset of Kelvin-Helmholtz rollers, as in Jiménez et al. (2001), García-Mayoral & Jiménez (2011) and Abderrahaman-Elena & García-Mayoral (2017). The analysis is restricted to spanwise-homogeneous modes, as Kelvin-Helmholtz rollers are predominantly spanwise coherent. Considering normal-mode solutions of the form $v^{\prime}=\hat{v}(y)\exp(i(\alpha x-\omega t))$, where the wavenumber $\alpha$ is real and the angular frequency $\omega=\omega_{r}+i\omega_{i}$ is complex, the Orr-Sommerfeld equation with a variable eddy viscosity in $y$ (Cess, 1958) can be solved. Modes are unstable if $\omega_{i}$ is positive. At the substrate-channel interface, equation (3) for $\hat{u}$ and $\hat{v}$ was imposed, particularising for $\alpha_{z}=0$ and $\hat{w}=0$.

The presence of permeable substrates at the boundaries of the channel destabilises the otherwise stable mean flow, in agreement with Jiménez et al. (2001), Tilton & Cortelezzi (2008) and Abderrahaman-Elena & García-Mayoral (2017). The most unstable mode forms counterrotating rollers separated by a wavelength $\lambda^{+}\approx 70$, as shown in figure 6, which resemble Kelvin-Helmholtz rollers. García-Mayoral & Jiménez (2011) and Abderrahaman-Elena & García-Mayoral (2017) argued that the height at which the energy-producing term of the Orr-Sommerfeld equation, $d^{2}U/dy^{2}$, concentrates, $y_{c}^{+}\simeq 9$, sets the lengthscale for the instability. This height is essentially independent of the Reynolds number when measured in wall units, resulting in the optimum $\lambda^{+}\sim 2\pi y_{c}^{+}$ regardless of the topology of the substrate.

While the substrate topology does not significantly alter the wavelength of the most amplified mode, it determines its amplification. Gómez-de-Segura et al. (2018b) proposed a single, empirically-fitted parameter to capture the effect of the topology on the amplification, given for streamwise-preferential substrates by

[TABLE]

Figure 7 shows how the amplification for different substrates is essentially a function of $K_{Br}^{+}$ only. From an application point of view, we are interested in sufficiently deep and streamwise-preferential substrates, h^{+}\gtrsim\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}\gg\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}. If this is the case, the second hyperbolic tangent in equation (10) is approximately $1$. For \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}\gtrsim 5, the first hyperbolic tangent is also approximately $1$, and $K_{Br}^{+}$ becomes

[TABLE]

Hence, the amplification of the instability is mainly determined by $K_{y}^{+}$, and $h^{+}$ and $K_{x}^{+}$ have only a secondary effect.

Depending on the value of $K_{Br}^{+}$, Gómez-de-Segura et al. (2018b) hypothesised three regimes for the instability, as shown in figure 7: a low-permeability regime, \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}\lesssim 1, where the instability would be weak and not expected to emerge in the flow; a high-permeability regime, \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}\gtrsim 2.2, where the amplification approaches an asymptote and the instability would be fully developed; and an intermediate regime, where the instability would set in. García-Mayoral & Jiménez (2011) found that, for riblets, the instability sets in for amplifications of approximately half the maximum. Following this, Gómez-de-Segura et al. (2018b) hypothesised that the intermediate regime would occur for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx\sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}=1-2.2, and the linear drag reduction of equation (9) could only hold for lower values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}. This hypothesis will be re-assessed based on the present DNS results in §5.4.

3.5 Theoretical prediction of drag-reducing curves

Combining the information on the linear drag reduction of equation (9) and the range of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} for the onset of Kelvin-Helmholtz rollers, the trend of the drag reduction curves for anisotropic permeable substrates can be estimated (Abderrahaman-Elena & García-Mayoral, 2017; Gómez-de-Segura et al., 2018b). An optimal substrate should seek to maximise the difference \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}-\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}} to obtain a large slip effect, while maintaining \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} as low as possible to inhibit the appearance of Kelvin-Helmholtz rollers. Fibrous substrates as those proposed in figure 2(b) would conform such a material, for instance.

In this study, a substrate configuration will be represented by three dimensionless parameters; the anisotropy ratios $\phi_{xy}=\sqrt{K_{x}/K_{y}}$ and $\phi_{zy}=\sqrt{K_{z}/K_{y}}$, and the dimensionless thickness, h/\sqrt{K_{\;\;}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}. Given that both $K_{y}^{+}$ and $K_{z}^{+}$ have an adverse effect on the drag, in what follows we consider materials with preferential permeability in $x$ and equally low permeabilities in $y$ and $z$, $K_{x}^{+}>K_{z}^{+}=K_{y}^{+}$. This implies $\phi_{xy}>1$ and $\phi_{zy}=1$. In addition, we consider deep substrates with large h/\sqrt{K_{\;\;}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, so that the substrate thickness does not affect the overlying flow. In §5, we study substrates with \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}\lesssim 10, so a thickness $h^{+}\gtrsim 50$ would suffice. In practical aeronautic applications, for instance, this would imply permeable layers with sub-millimetre thickness.

For substrates with $\phi_{zy}=1$, the expression for $\Delta U^{+}$ in equation (9) becomes

[TABLE]

The drag reduction for a given substrate configuration (i.e. for a fixed $\phi_{xy}$) can then be expressed solely as a function of the wall-normal permeability lengthscale, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, which can be interpreted as a substrate Reynolds number, as sketched in figure 8(a). In a wind-tunnel experiment, for instance, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} could be changed by changing the friction velocity, while $\phi_{xy}$ remained unaltered for a given substrate.

From the present analysis, the resulting drag reduction curves would be analogous to those for riblets (Bechert et al., 1997; García-Mayoral & Jiménez, 2011). The curves would exhibit a linear increase of $\Delta U^{+}$ with \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, breaking down no later than in the shaded region in the figure, where the onset of the Kelvin-Helmholtz instability would be expected. According to equation (12), the slope in the linear regime is predicted to depend on $\phi_{xy}$, and the maximum $\Delta U^{+}$ for a given $\phi_{xy}$ would be determined by the intercept of the corresponding curve with the shaded region. The exact value of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} for the breakdown, as well as the form of the curves in its proximity and for larger values will be obtained from the DNSs presented in §5.

The ideas illustrated in figure 8(a) for a few substrate configurations can be summarised for a wide range of anisotropy ratios, as is done in figure 8(b). Following a drag reduction curve as \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} increases in figure 8(a) would be equivalent to ascending vertically along a constant-$\phi_{xy}$ line in figure 8(b). The linear drag-reducing behaviour of equation (12) is expected to begin to fail in the shaded region. This shaded region represents the permeability values for which the drag-degrading Kelvin-Helmholtz rollers are expected to appear and is the same as in figure 8(a). It is determined by introducing the limiting values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} specified in §3.4, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{lim}\approx 1-2.2, into equation (12). Although additional adverse phenomena cannot be ruled out, figure 8(b) allows us to bound the parameter space for realisable drag reduction. This figure was used to select the region in the parametric space subsequently investigated in §5. Most cases studied are in the drag-reducing region, where equation (12) is expected to hold, and a few cases have been selected in the shaded region, to capture the breakdown. We have considered three substrate configurations, $\phi_{xy}=\sqrt{13}$, $\sqrt{30}$ and $\sqrt{130}$, represented by the three vertical lines of symbols in figure 8(b), and simulated them at different substrate Reynolds numbers, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, so that complete drag reduction curves could be obtained.

3.6 Change in drag in internal and external flows

The expressions for $\Delta U^{+}$ of equations (6) and (9) are valid only for external flows with mild or zero pressure gradients, where the flow near the wall is essentially driven by the overlying shear and the effect of the mean pressure gradient within the permeable substrate is negligible. We are mainly interested in vehicular applications, where the flow falls into that category, but for completion let us discuss the differences with internal flows. In the latter, the effect of the mean pressure gradient could be significant. This effect is essentially additive, so in the following discussion we will leave out turbulence for simplicity, and consider the laminar case.

Sketches of the mean velocity profiles in a boundary layer and in a channel are depicted in figures 9(a) and 9(b), respectively. In a boundary layer over a permeable substrate, there would be a slip velocity at the interface, $U_{Br}$, due solely to the formation of a Brinkman layer within the substrate. It follows from equation (8a) that, provided that the substrate is sufficiently deep, this slip velocity is U_{Br}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}. Compared with a smooth wall, the only change in the mean velocity profile would be a shift by $U_{Br}$, that is \Delta U^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}, and the drag reduction experienced would arise entirely from this slip effect.

In channel flows, there are two limiting forms of applying the permeable substrates to the reference smooth channel of height $2\delta$. They can substitute a layer of solid material, increasing the height to $2(\delta+h)$, or they can be placed on top of the reference smooth channel, reducing the free flow area. In the first case, depicted in figure 9(b), the mean pressure gradient acts on the region $2(\delta+h)$, which includes the permeable substrates. This produces two opposite effects on the drag: a positive effect due to an increment in the flow rate, not only within the substrate but also in the channel core, and a negative effect due to the pressure gradient being applied across a larger cross-section. In order to evaluate these two effects, we compare the friction coefficient for the permeable and the smooth channel under equal mean pressure gradient $P_{x}$. The integral force balance yields

[TABLE]

where $\tau_{w}$ accounts for the net force applied on the substrates. As we are now solely considering internal flows, we use the conventional bulk velocity, $U_{b}$, to define $c_{f}$,

[TABLE]

where the subscript ‘0’ refers to the smooth channel. The friction coefficient of the smooth channel is therefore defined as $c_{f0}=-2P_{x}\delta/U_{b0}$, and $q=2(\delta+h)U_{b}$ is the mass flow rate. The opposing effects of the increase in cross-section where the pressure gradient acts, $2h$, and the extra flow rate, $\Delta q$, are evidenced in equation (14). The result can be either a drag reduction or a drag increase depending on the values of $\Delta q$ and $h$.

The extra flow rate, $\Delta q$, can be expressed in terms of $K_{x}$. From figure 9(b), $\Delta q$ is

[TABLE]

where, in addition to the slip velocity caused by the overlying shear, $U_{Br}$, as in figure 9(a), there is an extra slip velocity caused by the mean pressure gradient, $U_{Darcy}$, and a resulting extra flow rate within the substrate, $q_{substrate}$. The former is obtained from Darcy’s law, $U_{Darcy}=-P_{x}K_{x}/\nu$, and $q_{substrate}$ is obtained using $U_{Darcy}$ and Brinkman’s velocity within the substrate, as defined by expression (55a). Substituting expression (15) into equation (14), the resulting change in $c_{f}$ depends on $\sqrt{K_{x}}/h$, $h/\delta$ and the Reynolds number. This dependency for a friction Reynolds number $\delta^{+}=180$ is illustrated in figure 10. The figure shows how the beneficial effect of adding a streamwise permeability is opposite to the deleterious effect of the increased area and, for certain substrate geometries, can even outweigh it, resulting in a net drag reduction. Note that in the turbulent case the effect of the spanwise, Brinkman contribution would also need to be included, as given by equation (9).

To better understand the relationship between $K_{x}$ and $h$, expression (14) can be simplified further for $\delta\gg h\gtrsim\sqrt{K_{x}}$. The extra flow rate is then dominated by $\Delta q\approx 2\delta U_{Br}\approx 2\delta\sqrt{K_{x}}\left.dU/dy\right|_{y=0}$, and in a first order approximation, equation (14) simplifies to

[TABLE]

It follows from this equation that, in $(h,\sqrt{K_{x}})$ parameter space, the isocontours of $DR$ are approximately oblique straight lines, as observed in figure 10, and specifically the neutral drag curve is $\sqrt{K_{x}}=3/2\,\,U_{b0}/(dU/dy|_{y=0})\,\,h$, which depends on the friction Reynolds number through $U_{b0}$ and $dU/dy|_{y=0}$. For $\delta^{+}=180$, the zero drag reduction line is given by $\sqrt{K_{x}}\approx 0.15h$, as indicated in figure 10, while for $\delta^{+}=5000$, $\sqrt{K_{x}}\approx 0.007h$.

The above analysis applies to channel flows where the permeable substrates substitute a layer of solid material and shows that, in this case, the drag can be either reduced or increased. If, on the other hand, the permeable coating was added on top of an existing smooth channel, the pressure gradient would still be applied over the whole height of $2\delta$, which includes the permeable coatings, and the resulting friction coefficient would be $c_{f}=c_{f0}/(1+\Delta q/q_{0})^{2}$. In this case, the flow rate would always decrease, $\Delta q<0$, resulting in an increase of drag independently of the values of $\sqrt{K_{x}}$ and $h$.

4 DNS setup

In this section, we present the numerical setup for direct numerical simulations of the domain sketched in figure 9(c), a channel of height $2\delta$ delimited by two identical anisotropic permeable substrates. The presence of the substrates is taken into account through the boundary conditions defined by equation (3), as in the stability analysis in §3.4.

4.1 The numerical method

The channel flow is governed by the incompressible Navier-Stokes equations,

[TABLE]

where the density has been assumed to be unity for simplicity, $p$ is the pressure, $\boldsymbol{u}=(u,v,w)$ the velocity vector and $Re$ the bulk Reynolds number defined as $Re=U_{b}\delta/\nu$, with $U_{b}$ being the bulk velocity in the channel region. The DNS code is adapted from García-Mayoral & Jiménez (2011) and Fairhall & García-Mayoral (2018) and was already used in Gómez-de-Segura et al. (2018b). It solves the incompressible Navier-Stokes equations (17) in a doubly-periodic channel of height $2\delta$, where $\delta=1$ is the distance between the substrate-channel interface and the centre of the channel. All simulations are conducted at a fixed friction Reynolds number $Re_{\tau}=u_{\tau}\delta/\nu=180$ by imposing a constant mean pressure gradient in $y\in[0,2\delta]$. The kinematic viscosity is $\nu=1/2870$ and we use a smooth-wall channel with the same mean pressure gradient as reference.

Although for convenience the present DNSs are conducted in channels, our scope of application is mainly external flows with mild pressure gradients. In a channel, in comparison, there would be an additional flow rate from Darcy’s contribution discussed in §3.6. To allow direct extrapolation to external flows, we simply do not include this contribution when implementing the boundary conditions on the mean flow, that is, mode $(0,0)$, which would be the only Fourier mode affected. This numerical artefact would be equivalent to applying the mean pressure gradient in the channel region only, as depicted in figure 9(c). The drag reduction in the present simulations results then entirely from the slip velocity due to an overlying shear, as in external flows.

The spatial discretisation is spectral in the wall-parallel directions $x$ and $z$, with $2/3$ rule de-aliasing, and uses second-order centred finite differences on a staggered grid in the wall-normal direction. The computational domain is of size $2\pi\times\pi\times 2$ in the streamwise, spanwise and wall-normal directions, respectively. A grid with $192\times 192\times 153$ collocation points with grid stretching is used, which in viscous units gives a resolution of $\Delta x^{+}\approx 5.9$, $\Delta z^{+}\approx 2.9$, and $\Delta y^{+}\simeq 0.3$ near the wall and $\Delta y^{+}\simeq 3$ in the centre of the channel. For the temporal integration, a Runge-Kutta discretisation is used, where every time-step is divided into three substeps, each of which uses a semi-implicit scheme for the viscous terms and an explicit scheme for the advective terms. Discretised this way, the Navier-Stokes equations in (17) result in

[TABLE]

where $\mathrm{L}$, $\mathrm{G}$ and $\mathrm{D}$ represent the discretised laplacian, gradient and divergence operators, respectively, and $\mathrm{N}$ represents the nonlinear, advective operator. The superscript $k=1,2,3$ denotes the Runge-Kutta substep. Hence, the velocities $\boldsymbol{u}^{0}$ and $\boldsymbol{u}^{3}$ correspond to the velocities at time-step $n$ and $n+1$, respectively. Additionally, $\alpha_{k}$, $\beta_{k}$, $\gamma_{k}$ and $\zeta_{k}$ are the Runge-Kutta coefficients for substep $k$ from Le & Moin (1991). In equation (18), the velocity at substep $k$ is expressed in terms of the velocities at the previous substeps, as well as the pressure at that same substep $k$. To solve it, a fractional step method is integrated in each substep (Le & Moin, 1991).

The presence of the permeable substrates is accounted for by the boundary conditions (3), and the coupling between the velocities and the pressure at the interface is implemented implicitly. Following Perot (1993, 1995), the discretised incompressible Navier-Stokes equations from (18) can be represented in matrix form,

[TABLE]

where $\mathbf{u}^{k}$ and $p^{k}$ are the discrete velocity and pressure unknowns, respectively. $\mathrm{A}$ is the operator containing the implicit part of the diffusive terms, which for the internal points of the domain equation (18) gives $\mathrm{A}=[I-\Delta t\frac{\beta_{k}}{\mathrm{Re}\mathrm{L}}]$, and the vector $r^{n}$ is the explicit right-hand side, which contains all the quantities from previous time-steps. The boundary conditions given by equations (3) are embedded in the block matrices in equation (19). The relationships between the three velocities and the shears $d\hat{u}/dy$ and $d\hat{w}/dy$ are embedded in $\mathrm{A}$, while the coupling between the velocities and pressure is embedded in $\mathrm{A}$ and $\mathrm{G}$. Taking then the LU decomposition of system (19) results in

[TABLE]

and the operations are solved in the following order

[TABLE]

where the variable $\mathbf{u^{*}}$ is an intermediate, non-solenoidal velocity. The Poisson equation in equation (21b) is computationally expensive, as it requires the inversion of matrix $\mathrm{A}$. For efficiency, $\mathrm{A^{-1}}$ is generally approximated to its first order term, $\approx\mathrm{I}$ (Perot, 1993). In the present work, we approximate the internal points in $\mathrm{A}$ by $\approx\mathrm{I}$, while keeping the rows of $\mathrm{A}$ that contain the boundary conditions unchanged, and then invert the resulting matrix to obtain $\mathrm{A}^{-1}$.

Statistics are obtained by averaging over approximately $100$ eddy-turnovers, once the statistically steady state has been reached. Statistical convergence was verified using the criterion of Hoyas & Jiménez (2008).

4.2 Validation

We validate the present Brinkman model with one of the cases studied by Breugem et al. (2006), where the authors used the VANS equations within the permeable substrate. We consider their case E80, here referred to as BB_E80, with a porosity $\epsilon=0.8$ – which refers to the ratio between the void volume and the total volume of the substrate – and an isotropic permeability $K^{+}\approx 1$. This permeability is of the same order of magnitude as our largest permeabilities $K_{y}^{+}$ and $K_{z}^{+}$ in the DNSs presented in §5. The case BB_E80 is compared to our Brinkman model, here referred to as BB_Br, using approximately the same value of permeability, $K^{+}\approx 1$.

To match the validation domain in Breugem et al. (2006), we use an asymmetric channel of height $2\delta$, delimited by an impermeable wall at the top and a permeable substrate at the bottom, as sketched in figure 11. The thickness of the permeable layer is $h=2\delta$.

Breugem et al. (2006) used a VANS approach to model the flow within the substrate. At the interface with the free flow, they let the porosity, and hence the permeability, to evolve gradually from the inner value $\epsilon=0.8$ to the free flow value $\epsilon=1$ over a thin interfacial layer of thickness $\delta_{i}$. This corresponded to the averaging volumes in VANS capturing varying proportions of free flow and substrate, so that the volumes centred at the top of the interfacial layer did not contain any substrate, and vice versa, as illustrated in figure 11. This model is consistent with applying VANS on a setup with a sharp interface half-way through the interfacial layer. We set our reference plane $y=0$ at this height, so for comparison we represent the results from Breugem et al. (2006) in the same frame of reference. Note that in Breugem et al. (2006) the reference plane was at the top of the interfacial region instead. For a consistent comparison, the results from Breugem et al. (2006) have been rescaled with the friction velocity measured at our $y=0$ and with the bulk velocity integrated between that plane and the top impermeable smooth wall.

Both BB_Br and the original simulation of Breugem et al. (2006) were run at a constant mass flow rate started from a smooth channel, at $Re=2750$ in the latter case and $Re=2832$ in ours. Defining viscous units using the friction velocity at $y=0$, the initial friction Reynolds numbers were $\delta^{+}=176$ and $\delta^{+}=180$, respectively, while in the final statistically-steady state they were $\delta^{+}\approx 204$ and $\delta^{+}\approx 203$, respectively. Results from Breugem & Boersma’s VANS approach (BB_E80) and Brinkman’s model under study (BB_Br) are compared in table 1 and figure 12, all showing good agreement.

This agreement between VANS and Brinkman’s approach could be expected, given the similarities between the models for the values of the parameters considered. VANS equations can be viewed as Brinkman’s equation with the addition of the advective and temporal terms, with $\epsilon$ playing in the former the role that $\nu/\tilde{\nu}$ plays in the latter. For small permeabilities, such as those under consideration, the advective terms can be neglected. In addition, an order of magnitude analysis shows that the temporal term can also be neglected. This term is of the order $\sim O[u_{c}/t_{c}]$, where $t_{c}$ and $u_{c}$ denote a characteristic time and velocity, respectively. When the substrate is isotropic and highly-connected (i.e. $\tilde{\nu}\approx\nu$), both the Brinkman and Darcy terms are of the same order of magnitude, as the penetration length in an isotropic medium is of order $\sim\sqrt{K}$. Comparing the temporal and the Brinkman terms then, we obtain

[TABLE]

For the temporal term to be negligible, the characteristic time should satisfy $t_{c}^{+}>K^{+}$. Considering that the fastest-evolving turbulent structures near the wall are the quasi-streamwise vortices, with a radius $r^{+}\sim 15$ and velocity $\sim u_{\tau}$, the smallest characteristic timescale would be $t_{c}^{+}\sim 15$, and given that $K^{+}\sim 1$, the condition $t_{c}^{+}>K^{+}$ is satisfied. The flow within the permeable medium can then be assumed to be quasi-steady. Additionally, for VANS and Brinkman’s equation to be equivalent, the value of the porosity $\epsilon$ should be equal to the ratio $\nu/\tilde{\nu}$. In these simulations, these values differ slightly, $\epsilon=0.8$ in BB_E80 and $\nu/\tilde{\nu}=1$ in our model. Nonetheless, Rosti et al. (2015) reported that, for porosity values beyond $\epsilon\gtrsim 0.6$, a further increase of $\epsilon$ had no significant effect on the overlying flow, and the permeability $K$ was the only key parameter. This justifies the similarities between the results of the two models, even if the values of $\epsilon$ and $\nu/\tilde{\nu}$ are not exactly matching.

5 Results and discussion for DNS

In this section, we present results from DNSs to investigate in detail the effect that permeable substrates have on the overlying flow and assess the validity of the predictions presented in §3. We study three substrate configurations, given by three different anisotropy ratios $\phi_{xy}\approx 3.6$, $5.5$, $11.4$. For our main set of simulations, the substrates have thickness h=100\sqrt{K_{\;\;}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, large enough for the problem to become independent of it, and the same permeabilities in $y$ and $z$, $\phi_{zy}=1$. An additional subset of simulations was conducted to explore the effect of a finite $h$ on the substrate performance. For a given configuration (i.e. a fixed $\phi_{xy}$ and h/\sqrt{K_{\;\;}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}), we vary proportionately the permeabilities in viscous units, $K_{x}^{+}$, $K_{z}^{+}$ and $K_{y}^{+}$, which is equivalent to varying the viscous length. For each configuration, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} varies between $0.7-11$. The simulations under study are summarised in table 2, where each case is labelled with a letter and a number. In the main set of simulations, the letter refers to the anisotropy ratio $\phi_{xy}$ of the substrate and the number to the specific substrate, with fixed permeabilities in viscous units. In the secondary set, an additional subscripts ′, ′′ and ′′′ indicate decreasing substrate depth.

The virtual-origin model presented in §3 is based on the idea that the near-wall cycle remains smooth-wall-like, other than by being displaced a depth $\ell_{T}$ towards the substrate. Given that the origin perceived by turbulence is expected to be at $y=-\ell_{T}\approx-\sqrt{K_{z}}$, throughout this section results are scaled taking that as the reference for the wall-normal height. The friction velocity is obtained by extrapolating the total stresses to that height, $u_{\tau}=u_{\tau_{y=0}}\left(1+\sqrt{K_{z}}/\delta\right)^{1/2}$, and the effective half-height channel becomes $\delta^{\prime}=\delta+\sqrt{K_{z}}$ (García-Mayoral et al., 2019), although the effect is negligible for the small values of $\sqrt{K_{z}}/\delta$ considered here. Beyond the breakdown of the linear regime, the virtual-origin model begins to fail and the effect of the substrates can no longer be solely ascribed to a shift in origins. Nonetheless, for the cases lying in the degraded regime, we still use the virtual origin that would be valid in the linear regime, $y=-\sqrt{K_{z}}$, to measure $u_{\tau}$. In this framework, any further effect can be interpreted as additive. The values of $\Delta U^{+}$ have been obtained using this $u_{\tau}$ and comparing with a smooth-wall velocity profile with the origin shifted to $y=-\sqrt{K_{z}}$, although the effect of the shift on $\Delta U^{+}$ is also negligible.

5.1 Drag reduction curves

The drag reduction curves obtained from the main set of DNSs are shown in figure 13. For small permeabilities, a linear drag-reduction regime is observed. In §3.3, we predicted $\Delta U^{+}$ in this regime to be equal to the difference between the virtual origin for the mean flow, $\ell_{U}^{+}$, and that perceived by turbulence, $\ell_{T}^{+}$. For the substrates under consideration, these would be \ell_{U}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} and \ell_{T}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}, as given by equation (9). This prediction agrees well with the DNS results, and the three substrate configurations exhibit roughly the same initial unit slope in figure 13(b). The breakdown of the linear drag reduction, however, occurs for different values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}-\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}} depending on the substrate.

In contrast, when the lengthscale is represented using \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} – the parameter predicted in §3.4 to trigger the Kelvin-Helmholtz instability – the location of the breakdown coincides for all the curves, as shown in figure 13(c). For all substrate configurations, the drag reduction is maximum for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{opt}\approx 0.38 and the drag becomes greater than for a smooth wall for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\gtrsim 0.6.

The common linear drag reduction behaviour, observed in figure 13(b), and its common breakdown, observed in figure 13(c), are condensed in figure 13(d). This is done by dividing $\Delta U^{+}$ from figure 13(c) by the slope for each curve expected from equation (12), $\phi_{xy}-1$. Given that in this equation $\Delta U^{+}$ depends only on $\phi_{xy}$ and \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, the general collapse suggested by this figure could be used to predict the performance of permeable substrates different to those explored in this work. Considering that the maximum $\Delta U^{+}$ in figure 13(d) occurs for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{opt}\approx 0.38 and is approximately $80\%$ of that estimated by equation (12), the maximum $\Delta U^{+}$ would depend only on the anisotropy ratio,

[TABLE]

For substrates with different cross permeabilities, $\phi_{zy}\neq 1$, it follows from equation (9) that the maximum $\Delta U^{+}$ would be $\Delta U^{+}_{max}\approx 0.8\times 0.38\times(\phi_{xy}-\phi_{zy})$.

The secondary set of simulations aims to explore the effect of the substrate depth on $\Delta U^{+}$, and to test if the performance could be improved by reducing the depth enough for it to become a parameter in the problem. For this, the same substrate of cases C1-C7 is studied with depths $h/\sqrt{K_{x}}=1.5$, $1.0$, and $0.5$. From equations (8), we can expect shallower substrates to have smaller $\ell_{U}^{+}$ and $\ell_{T}^{+}$, as the hyperbolic tangent terms become smaller than unity. This would reduce the slope of the $\Delta U^{+}$ curve in the linear regime and be an adverse effect. However, a reduced depth would also have the beneficial effect of making the substrate more robust to the onset of Kelvin-Helmholtz-like rollers, as at a given Reynolds number (i.e. \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}) equation (10) would predict a smaller \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}. Note also that \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}} is a parameter empirically fitted to the results from the linear stability model, and that the actual results in §3.4 show that shallower substrates have in fact a delayed onset in terms of \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}, as shown in figure 7.

The results for $\Delta U^{+}$ for the shallow substrates of the secondary set of simulations are portrayed in figure 14, compared with the corresponding deep substrate from the main set, cases C1-C7. Given that all our substrates have higher permeability in $x$, the first terms to experience the effect of a finite $h$ in equations (8) and (10) are those where $h$ appears scaled with $\sqrt{K_{x}}$. Note that if we had considered values of $h$ small enough for $h/\sqrt{K_{z}}$ to be also small, we would have $\ell_{U}^{+}\approx\ell_{T}^{+}\approx h^{+}$, which would yield no drag-reducing effect. For the values of $h/\sqrt{K_{x}}$ studied, we have $h/\sqrt{K_{y}}=h/\sqrt{K_{z}}=6$, $11$ and $17$, so the corresponding hyperbolic tangent terms in equations (8) and (10) are still essentially unity. This can be appreciated for instance in figure 14(a), where the predicted slope in the linear regime has been adjusted for the effect of $h^{+}$ on the streamwise slip, \ell_{U}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}\tanh(h^{+}/\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}), but the spanwise slip remains \ell_{T}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}. Figure 14(b), however, shows that \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} is no longer adequate to parametrise the onset of the drag degradation. Panel (c), in turn, suggests that a suitable alternative is \sqrt{K^{\prime}\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}=\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\tanh(h^{+}/(9\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}})), and that the optimum value is still \sqrt{K^{\prime}\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}\approx 0.38, as in figure 13. All the curves can be once more collapsed by reducing $\Delta U^{+}$ with its predicted slope in the linear regime and expressing the Reynolds number in terms of \sqrt{K^{\prime}\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}, as is done in panel (d). This suggests that the optimum performance for shallow substrates can also be predicted and would be \Delta U^{+}_{max}\approx 0.8\times 0.38\times[\phi_{xy}\tanh(h/\sqrt{K_{x}})-\phi_{zy}]/\tanh(h/9\sqrt{K_{\;\;}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}). Note, however, that $\Delta U^{+}_{max}$ decreases slightly as the substrate depth is reduced, as can be appreciated in panel (a), and that even if there is a delay in the critical \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} in absolute terms, as observed in panel (b), any gain in the relative width of the ‘drag bucket’ region – the near-optimal range – is insignificant, as is clear from panel (d).

5.2 Flow statistics

To explore the underlying mechanisms for the behaviour observed in the drag reduction curves, let us focus on a fixed substrate configuration, that is on one of the curves in figure 13. Let us take the one with the anisotropy ratio $\phi_{xy}\approx 11.4$, that is, simulations C1-C7. The corresponding data for the other two substrate configurations can be found in Appendix B. To illustrate how the overlying turbulence is modified at different points along the drag reduction curve, figure 15 shows instantaneous realisations of $u$, $v$ and $p$ in an $x$-$z$ plane immediately above the substrate-channel interface. For small \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, the flow field resembles that observed over a smooth wall. This is shown in panels (a-c) and (d-f), where the $u$-field displays the signature of near-wall streaks, and the $v$-field that of quasi-streamwise vortices. As \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} increases beyond the linear regime, the flow begins to be altered, as shown in panels (g-l). Some spanwise coherence emerges, becoming more prevalent for larger \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}. Eventually, the flow becomes strongly spanwise-coherent and no trace of the near-wall cycle remains, as shown for a drag-increasing case in panels (m-o).

To assess quantitatively to what extent turbulence differs from that over smooth walls, we first focus on the one-point statistics resulting from the DNSs, portrayed in figures 17 and 18. The former shows the mean velocity profiles. In panel (a) the results are represented with the origin for the wall-normal height at the substrate-channel interface, $y^{+}=0$, as is typically done in the literature. In this representation, the non-zero slip velocity at the interface, $U_{slip}^{+}$, is apparent at $y^{+}=0$, while far away from the wall the adverse effect of $\ell_{T}^{+}$ and the ‘roughness-like’ shape of the profile, that is the deviation from a smooth-wall-like shape, combine with $U_{slip}^{+}$ to yield the net velocity offset. In this framework, the effect of $\ell_{T}^{+}$ and the deviation from the shape of a smooth-wall profile cannot be easily disentangled.

If the velocity profiles are represented with the origin for the wall-normal height at $y^{+}=-\ell_{T}^{+}$ and if turbulence remained smooth-wall like, the profiles could then be expected to be like those for smooth walls, save for the offset given by equation (6). Subtracting that offset would then give a collapse of all the velocity profiles, and any deviation can then be separately attributed to modifications in the turbulence (García-Mayoral et al., 2019). In figure 18(b), the profiles are portrayed with the origin at $y^{+}=-\ell_{T}^{+}$ and with the offset subtracted from the velocities. For cases C1-C3, which lie in the linear regime, the resulting collapse is indeed good, but beyond this regime the profiles deviate from the smooth-wall behaviour increasingly. Let us note that defining $u_{\tau}$ at $y^{+}=-\ell_{T}^{+}$ implies that the wall-normal gradient of the mean profile at the interface is no longer necessarily $dU^{+}/dy^{+}|_{y^{+}=0}=1$. This is because the stresses at that height in viscous units sum slightly less than one, and more specifically, because a non-zero transpiration gives rise to a Reynolds stress at the interface, so the viscous stress is no longer the only contribution to the total. As a result, $U_{slip}^{+}$ and $\ell_{x}^{+}$ do not strictly have equal value and cannot be used interchangeably. For small values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, the Reynolds stress at the substrate-channel interface is negligible, so this effect is small and $U_{slip}^{+}\approx\ell_{x}^{+}$. This is the case for the substrates lying on the linear regime. However, as \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} increases, the Reynolds stress at the interface ceases to be negligible, and $U_{slip}^{+}=dU^{+}/dy^{+}|_{y^{+}=0}\,\ell_{x}^{+}<\ell_{x}^{+}$. This discrepancy between $U_{slip}^{+}$ and \ell_{x}^{+}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} for the substrates under consideration is shown in figure 17. The effect is small for the substrate of simulations C1-C7, but is significant for the substrates of B1-B7 and A1-A8, with results portrayed in Appendix B. The effect is particularly intense for the latter substrate, which reaches \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 3 and experiences significant transpiration. Although $U_{slip}^{+}$ and $\ell_{x}^{+}$ represent essentially the same concept, the quantitative effect of the streamwise slip is carried more accurately by $U_{slip}^{+}$, so the latter has been used for the velocity offset in figure 18(b). Notice that this effect is negligible in slip-only simulations or other idealised surfaces where zero transpiration is assumed (Fairhall et al., 2019).

The observations on the agreement or deviation from smooth-wall data in the mean velocity profiles extend also to the rms velocity fluctuations and streamwise vorticity, as well as the Reynolds shear stress, portrayed in figures 18(a-e). For the cases in the linear regime, the agreement with smooth-wall data is good. The only difference is a small deviation in the profile of $u^{\prime+}$ in the region immediately above the interface. This deviation is caused by the streamwise velocity effectively tending to zero at $y^{+}=-\ell_{U}^{+}$, below the reference height $y^{+}=-\ell_{T}^{+}$, and essentially does not alter near-wall dynamics (Gómez-de-Segura et al., 2018a). Beyond the linear regime, the fluctuations of the streamwise velocity decrease in intensity, while those of the transverse components increase. For rough surfaces, this is often associated with a decreased anisotropy of the fluctuating velocity (Orlandi & Leonardi, 2006). The Reynolds stress behaves analogously, and the rms streamwise vorticity also becomes more intense, but experiences a significant drop for the final case, C7. The snapshots of figure 15 could suggest that this is caused by the eventual annihilation of the quasi-streamwise vortices of the near-wall cycle, as the spanwise-coherent structures become prevalent.

In the models proposed in §3, the streamwise, spanwise and wall-normal permeabilities have separate effects. These models capture leading-order features, but in equations (3) the effect of the three permeabilities is coupled. This manifests in the DNS results and, although the coupled effects are secondary, they become increasingly important for large permeabilities.

The leading-order effect of the substrate on the overlying turbulence is, as discussed above, set by the transverse permeabilities. Although in the present study they are equal, it could be expected that \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}} governed the virtual-origin effect, while \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} governed the onset of spanwise-coherent dynamics. However, once \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} becomes sufficiently large, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} plays a secondary role by indirectly modulating the transpiration. Quantitatively, this influence is embedded in equations (3). In essence, the wall-normal flow that penetrates into the substrate is in a first instance impeded by \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, but from continuity it eventually needs to traverse the substrate tangentially, being then impeded by \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}, before it leaves through the interface elsewhere. Thus, a large \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} amplifies the transpiration effect of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} or, rather, a small \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} limits it. This can be observed by comparing the three substrates studied at roughly equal values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}. As they have different anisotropy ratios, for the same \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} they have different \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}. Examples are shown in figure 19. The values \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.2, $0.4$ and $1.0$ have been chosen to observe the secondary effect of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} in the linear regime, near the optimum drag reduction, and in the fully degraded regime, respectively. In the first case, the effect of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} is negligible. The only effect is essentially that of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}} setting the virtual origin, and all the one-point statistics show good agreement with smooth wall data. The effect is still small near the optimum, for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.4, but the modulation by \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} begins to manifest, amplifying the effects of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} already discussed above, such as the decreased anisotropy of the velocity fluctuations. Nevertheless, the Reynolds stress curve, and thus the shape of the mean velocity profile, remain close to those in the linear regime and for smooth walls. In the fully-degraded regime, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 1.0, the modulating effect of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} becomes stronger and results in a further degradation of the Reynolds stress, the mean profile and the drag. The near-wall cycle is severely disrupted in this regime, and the main effect of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} on the velocity fluctuations is on $u^{\prime+}$ near the wall, directly through the increased streamwise slip.

In turn, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} also has a secondary effect on the streamwise slip, through the non-zero Reynolds stress at the interface discussed above. Figure 17 illustrates how, for the same \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}, which governs $U_{slip}^{+}$ to leading-order, substrates with larger \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} have a smaller slip velocity.

While the analysis of the one-point statistics reveals variations in average intensities at different heights, it cannot provide information on whether those variations are caused by contributions from lengthscales that are not active over smooth walls, or from a change in the intensity of the typical lengthscales of canonical wall turbulence. To investigate this, we analyse the spectral energy distribution of the fluctuating velocities.

As an example, spectral density maps of $u^{2}$, $v^{2}$, $w^{2}$ and $uv$ are represented at a height of roughly 15 wall units above the virtual origin for turbulence in figure 20. For substrates in the linear regime, such as C2 in panels (a-d), the agreement in spectral distribution with smooth-wall flows is excellent, as it was for the rms values, further supporting the idea that near-wall turbulence remains essentially canonical. For substrate C4, which is just past the linear regime and has a near-optimum \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.32, differences begin to appear in the spectral distributions, like additional energy at slightly shorter streamwise wavelengths, but most notably the emergence of a spectral region with high $v^{2}$ for large spanwise wavelengths, $\lambda_{z}^{+}\approx 200-\infty$ and streamwise wavelengths $\lambda_{x}^{+}\approx 100-200$. This feature is consistent with the onset of spanwise-coherent structures observed in figure 15, and was also observed previously for riblets and connected to the presence of Kelvin-Helmholtz-like rollers (García-Mayoral & Jiménez, 2011). The above effects become more intense for cases C6 and C7. For C6, which lies in the degraded regime but still yields a net reduction in drag, energy appears in wavelengths as short as $\lambda_{x}^{+}\approx 50$, and the spanwise-coherent region spans a wider set of streamwise wavelengths, $\lambda_{x}^{+}\approx 60-350$, although there is still a trace of the spectral densities of smooth-wall flow for long wavelengths, $\lambda_{x}^{+}\gtrsim 500$, specially for $v^{2}$ and $w^{2}$. For case C7, which gives a net drag increase, any residual trace of the spectral distribution for smooth-wall turbulence has disappeared, and the range $\lambda_{x}^{+}\approx 60-350$ becomes dominant in $v^{2}$.

5.3 Contributions to $\Delta U^{+}$

The degradation of the drag reduction curves in figure 13 and the lack of collapse of the mean velocity profiles in figure 17(b) show that there is an additional contribution to $\Delta U^{+}$ beyond the virtual-origins effect predicted in §3.3. To investigate this, we obtain an expression for $\Delta U^{+}$ by integrating the mean streamwise momentum equation for a permeable channel and comparing it with that for a smooth channel. This procedure follows closely MacDonals et al. (2016), Abderrahaman-Elena et al. (2019) and Fairhall et al. (2019), and is similar to that followed by García-Mayoral & Jiménez (2011). The streamwise momentum equation is averaged in time and in the streamwise and spanwise directions, and integrated in the wall-normal direction,

[TABLE]

where the virtual origin of turbulence is taken as the reference for the wall-normal coordinate, i.e. $y^{\prime+}=y^{+}+\ell_{T}^{+}$, and it is also the height where $u_{\tau}$ is measured. The effective half-channel height or the effective friction Reynolds number is then $\delta^{\prime+}=\delta^{+}+\ell_{T}^{+}$, as previously defined. In equation (24), $\overline{u^{\prime}v^{\prime}}^{+}$ is the Reynolds stress, $dU^{+}/dy^{+}$ the viscous stress and the right-hand side represents the total stress. These three terms are represented in figure 21(a).

Integrating again between two heights, the viscous stress term gives the velocity $U^{+}$ at those two heights, which can be compared to the corresponding equation for a smooth channel to obtain an expression for $\Delta U^{+}$. The upper integration limit is then taken at an arbitrary height in the logarithmic region, $y^{\prime+}=H^{+}$, so that the difference in $U^{+}$ yields $\Delta U^{+}$. For the lower limit, we set it at $y^{\prime+}=\ell_{T,P}^{+}$, where $\ell_{T,P}^{+}$ refers to the virtual origin of turbulence for the permeable case, since for that layout equation (24) is defined only above that height. Integrating equation (24) from $y^{\prime+}=\ell_{T,P}^{+}$, to an arbitrary height in the logarithmic region, $y^{\prime+}=H^{+}$, yields

[TABLE]

This equation applies not only to a permeable channel, but also to a smooth channel at the same Reynolds number, $\delta^{\prime+}$, as depicted in figure 21(b). Note that for a smooth channel $y^{\prime+}=y^{+}$, since the origin of turbulence is at the wall, but the lower integration limit can still be set at some height above the wall, $y^{\prime+}=\ell_{T,P}^{+}$, with $\ell_{T,P}^{+}$ referring to the origin of the permeable case being compared.

Subtracting equation (25) for the permeable case and for the smooth channel, the resulting expression for $\Delta U^{+}$ is,

[TABLE]

where subscript ‘$P$’ denotes the permeable channel, and subscript ‘$S$’ the reference smooth channel at the same friction Reynolds number $\delta^{\prime+}$. Equation (26) shows that $\Delta U^{+}$, defined as the difference in $U^{+}$ between a permeable and smooth channel measured at the same distance from their respective origins of turbulence, consists of the sum of three terms.

The first term, is the slip velocity of the permeable case at the substrate-channel interface, $U_{slip}^{+}$. This is a drag-reducing term, and for the cases lying in the linear regime it can be approximated to the virtual origin of the mean flow, $\ell_{U}^{+}$, since $dU_{P}^{+}/dy^{+}|_{y^{\prime+}=\ell_{T,P}^{+}}\approx 1$. The second term, $U_{S}^{+}(\ell_{T,P}^{+})$, is the mean velocity of the smooth channel measured at $y^{\prime+}=\ell_{T,P}^{+}$. It is a drag-increasing term, and if $\ell_{T,P}^{+}\lesssim 5$, it can be accurately approximated as $U_{S}^{+}(\ell_{T,P}^{+})\approx\ell_{T,P}^{+}$. This is essentially the same as the spanwise protrusion height of Luchini et al. (1991) and Luchini (1996), and the spanwise slip of superhydrophobic surfaces (Min & Kim, 2004; Busse & Sandham, 2012). The offset between these terms is then $U_{slip}^{+}-U_{S}(\ell_{T,P}^{+})\approx\ell_{U}^{+}-\ell_{T,P}^{+}$ and represents the virtual-origin effect discussed throughout this paper. Note however that the exact contribution to $\Delta U^{+}$ involves velocities and not virtual origins as pointed out before. The contribution of the offset between these two terms to $\Delta U^{+}$ is shown in figure 22, where we can appreciate that the virtual origin approximation $\ell_{U}^{+}-\ell_{T}^{+}$ is valid not only in the linear regime, but even slightly beyond the optimum.

The third term, $\mathcal{T}_{uv}$, represents the additional Reynolds stress induced by the permeable substrate. It is a drag-increasing term and its contribution to $\Delta U^{+}$ is also shown in figure 22. For the substrates lying in the linear regime, the Reynolds stress is smooth-wall-like, except for the displacement $\ell_{T}^{+}$ towards the interface, and the term $\mathcal{T}_{uv}$ is therefore zero. The contribution of this term begins to be significant at the breakdown \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{opt}, and increases with \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} in the degraded region. An increase in Reynolds stress is therefore responsible for the degradation of the drag-reducing behaviour of permeable substrates.

The spectral energy distribution of the wall-normal velocity in figure 20 shows the appearance of a new spectral region for large spanwise wavelengths centred around $\lambda_{x}^{+}\approx 150$, which is associated to the large spanwise coherent structures observed in figure 15. To explore whether the additional Reynolds stress accounted for by $\mathcal{T}_{uv}$ is due to the energy accumulated in this spectral region, we define a spectral box with $\lambda_{x}^{+}\approx 70-320$ and $\lambda_{z}^{+}\gtrsim 120$, as that depicted in figure 20, and quantify its contribution to the additional Reynolds stress, as in García-Mayoral & Jiménez (2011). The values are also included in figure 22, showing a close agreement with the whole $\mathcal{T}_{uv}$. This suggests that the new spanwise-coherent structures are indeed responsible for the degradation of the drag, as it was also observed for riblets in García-Mayoral & Jiménez (2011). In essence, these structures increase the turbulence mixing, increasing the local Reynolds stress, and consequently the global drag.

Note that equation (26) compares a permeable channel with a smooth one at the same friction Reynolds number. Often, however, a reference smooth channel at exactly the same Reynolds number is not available. This is for instance the case for the simulations presented in this paper, where all the permeable cases are compared to the same smooth channel at a slightly different friction Reynolds number. When the Reynolds numbers match exactly, the total stress, and thus the Reynolds stress, collapse sufficiently far away from the surface, as they approach zero-value at the centre of the channel. The contribution $\mathcal{T}_{uv}$ can then be entirely ascribed to wall effects, that is to the presence of the substrate. If the Reynolds numbers differ, however, there may be a significant contribution to $\mathcal{T}_{uv}$ far from the surface, which is a Reynolds-number effect, rather than a direct effect of the surface. The same effect appears when comparing smooth channels at different friction Reynolds numbers, as illustrated in figure 21(c). To quantify this effect, we compare the smooth channel at Reynolds number $\delta^{\prime+}$ used for equation (26), represented by a subscript ‘$S$’, with another at a different Reynolds number $\delta^{+}$, represented by a subscript ‘$S0$’. Subtracting the two integrated mean streamwise momentum equations, the universality of the near-wall mean velocity profile over smooth walls gives $U_{S}^{+}(\ell_{T,P}^{+})=U_{S0}^{+}(\ell_{T,P}^{+})$ and $U_{S}^{+}(H^{+})=U_{S0}^{+}(H^{+})$, yielding

[TABLE]

When the break-up of equation (26) is applied to DNS results from a complex surface, $P$ in our case, and a smooth wall at a different Reynolds number, $S0$, the integral of the difference in Reynolds stresses would include both the surface and the Reynolds number effects. These, however, can be easily separated as

[TABLE]

Note that, from equation (27), $\mathcal{T}_{Re}$ can be easily calculated a priori as the area of the trapezoid formed between the total stress lines for $\delta^{+}$ and $\delta^{\prime+}$, as highlighted in figure 21(c). $\mathcal{T}_{uv}$ can subsequently be obtained by subtracting $\mathcal{T}_{Re}$ from the integral of the difference in Reynolds stresses of cases $P$ and $S0$, as given by equation (28), so that it only includes the effect of the surface. This has been the procedure used to obtain the results shown in figure 22, even though for the small values of $\ell_{T}^{+}$ considered, the Reynolds number effect, $\mathcal{T}_{Re}$, is negligible.

5.4 Adjustment of the theoretical models

In §3, we presented theoretical models to estimate the drag-reducing behaviour for anisotropic permeable substrates, specifically, a linear drag-reduction model for small permeabilities given by equation (9) and a threshold for the degradation of this linear regime based on the onset of Kelvin-Helmholtz rollers. The information obtained from the present DNSs can be used to assess the validity of the theoretical models summarised in §3, and if necessary adjust them, so that more accurate predictions can be made.

The drag reduction curves in figure 13 show that the linear regime is accurately represented by the offset between the virtual origins of the mean flow and that of turbulence, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}-\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}, as predicted by equation (9). As discussed above, $\Delta U^{+}$ in this regime would be more precisely given by the difference $U_{slip}^{+}-U_{S}^{+}(\ell_{T}^{+})$, but the differences between \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}} and $U_{slip}^{+}$, and between \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}} and $U_{S}^{+}(\ell_{T}^{+})$ only become significant for larger permeabilities – beyond the linear regime, as shown in figure 22.

The DNS results and the discussion in §5.3 also support the idea that the degradation of the drag-reducing behaviour is caused by the formation of spanwise-coherent structures. These are generally associated to a Kelvin-Helmholtz-like instability, as discussed in §3. In that section, we predicted that the onset of these structures was governed by \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, the leading order term of \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}} from equation (10), as shown in figure 7. From this figure, we estimated an a priori threshold for the onset of Kelvin-Helmholtz-like rollers in the range \sqrt{K\raisebox{4.30554pt}{\scriptsize+}_{\;}}\raisebox{-2.58334pt}{\makebox[0.0pt][l]{\scriptsizeBr}}\approx\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 1-2.2, beyond which equation (9) would no longer be valid. The drag reduction curves in figure 13(c), however, show that the degradation sets in for lower values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} than initially hypothesised. The optimum value of $\Delta U^{+}$ occurs at \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{opt}\approx 0.38, after which performance degrades, and drag becomes greater than that for smooth walls for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}|_{\Delta U^{+}=0}\approx 0.6. Adjusting figure 7 to account for these observed values, we obtain figure 23(a), which shows that the onset occurs as soon as the predicted amplification of the instability becomes positive.

In §3, combining the equation for the linear regime with the limiting values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}, allowed us to design the parameter space for realisable drag reduction shown in figure 8(b), which later served to select the DNS cases studied in §5. Using the limiting values of \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}} observed in the DNSs (\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.38-0.6), the adjusted prediction map for $\Delta U^{+}$ is shown in figure 23(b), where the actual values of $\Delta U^{+}$ measured from DNSs are also shown. This figure illustrates how the theoretical predictions compare to the actual results obtained from DNS. In the linear regime, $\Delta U^{+}$ is well predicted by the theory. At the optimum $\Delta U^{+}$ line, \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.38, the exact value of $\Delta U^{+}$ is given by equation (23), that is, it is roughly $80\%$ of the linear-regime prediction. Beyond this line, the performance degrades, and for the line \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.6, the drag reduction is fully negated. Assuming that this behaviour holds for substrates with anisotropy ratios different to those studied in this work, figure 23(b), which essentially contains the same information of figure 14(d), can be used to estimate their performance.

6 Conclusions

We have explored the ability of anisotropic permeable substrates to reduce turbulent skin friction. We have examined the effect of the streamwise, wall-normal and spanwise permeabilities in highly-connected substrates, and showed that streamwise-preferential substrates can reduce drag.

We have conducted a series of DNSs of turbulent channels delimited by permeable substrates, where the flow within the substrates was modelled using Brinkman’s equation. The resulting drag reduction curves obtained for different substrate configurations (different anisotropy ratios) are similar to the classical curves for riblets: they exhibit a linear drag reduction regime followed by a degradation of performance, eventually leading to an increase of drag.

We have observed that, in the linear regime of small permeabilities, the drag reduction is proportional to the difference between the virtual origin perceived by the mean flow and that perceived by turbulence, which for permeable substrates gives $\Delta U^{+}\approx\sqrt{K_{x}^{+}}-\sqrt{K_{z}^{+}}$. The drag-reducing ability of this technology results therefore from the streamwise-preferential configuration of the substrates, as in other complex surfaces (García-Mayoral et al., 2019). In this regime, the overlying turbulence remains smooth-wall-like, but shifted towards the substrate-channel interface by the origin perceived by turbulence, i.e. \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizez}}.

As permeabilities increase, the linear regime eventually breaks down. We observe that the breakdown is essentially governed by the wall-normal permeability, $K_{y}^{+}$, and occurs for \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.38, independently of the substrate anisotropy. The breakdown can be attributed to the appearance of spanwise-coherent structures, associated to a Kelvin-Helmholtz-like instability. These structures appear to disrupt the near-wall cycle and modify the near-wall turbulence, increasing the Reynolds stress, and consequently, the drag. As permeabilities increase, the drag-increasing, spanwise-coherent structures become prevalent in the flow, outweighing the drag-reducing effect of the streamwise slip and eventually leading to an increase of drag.

In order to predict the drag-reducing behaviour of anisotropic permeable substrates, we have established some theoretical models, which agree well with the behaviour observed from DNS results. The linear regime is accurately described by the expression derived by Abderrahaman-Elena & García-Mayoral (2017), where $\Delta U^{+}=\sqrt{K_{x}^{+}}-\sqrt{K_{z}^{+}}$. This assumes that the permeable medium is highly connected and that the substrate is sufficiently deep for the overlying turbulence not to perceive that its depth is finite, h^{+}\gtrsim 2\sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizex}}. Beyond \sqrt{K\raisebox{4.52083pt}{\scriptsize+}}\raisebox{-1.1625pt}{\makebox[0.0pt][l]{\scriptsizey}}\approx 0.38, the formation of drag-increasing, Kelvin-Helmholtz rollers can be captured with a linear stability analysis. These models provide design guidelines to produce a drag-reducing permeable substrate and give quantitative estimates as to how much drag reduction could be expected.

Further work is nevertheless required to confirm these findings. Direct numerical simulations fully resolving the microstructure of the permeable substrates need to be conducted in order to set the region of validity of the current models, and to gain full understanding on the effect that these substrates have on the overlying turbulence.

GG was supported by an educational grant from Fundación Bancaria ‘la Caixa’, Amelia Earhart Fellowship and an award from The Cambridge Commonwealth, European and International Trust. Some simulations were run using the computational resources provided under EPSRC-UK Tier-2 grant EP/P020259/1 by CSD3, Cambridge.

Appendix A Analytic solution of Brinkman’s equation

The flow within the porous medium is approximated using Brinkman’s equation (1), where $K_{x}$, $K_{y}$ and $K_{z}$ are the principal directions of the permeability tensor and are considered to be different. Together with the continuity equation, the system of equations is

[TABLE]

which can be solved analytically. Here we restrict ourselves to the permeable substrate at the bottom of the channel, which extends from $y=-h$ to $y=0$ and we neglect the influence of a mean pressure gradient within the substrate, as discussed in §3.6.

In order to solve equation (29), we reduce this system of partial differential equation (PDE) with three dependent variables into a single equation with a single dependent variable. We start by taking the divergence of the Brinkman equation (29a)-(29c) and use the continuity equation (29d) to simplify, which yields

[TABLE]

We then take the $y$-derivative of equation (30) and replace $\partial p/\partial y$ from equation (29b) to eliminate the pressure term. Using continuity again to remove the terms in $w$, the following equation is obtained

[TABLE]

which has terms in $v$ and $u$ alone. To remove $u$, we take the $y$-derivative of equation (29a) and subtract the $x$-derivative of (29b). The obtained expression is then differentiated with respect to $x$, yielding

[TABLE]

Substituting for $\partial^{2}u/\partial x\partial y$ from equation (31), a single equation for $v$ is obtained. This equation can be solved by expanding in Fourier series along $x$ and $z$, so that $v(x,y,z)=\hat{v}(y)\mathrm{e}^{i\alpha_{x}x}\mathrm{e}^{i\alpha_{z}z}$, where $\alpha_{x}$ and $\alpha_{z}$ are the wavenumbers, $i$ the imaginary unit, $i=\sqrt{-1}$. Differentiating in $x$ and $z$ becomes then multiplying by $i\alpha_{x}$ and $i\alpha_{z}$, respectively, leading to the following ordinary differential equation (ODE)

[TABLE]

where $D$ denotes $\partial/\partial y$ and $\alpha^{2}=\alpha_{x}^{2}+\alpha_{z}^{2}$. This is a sixth order equation, where all the derivatives are even, and the corresponding characteristic equation is a bicubic equation

[TABLE]

where

[TABLE]

This equation can be reduced to a cubic equation and then solved algebraically. If the discriminant of equation (34) is non-zero, i.e. $\Delta=18a_{3}a_{2}a_{1}a_{0}-4a_{2}^{3}a_{0}+a_{2}^{2}a_{1}^{2}-4a_{3}a_{1}^{3}-27a_{3}^{2}a_{0}^{2}\neq 0$, there are 6 different roots. The roots of the original equation (34) are denoted as $\pm r_{1}$, $\pm r_{2}$ and $\pm r_{3}$ and the general solution for $\hat{v}$ is then

[TABLE]

The constants $A$, $B$, $C$, $D$, $E$ and $F$ are determined once the boundary conditions are imposed and are a function of the geometry and the wavenumbers, $\alpha_{x}$ and $\alpha_{z}$. Similar expressions for the pressure and the streamwise and spanwise velocities can be obtained by substitutions into equations (29b), (31), and the continuity equation (29d), respectively,

[TABLE]

[TABLE]

[TABLE]

To obtain $A$, $B$, $C$, $D$, $E$ and $F$, the boundary conditions need to be considered. The permeable substrate is delimited by an impermeable solid wall at the bottom, where no-slip and impermeability conditions are imposed, and by the free channel flow at the top, where continuity of the normal and tangential stresses holds, together with the continuity of the three velocity components. The boundary conditions at the substrate-channel interface have already been introduced in equation (2). Expanding these boundary conditions in Fourier space, and assuming $\tilde{\nu}\approx\nu$, the continuity of the normal and tangential stresses at the interface simplifies to the continuity of pressure and wall-normal shear ($d\hat{u}/dy$ and $d\hat{w}/dy$), respectively. Thus, the boundary conditions for the permeable substrates are

[TABLE]

where, at $y=0$, the plus and minus signs correspond to the substrate and fluid sides of the interface, respectively.

By applying the above boundary conditions to equations (35),  (37),  (38) and  (36), and particularising the solution at the substrate-channel interface, the velocities at the interface are

[TABLE]

where the coefficients $\mathcal{C}_{ij}(\alpha_{x},\alpha_{z})$ are a function of the wavenumbers, $\alpha_{x}$ and $\alpha_{z}$, and of the geometry of the substrate, i.e. $K_{x}$, $K_{y}$, $K_{z}$ and $h$. An equivalent analysis can be carried out for the upper permeable substrate. Considering the symmetry properties for each variable, this yields

[TABLE]

When $\alpha_{x}=0$ or $\alpha_{z}=0$, Brinkman’s equation simplifies and so does its solution. These cases are solved separately in §A.1, A.2 and A.3.

A.1 Modes $\alpha_{x}\neq 0$, $\alpha_{z}=0$

When $\alpha_{z}=0$, the $z$-derivatives become zero and the Brinkman equation for $w$, i.e. equation (29c), decouples from the other two. The original system of equations simplifies then to

[TABLE]

The velocities $u$ and $v$ can be solved with a procedure similar to that described above, while $w$ can be solved separately.

Taking the two-dimensional divergence of equations (42a) and (42b) in the $(x,y)$ plane and using continuity yields

[TABLE]

Taking the $y$-derivative of equation (42b) and substituting $\partial v/\partial y$ from (43) yields an equation in $\hat{p}$ alone. Taking the Fourier transform in $x$ leads to

[TABLE]

The corresponding characteristic equation is biquadratic,

[TABLE]

Rewriting it as a second order equation, the roots of the characteristic equation are

[TABLE]

Except for the case in which $m_{1}=m_{3}$, i.e. $\frac{K_{x}}{K_{y}}=\frac{1+4\alpha^{2}K_{x}}{4\alpha^{2}K_{x}}$, the expression for $\hat{p}$ becomes:

[TABLE]

where $A^{\prime}$, $B^{\prime}$, $C^{\prime}$ and $D^{\prime}$ depend on the wavenumber $\alpha_{x}$ and the geometrical properties of the permeable medium, and are determined by imposing the boundary conditions – impermeability and no slip conditions at $y=-h$, and continuity of pressure and $d\hat{u}/dy$ at $y=0$. The general solutions for $\hat{v}$ and $\hat{u}$ can be obtained from equations (43) and (42d), respectively.

In contrast, solving equation (42c) for $\hat{w}$ is straightforward. Expanding it in Fourier series gives

[TABLE]

whose solution is

[TABLE]

where $L_{w}=1/\sqrt{\alpha_{x}^{2}+1/K_{z}}$. Applying now the boundary conditions for $\hat{w}$, $\hat{w}=0$ at the impermeable wall and continuity of $d\hat{w}/dy$ at the interface, leads to

[TABLE]

Comparing to the expressions presented in (40), $\mathcal{C}_{ww}(\alpha_{x},0)$ is the proportionality term in equation (49) between $\hat{w}$ and its gradient, whereas $\mathcal{C}_{wu}(\alpha_{x},0)=\mathcal{C}_{wp}(\alpha_{x},0)=0$. Also, from the general solutions for $\hat{u}$ and $\hat{w}$, we observe that $\mathcal{C}_{uw}(\alpha_{x},0)=\mathcal{C}_{vw}(\alpha_{x},0)=0$, which was expected, as there is no coupling between $\hat{w}$, and the other two velocities, $\hat{u}$ and $\hat{v}$, for modes $(\alpha_{x},0)$. Hence, in this case the 9 coefficients presented for the general interface conditions (40) are reduced to only 5.

A.2 Modes $\alpha_{x}=0$, $\alpha_{z}\neq 0$

In this case, Brinkman’s equation for $u$ decouples from the other two. For cases with the same permeability in $y$ and $z$ directions, $K_{z}=K_{y}$, the solution simplifies even more, to

[TABLE]

Taking the two-dimensional divergence of equations (50b) and  (50c) in the $(y,z)$ plane leads to a Laplace equation for the pressure. We then take the Fourier transform with respect to $z$ (i.e. $p(y,z)=\hat{p}(y)\mathrm{e}^{i\alpha_{z}z}$) to get

[TABLE]

The general expressions for $\hat{v}$ and $\hat{w}$ can then be derived from the equations (50b) and  (50d), respectively.

The streamwise velocity is solved similarly to $w$ in §A.1. We take the Fourier transform of equation (50a) in $z$, which gives

[TABLE]

The solution, after applying the boundary conditions for $\hat{u}$, is

[TABLE]

where $L_{u}=1/\sqrt{\alpha_{z}^{2}+1/K_{x}}$. The proportionality coefficient relating $\hat{u}$ with its gradient is the coefficient $\mathcal{C}_{uu}(0,\alpha_{z})$, i.e.

[TABLE]

and $\mathcal{C}_{up}(0,\alpha_{z})=\mathcal{C}_{uw}(0,\alpha_{z})=\mathcal{C}_{wu}(0,\alpha_{z})=\mathcal{C}_{vu}(0,\alpha_{z})=0$.

A.3 Mode $\alpha_{x}=0$, $\alpha_{z}=0$

Although the coefficients for the mean can be directly obtained from the expressions derived in §A.1 and §A.2, this case deserves further discussion. When $\alpha_{x}=\alpha_{z}=0$, the equations for the three velocities $u$, $v$ and $w$ decouple from each other and the velocities for mode zero become

[TABLE]

Equations (55a) and (55b) are obtained from particularising equations (49) and (53) for $\alpha_{z}=0$ and $\alpha_{x}=0$, respectively, while equation (55b) is obtained from continuity, after applying the boundary condition that $\hat{v}=0$ at $y=0$. Particularising at $y=0$ and comparing to the general boundary conditions introduced in equation (40), we have

[TABLE]

where all the coefficients in equation (40) are zero except for $\mathcal{C}_{uu}$ and $\mathcal{C}_{ww}$, which relate the tangential velocities to their wall-normal gradient. These are the mean slip lengths $\ell_{x}^{+}$ and $\ell_{z}^{+}$ derived by Abderrahaman-Elena & García-Mayoral (2017).

Appendix B Turbulence statistics for permeable substrates

In §5 results for only the permeable substrates with $\phi_{xy}\approx 11.4$ are discussed. In this appendix, the flow statistics for the other two substrate configurations are presented. The mean velocity profiles and the turbulence fluctuations for configurations with $\phi_{xy}\approx 5.5$ and $\phi_{xy}\approx 3.6$ are compiled in figure 24.