Drag Reduction by Herringbone Riblet Texture in Direct Numerical Simulations of Turbulent Channel Flow

TL;DR

This study investigates a bird-inspired herringbone riblet texture for turbulent drag reduction using direct numerical simulations, revealing that certain patterns can slightly reduce drag while others increase it significantly due to secondary flows.

Contribution

It extends the drag decomposition method to textured walls and analyzes how herringbone riblet patterns affect turbulent drag, providing insights into their effectiveness.

Findings

Herringbone riblet texture with large spanwise wavelength reduces drag by about 2%.

Small spanwise wavelength increases drag up to 73% due to secondary flows.

Convergent/divergent riblet patterns are generally detrimental to drag reduction.

Abstract

A bird-feather-inspired herringbone riblet texture was investigated for turbulent drag reduction. The texture consists of blade riblets in a converging/diverging or herringbone pattern with spanwise wavelength $\Lambda_f$. The aim is to quantify the drag change for this texture as compared to a smooth wall and to study the underlying mechanisms. To that purpose, Direct Numerical Simulations of turbulent flow in a channel with height $L_z$ were performed. The FIK-identity for drag decomposition was extended to textured walls and was used to study the drag change mechanisms. For $\Lambda_f/L_z \gtrsim O(10)$, the herringbone texture behaves similarly to a conventional parallel-riblet texture in yaw: the suppression of turbulent advective transport results in a slight drag reduction of 2%. For $\Lambda_f/L_z \lesssim O(1)$, the drag increases strongly with a maximum of 73%. This is attributed to enhanced mean and turbulent advection, which results from the strong secondary flow that forms over regions of riblet convergence/divergence. Hence, the employment of convergent/divergent riblets in the texture seems to be detrimental to turbulent drag reduction.

Figures

Tables

Table 1. Table 1: Validation of drag and flow statistics of smooth-wall simulations. The relative deviation of R ​ e τ 𝑅 subscript 𝑒 𝜏 Re_{\tau} from that predicted by Pope’s relation is given by Δ ​ R ​ e τ Pope Δ 𝑅 superscript subscript 𝑒 𝜏 Pope \Delta Re_{\tau}^{\mathrm{Pope}} ( Equation 8 ). Statistics of U 𝑈 U , u rms subscript 𝑢 rms u_{\mathrm{rms}} , v rms subscript 𝑣 rms v_{\mathrm{rms}} , w rms subscript 𝑤 rms w_{\mathrm{rms}} and p rms subscript 𝑝 rms p_{\mathrm{rms}} are compared with simulations by Vreman and Kuerten [ 84 , 85 ] at R ​ e τ ref = 180 𝑅 superscript subscript 𝑒 𝜏 ref 180 Re_{\tau}^{\mathrm{ref}}=180 and 590 590 590 . The relative difference between current and reference flow statistic Q 𝑄 Q is measured with the root-mean-square relative deviation δ ​ Q 𝛿 𝑄 \delta Q ( Equation 9 ).

$R{e}_b$	$R{e}_{\tau }$	$\Delta R{e}_{\tau }^{\mathrm{Pope}}[\%]$	$R{e}_{\tau }^{\mathrm{ref}}$	$\delta U[\%]$	$\delta u_{\mathrm{rms}}[\%]$	$\delta v_{\mathrm{rms}}[\%]$	$\delta w_{\mathrm{rms}}[\%]$	$\delta p_{\mathrm{rms}}[\%]$
5500	175.0	-0.6	180	0.58	1.09	1.39	1.52	3.98
11000	320.9	-1.0	-	-	-	-	-	-
22000	587.4	-1.5	590	0.55	1.85	0.82	0.88	0.94

Table 2. Table 2: Contributions to the total skin friction for smooth-wall turbulent channel flow. The terms are as given in Equation 19 and divided by d tot superscript 𝑑 tot d^{\mathrm{tot}} . The rest term equalizes the left-hand and right-hand side of that equation.

$R{e}_b$	$d^{\mathrm{tot}}[\%]$	$d^{\mathrm{bulk}}[\%]$	$d^{\text{mean advection}}[\%]$	$d^{\text{turbulent advection}}[\%]$	$d^{\mathrm{IBM}}[\%]$	$d^{\mathrm{rest}}[\%]$
5500	100	26.9	$4.3\cdot {10}^{-17}$	73.1	0	-0.013
11000	100	16.0	$1.7\cdot {10}^{-17}$	83.9	0	0.086
22000	100	9.6	$1.4\cdot {10}^{-16}$	90.3	0	0.10

Table 3. Table 3: Parameters of all Direct Numerical Simulations presented in this paper. The simulations are grouped based on the type of texture. The short name indicates which parameters have been varied for a certain texture. The addition (sp. res.) refers to a case with double spanwise resolution; (res.) indicates a double streamwise and spanwise resolution. The parallel riblet variant with N groove = 0.5 subscript 𝑁 groove 0.5 N_{\mathrm{groove}}=0.5 refers to the conventional parallel-riblet texture. The averaging time T 𝑇 T is normalized with δ / u τ 𝛿 subscript 𝑢 𝜏 \delta/u_{\tau} derived from smooth-wall flow at the same bulk Reynolds number.

short name	$s^{+}$	$\alpha$ [∘]	$h/s$	$N_{\mathrm{groove}}$	$L_x$	$L_y$	$N_x$	$N_y$	$N_z$	$\Delta x^{+}$	$\Delta y^{+}$	$\Delta z_w^{}{}_{}{}^{+}$	$\Delta z_c^{}{}_{}{}^{+}$	$N_{cg}$	$R{e}_b$	$T{u}_{\tau }/\delta$	$d^{\mathrm{tot}}\cdot {10}^3$	$DC$ [%]
Smooth Wall
Re5500	-	-	-	-	5.8	2.9	512	512	320	4.0	2.0	0.50	1.7	-	5500	136	8.10 $\pm$ 0.02	0.0 $\pm$ 0.4
Re11000	-	-	-	-	4.9	2.5	800	800	512	4.0	2.0	0.50	2.0	-	11000	66	6.81 $\pm$ 0.02	0.0 $\pm$ 0.5
Re22000	-	-	-	-	5.4	2.7	1600	1600	1024	3.9	2.0	0.49	1.8	-	22000	14	5.70 $\pm$ 0.02	0.0 $\pm$ 0.6
Parallel Riblets
splus10	9.9	0	0.5	-	4.1	2.7	360	960	320	4.0	1.0	0.50	1.8	10	5500	95	7.54 $\pm$ 0.04	-6.9 $\pm$ 0.5
splus17	16.9	0	0.5	-	4.1	2.9	360	1020	320	4.0	1.0	0.50	1.8	17	5500	120	7.49 $\pm$ 0.05	-7.6 $\pm$ 0.7
splus24	23.9	0	0.5	-	4.1	2.7	360	960	320	4.0	1.0	0.50	1.9	24	5500	113	8.34 $\pm$ 0.04	3.0 $\pm$ 0.6
splus24 (sp. res.)	23.9	0	0.5	-	4.1	2.7	360	1920	320	4.0	0.5	0.50	1.9	48	5500	69	8.38 $\pm$ 0.04	3.4 $\pm$ 0.6
splus17 Re11000	16.8	0	0.5	-	4.3	2.9	704	1904	512	4.0	1.0	0.50	2.1	17	11000	39	6.17 $\pm$ 0.03	-9.3 $\pm$ 0.5
splus24 Re11000	23.8	0	0.5	-	4.2	2.7	680	1728	512	4.0	1.0	0.50	2.2	24	11000	50	6.78 $\pm$ 0.04	-0.5 $\pm$ 0.7
splus24 Re22000	24.0	0	0.5	-	4.1	2.6	1200	3072	1024	4.0	1.0	0.50	1.9	24	22000	13	5.61 $\pm$ 0.04	-1.6 $\pm$ 0.8
Parallel Riblets in Yaw
alpha10	16.9	10	0.5	-	4.4	2.5	384	1200	320	4.1	0.7	0.50	1.8	24	5500	124	7.77 $\pm$ 0.05	-4.0 $\pm$ 0.7
alpha15	16.9	15	0.5	-	4.5	2.5	384	800	320	4.1	1.1	0.50	1.8	16	5500	124	8.02 $\pm$ 0.03	-0.9 $\pm$ 0.4
alpha15 (res.)	16.9	15	0.5	-	4.5	2.5	768	1600	320	2.0	0.5	0.50	1.8	32	5500	133	8.01 $\pm$ 0.04	-1.2 $\pm$ 0.5
alpha20	16.9	20	0.5	-	4.5	2.6	512	800	320	3.1	1.1	0.50	1.8	16	5500	124	8.40 $\pm$ 0.04	3.7 $\pm$ 0.6
Herringbone Riblets: $\alpha ={\mathrm{𝟏𝟓}}^{\circ }$
ngroove1	17.0	15	0.5	1	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	77	12.46 $\pm$ 0.06	53.8 $\pm$ 0.8
ngroove4	17.0	15	0.5	4	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	82	14.0 $\pm$ 0.1	73.4 $\pm$ 1.4
ngroove4 (res.)	17.0	15	0.5	4	4.7	3.2	800	2048	320	2.1	0.5	0.50	1.8	32	5500	46	14.24 $\pm$ 0.09	75.7 $\pm$ 1.3
ngroove16	17.0	15	0.5	16	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	72	10.86 $\pm$ 0.05	34.1 $\pm$ 0.7
ngroove128	17.0	15	0.5	128	4.7	12.9	400	4096	320	4.1	1.1	0.50	1.8	16	5500	77	8.29 $\pm$ 0.02	2.3 $\pm$ 0.4
Herringbone Riblets: $\alpha ={\mathrm{𝟏𝟔𝟓}}^{\circ }$
ngroove1	17.0	165	0.5	1	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	77	13.1 $\pm$ 0.1	61.2 $\pm$ 1.5
ngroove4	17.0	165	0.5	4	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	75	11.98 $\pm$ 0.07	47.9 $\pm$ 1.0
ngroove16	17.0	165	0.5	16	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	77	8.71 $\pm$ 0.04	7.5 $\pm$ 0.5
ngroove32	17.0	165	0.5	32	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	73	8.27 $\pm$ 0.04	2.1 $\pm$ 0.5
ngroove128	17.0	165	0.5	128	4.7	12.9	400	4096	320	4.1	1.1	0.50	1.8	16	5500	77	7.94 $\pm$ 0.02	-2.0 $\pm$ 0.4
Herringbone Riblets: $\alpha ={\mathrm{𝟏𝟔𝟓}}^{\circ }$, shifted variant
ngroove1	17.0	165	0.5	1	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	76	13.0 $\pm$ 0.1	60.3 $\pm$ 1.3
ngroove4	17.0	165	0.5	4	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	76	11.9 $\pm$ 0.1	47.2 $\pm$ 1.4
ngroove16	17.0	165	0.5	16	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	77	8.79 $\pm$ 0.03	8.5 $\pm$ 0.5
ngroove32	17.0	165	0.5	32	4.7	3.2	400	1024	320	4.1	1.1	0.50	1.8	16	5500	65	8.27 $\pm$ 0.05	2.1 $\pm$ 0.7
Herringbone Riblets: $\alpha ={\mathrm{𝟎}}^{\circ }$, parallel variant
ngroove0.5	17.0	0	0.5	0.5	4.6	3.1	400	1024	320	4.0	1.1	0.50	1.8	16	5500	81	7.49 $\pm$ 0.03	-7.6 $\pm$ 0.5
ngroove1	17.0	0	0.5	1	4.6	3.1	400	1024	320	4.0	1.1	0.50	1.8	16	5500	71	8.06 $\pm$ 0.03	-0.5 $\pm$ 0.5
ngroove4	17.0	0	0.5	4	4.6	3.1	400	1024	320	4.0	1.1	0.50	1.8	16	5500	75	7.77 $\pm$ 0.04	-4.1 $\pm$ 0.6
ngroove16	17.0	0	0.5	16	4.6	3.1	400	1024	320	4.0	1.1	0.50	1.8	16	5500	78	7.64 $\pm$ 0.08	-5.7 $\pm$ 1.0

Equations

$$\frac{\partial u_i}{\partial t}+\frac{\partial u_i{u}_j}{\partial x_j}$$

$$\frac{\partial u_j}{\partial x_j}$$

$$\overline{\varphi }$$

$${⟨\varphi ⟩}_{x_i}$$

$${⟨\varphi ⟩}_V$$

$$DC$$

$$\equiv \frac{d_{\mathrm{textured}}^{\mathrm{tot}}}{d_{\mathrm{smooth}}^{\mathrm{tot}}}-1.$$

$$u_{DC}^2={\left(\frac{u_{d_t}}{d_s}\right)}^2+{\left(\frac{d_t{u}_{d_s}}{d_s^2}\right)}^2.$$

$$\Delta R{e}_{\tau }^{\mathrm{Pope}}=\frac{R{e}_{\tau }}{R{e}_{\tau }^{\mathrm{Pope}}}-1.$$

$$\delta Q=\sqrt{{⟨{\left(\frac{Q(z)-Q^{\mathrm{ref}}(z)}{Q^{\mathrm{ref}}(z)}\right)}^2⟩}_z}$$

$$D{R}_{\mathrm{loss}}(\alpha ,s^{+})=\frac{D{R}_{\max }-DR(\alpha ,s^{+})}{D{R}_{\max }},$$

$$\frac{h_r(y)}{h}=1-\left|\frac{y-y_s}{{\Lambda }_f/2}\right|,\text{for }-\frac{1}{2}\le \frac{y-y_s}{{\Lambda }_f}\le \frac{1}{2}.$$

$${\Gamma }_s$$

$${\Gamma }_m$$

$${\Gamma }_t$$

$$0=$$

$$0=\underset{\text{adv z}}{\underbrace{-{⟨\frac{\partial \overline{uw}}{\partial z}⟩}_{xy}}}+\underset{\text{dif z}}{\underbrace{{⟨\frac{1}{R{e}_b}\frac{{\partial }^2\overline{u}}{\partial z^2}⟩}_{xy}}}+\underset{\text{driving}}{\underbrace{{⟨{\gamma }_f\overline{f^{\mathrm{driving}}}⟩}_{xy}}}+\underset{\text{IBM}}{\underbrace{{⟨\overline{f_x^{\mathrm{IBM}}}⟩}_{xy}}}.$$

$$\frac{1}{12}\overline{f^{\mathrm{driving}}}=\frac{1}{R{e}_b}+{\int }_0^1\left(\frac{1}{2}-z\right){⟨-\overline{u^{\prime }{w}^{\prime }}⟩}_{xy}\mathrm{d}z,$$

$$d^{\mathrm{tot}}$$

$$={⟨\frac{1}{R{e}_b}{\frac{\partial \overline{u}}{\partial z}|}_{z=0}⟩}_{xy}-{⟨\frac{1}{R{e}_b}{\frac{\partial \overline{u}}{\partial z}|}_{z=1}⟩}_{xy}-{⟨\overline{f_x^{\mathrm{IBM}}}⟩}_V.$$

$$0=\sum _l{\int }_0^1{\int }_0^z{\int }_0^{\tilde{z}}\frac{1}{2}\left({f^{(l)}|}_{\hat{z}}+{f^{(l)}|}_{1-\hat{z}}\right)\mathrm{d}\hat{z}\mathrm{d}\tilde{z}\mathrm{d}z.$$

$$\overline{f^{\mathrm{driving}}}\left[{\int }_0^1\left\{\frac{1}{2}z\left(1-z\right)\right\}{⟨{\gamma }_f⟩}_{xy}\mathrm{d}z\right]$$

$$\frac{1}{R{e}_b}+{\int }_0^1\left\{\frac{1}{2}-z\right\}{⟨-\overline{uw}⟩}_{xy}\mathrm{d}z$$

$$d^{\mathrm{tot}}$$

$$=\sum _l{d}^{(l)}.$$

$$DC$$

$$=\sum _l\left\{\frac{d^{(l)}-d_{\mathrm{smooth}}^{(l)}}{d_{\mathrm{smooth}}^{\mathrm{tot}}}\right\}\equiv \sum _{l}d{c}^{(l)}.$$

$$d^{(l)}$$

$$={⟨{\int }_0^{1}g{⟨\overline{f}⟩}_x\mathrm{d}z⟩}_y\equiv {⟨d_{sp}^{(l)}⟩}_y.$$

$$d{c}^{(l)}$$

$$={⟨\frac{d_{sp}^{(l)}-d_{\mathrm{smooth}}^{(l)}}{d_{\mathrm{smooth}}^{\mathrm{tot}}}⟩}_y\equiv {⟨d{c}_{sp}^{(l)}⟩}_y,$$

$$\frac{\partial u_i}{\partial t}$$

$$u_i^{(s)}+{\alpha }_{(s)}\Delta t\frac{\partial {\tilde{p}}^{(s)}}{\partial x_i}$$

$${\alpha }_{(s)}\Delta t{\gamma }_f{f}^{\mathrm{driving}(s)}{\delta }_{i1}$$

$$\equiv u_i^{\ast (s)}+{\alpha }_{(s)}\Delta t{\gamma }_f{f}^{\mathrm{driving}(s)}{\delta }_{i1}$$

$$\equiv u_i^{\ast \ast (s)},$$

$${⟨u^{(s)}⟩}_V+{\alpha }_{(s)}\Delta t{⟨\frac{\partial {\tilde{p}}^{(s)}}{\partial x}⟩}_V={⟨u^{\ast (s)}⟩}_V+{\alpha }_{(s)}\Delta t{⟨{\gamma }_f⟩}_V{f}^{\mathrm{driving}(s)}.$$

$$f^{\mathrm{driving}(s)}=\frac{1-{⟨u^{\ast (s)}⟩}_V}{{\alpha }_{(s)}\Delta t{⟨{\gamma }_f⟩}_V}.$$

Full text

Drag Reduction by Herringbone Riblet Texture in Direct Numerical Simulations of Turbulent Channel Flow

H.O.G. Benschop*∗*   and W.P. Breugem

*Laboratory for Aero and Hydrodynamics, Delft University of Technology, Delft,

The Netherlands

∗Corresponding author. Email: H.O.G.Benschop@tudelft.nl

A bird-feather-inspired herringbone riblet texture was investigated for turbulent drag reduction. The texture consists of blade riblets in a converging/diverging or herringbone pattern with spanwise wavelength $\Lambda_{f}$. The aim is to quantify the drag change for this texture as compared to a smooth wall and to study the underlying mechanisms. To that purpose, Direct Numerical Simulations of turbulent flow in a channel with height $L_{z}$ were performed. The FIK-identity for drag decomposition was extended to textured walls and was used to study the drag change mechanisms. For $\Lambda_{f}/L_{z}\gtrsim O(10)$, the herringbone texture behaves similarly to a conventional parallel-riblet texture in yaw: the suppression of turbulent advective transport results in a slight drag reduction of 2%. For $\Lambda_{f}/L_{z}\lesssim O(1)$, the drag increases strongly with a maximum of 73%. This is attributed to enhanced mean and turbulent advection, which results from the strong secondary flow that forms over regions of riblet convergence/divergence. Hence, the employment of convergent/divergent riblets in the texture seems to be detrimental to turbulent drag reduction.

Keywords: Drag Reduction; Riblets; Direct Numerical Simulations

1 Introduction

Drag-reducing techniques can be of great value for fuel consumption reduction, as a significant part of the fuel used for transportation arises from drag in turbulent flows. Drag reduction (DR) in fluid flows can be obtained by active and passive methods. Active methods include the use of additives [90], gas injection [9], wall mass transfer [43], wall cooling or heating [92, 44, 78], wall motion [62], wall deformation [77] and electromagnetic forcing [72]. Passive methods include the change of surface chemistry (e.g. superhydrophobicity [66]), elasticity (e.g. compliant walls [16]), shape (e.g. airfoil or ship hull shape) and texture (e.g. riblets [19]). Surface roughness and surface textures have been successfully used for drag reduction by transition delay in laminar flow [26], separation delay in turbulent flow over curved surfaces [15, 74] and turbulence modification in turbulent flow over flat surfaces.

Of the investigated flat plate textures - such as sifted sand grains [2], transverse square grooves [86], dimples [1] and V-shaped protrusions [67] - the well-studied riblet texture seems most promising for turbulent DR. This texture has been found on the scales of some fast-shark skins [20] and consists of ridges or riblets aligned with the mean flow direction. With a simplified riblet geometry, a maximum DR of almost 10% has been obtained [5]. In the search for even higher values of DR, many variations on the standard riblet geometries have been investigated [19], such as hierarchical or compound riblets [91], riblets on a spanwise traveling surface wave [47], oscillating riblets [88, 37, 83], riblets in a wave-like pattern (either in phase [36] or out of phase [69]) and riblets combined with drag-reducing polymers [12]. The rationale behind these alternatives is to further reduce drag by somehow incorporating other drag-reducing methods, such as oscillating walls or polymer addition.

A little-studied alternative to the standard riblets is the herringbone riblet texture that has been found on bird flight feathers. Feathers serve several functions, including flight, thermal insulation, waterproofing and coloration (e.g. for camouflage or visual signals). The following description of the feather morphology is derived from Chen et al. [10, 11, 13]. A feather consists of two vanes (anterior and posterior), separated by a central supporting shaft. Each vane consists of parallel barbs attached obliquely to the shaft. The barbs are linked together by a set of finer barbs, called barbules. A microgroove is formed between neighbouring barbs. The angle between shaft and barbs is typically $\alpha\approx 30^{\circ}$. The groove spacing $s$ remains approximately constant with $s^{+}\approx 20$. The groove depth decreases gradually away from the shaft. A wing is formed when several feathers are positioned next to each other with almost parallel shafts. The feather shafts are approximately parallel to the flying direction for steady forward glide.

The study of Chen et al. [11] is one of the first to investigate the herringbone riblet texture for turbulent drag reduction. Inspired by bird flight feathers, they designed a riblet texture with two typical features that differ from the shark-skin riblets. First, the riblets were arranged in a converging/diverging or herringbone pattern. Second, the riblet height or groove depth changed gradually. Using laser engraving and replica molding, the researchers manufactured such bio-inspired herringbone-riblet skins. They covered the inner wall of a test pipe with these skins and measured a DR of up to 20%, twice that of optimal standard riblets. Although this seems promising, it has not been reproduced yet, and it is unclear for which texture and flow parameters drag is maximally reduced.

Sagong et al. [67] investigated a rather comparable geometry, namely the V-shaped protrusions that were found on the sailfish skin. In a comprehensive experimental and numerical study, they found a few cases for which drag was decreased slightly ($\sim 1$%), although within the experimental uncertainty.

The herringbone riblet texture also resembles the vortex generators that have been used to reduce drag by flow-separation delay [49] or transition delay [71]. Furthermore, several studies confirm that roughness on a bird wing contributes to separation control [8, 48, 79].

Herringbone riblet textures have been studied for several other reasons, such as for heat-transfer enhancement [29, 24] and mixing of laminar flows in microchannels [76]. The hot-wire study by Koeltzsch et al. [45] was initiated by the observation of convergent and divergent riblet patterns on the shark skin near sensory organs, possibly used for local flow-noise reduction. Nugroho et al. [56] revisited the same texture and conducted a parametric study using the hot-wire technique. They were motivated by the potential use for passive flow control and disruption of large-scale coherent motions. They proposed that the herringbone riblets induce large-scale counter-rotating vortices, giving rise to large-scale spanwise periodicity in the boundary layer. The suggested streamwise counter-rotating vortices have recently been visualized in a flat-plate laminar boundary layer over convergent riblets [54]. So, the motivation for herringbone riblet studies has been drag reduction, heat-transfer enhancement, mixing improvement, flow control and large-scale flow structures manipulation.

The study by Nugroho et al. [56] is a manifestation of the recent interest in surfaces that manipulate the whole boundary layer. There is much evidence that DR techniques that rely on near-wall flow manipulation are less effective at higher Reynolds numbers: $DR$ degrades with increasing $Re$ [40, 41, 75, 33]. That has partially motivated the study of rough or textured surfaces that also impact the outer part of the turbulent boundary layer by means of large-scale secondary flows. These secondary flows have been observed over spanwise-varying longitudinal bedforms [87], a regular array of cubes [65], a bed with two parallel lanes of different roughnesses [81], the irregular surface of a replica of a damaged turbine blade [4], streamwise strips of elevated roughness [80] and the aforementioned herringbone riblet texture [45, 56]. Interestingly, Schoppa and Hussain [70] obtained 20% DR in Direct Numerical Simulations of turbulent channel flow with imposed large-scale counter-rotating streamwise vortices.

The aim of this paper is to quantify how and why drag is changed by a herringbone texture as compared to a smooth wall, using Direct Numerical Simulations (DNSs). The numerical methods are described in section 2. Section 3 explains how drag reduction is quantified. Section 4 validates the simulated drag of smooth walls, parallel riblets aligned with the mean flow, and parallel riblets in yaw. Section 5 examines the herringbone texture with use of a parametric study. A drag decomposition is derived and applied in section 6. A discussion of the results is presented in section 7, followed by the main conclusions and an outlook in section 8.

2 Numerical Methods

In this study, Direct Numerical Simulations (DNSs) of incompressible turbulent flow were performed. An overview of all simulations with the corresponding parameters can be found in Appendix C. This section outlines the method that was used, describing successively the notation conventions, flow domain, texture, grid, numerical code and averaging.

Notation In what follows, dimensional variables are denoted by an asterisk ∗. Variables without that asterisk are nondimensionalized using the domain height $L_{z}^{*}$ and the bulk velocity $U_{b}^{*}$, such that $L_{z}=1$ and $U_{b}=1$. Note that $U_{b}^{*}$ is a constant, as simulations were performed at fixed mass flow rate. The constant bulk Reynolds number is defined as $Re_{b}=U_{b}^{*}L_{z}^{*}/\nu^{*}$, with kinematic viscosity $\nu^{*}$. The superscript + is used for nondimensionalization with $\nu$ and $u_{\tau}=\sqrt{\tau_{w}/\rho}$, with wall shear stress $\tau_{w}$ and fluid density $\rho$. Nondimensionalization for textured walls uses the viscous wall units derived from the smooth-wall simulation with the same $Re_{b}$.

Domain The flow domain is a plane channel, bounded by two horizontal walls. At the channel walls, no-penetration and no-slip boundary conditions are applied, whereas periodic boundary conditions are used in the streamwise and spanwise directions. The streamwise, spanwise and wall-normal coordinates are denoted by $x$, $y$ and $z$ with the corresponding velocity components $u$, $v$ and $w$.

The domain is specified by its length $L_{x}$, width $L_{y}$ and height $L_{z}$. For a good comparison, the domain size should ideally be the same for all simulations. However, slight size variation was needed to fit an integer number of texture periods in the streamwise and spanwise directions, or to ensure that the number of grid cells complies with the parallel-computing algorithm. In general, all domains are approximately of size $(4.0\times 2.5\times 1)$, which is considered to be large enough to obtain reliable statistics for several reasons. First, it is comparable to the domain size used by other researchers [53, 7, 58, 85]. In addition, it is full-span, as opposed to the recently reconsidered minimal-span channels [17]. Finally, it is larger than the moderate box of size $(\pi\times\pi/2\times 1)$ that is large enough to reproduce the one-point statistics of larger boxes [50].

Texture Textures are applied to the inside of both channel walls to enforce symmetry in the mean flow [31]. Unless stated otherwise, the top wall texture is the bottom wall texture mirrored in the centerline plane. To simulate flow over a non-smooth surface, two methods can be adopted: coordinate transformation or the Immersed Boundary Method (IBM) [58]. Both methods have been applied to simulate turbulent flow over riblet walls [14, 35]. We used an IBM similar to the one employed by Breugem and Boersma [7] and Pourquie et al. [61], which is based on Fadlun et al. [23]. The IBM forcing is direct, i.e. a forcing term is added to the discretized equations. Appendix B provides case-specific details.

Grid The chosen IBM allows the use of a simple staggered Cartesian grid. The number of grid cells in the three Cartesian directions is denoted by $N_{x}$, $N_{y}$ and $N_{z}$. The grid is uniform in the horizontal directions, so the grid spacings $\Delta x$ and $\Delta y$ are constant. In the wall-normal direction, three zones are distinguished, namely the roughness regions near the two walls and the remaining part of the channel. In the roughness region, which extends from the wall to one grid cell above the maximum texture height, the vertical grid spacing $\Delta z_{w}$ is constant. In the remaining part of the channel, grid stretching is applied using a cosine function that is symmetrical with respect to the channel centerline. The maximum vertical grid spacing occurs at the centerline and is called $\Delta z_{c}$.

The grid-cell size is important to correctly resolve small-scale fluid motions. The recommendations of Vreman and Kuerten [84] for finite difference codes were followed, namely ${\Delta x}^{+}=4.4$, ${\Delta y}^{+}=2.9$, ${\Delta z_{w}}^{+}=0.49$ and ${\Delta z_{c}}^{+}=2.2$ as maximum grid spacings for smooth-wall turbulent flows. For textured walls, the spanwise grid spacing was reduced to ${\Delta y}^{+}\approx 1.0$. For all investigated textures, one simulation at a double spanwise and/or streamwise resolution was performed. Particular attention was paid to the resulting drag change, which was marginal in all cases. As only one simulation at a higher resolution was performed for each texture, grid independence of the results cannot be claimed. However, as the grid resolution is relatively high and about the same in all cases, comparison of results is still justified.

Code The incompressible Navier-Stokes equations and continuity equation were solved at fixed bulk velocity:

[TABLE]

where the Einstein summation convention for repeated indices is used. Here $u_{i}$ represents one component of the velocity vector, $t$ time, $x_{j}$ a spatial coordinate, $p$ the pressure and $f^{\mathrm{driving}}$ the spatially uniform forcing term to obtain a constant bulk velocity. The phase-indicator function $\gamma_{f}$ is defined at grid points of the streamwise velocity. It equals 1 in fluid and 0 in solid obstacle volume to ensure that only fluid experiences the bulk forcing that drives the flow. The Kronecker-delta function $\delta_{i1}$ guarantees that fluid is driven in the streamwise direction with $i=1$. Appendix A describes how $f^{\mathrm{driving}}$ is calculated. The IBM forcing $f^{\mathrm{IBM}}_{i}$ is a body force that models the (drag) force that the texture exerts on the flow (see Appendix B for more details).

These equations are discretized using the finite-volume method combined with a pressure-correction scheme [25]. Fluxes or stresses at the cell faces are evaluated using linear interpolation, i.e. a central-differencing scheme is used. Time-integration is performed using a fractional-step method that consists of three steps. For the pressure, which is staggered in time with respect to the velocities, a Crank-Nicolson scheme is used. All other terms are advanced in time using a three-step Runge-Kutta method [89]. This discretization procedure yields a Poisson equation, which is solved using a non-iterative FFT-based solver. Specifically, FFTs are applied to the horizontal directions and the resulting tridiagonal system is solved using Gaussian elimination. More details about the time advancement at fixed bulk velocity can be found in Appendix A.

Averaging For computation of flow statistics, a unit-cell average was stored each 100 timesteps. Like in crystallography, a unit cell is the smallest unit of volume that builds up the entire texture by translation. It extends vertically from bottom to top wall. For smooth walls, its size in grid cells is $1\times 1\times N_{z}$. In a unit-cell average, the data of all unit cells are reduced to an average in one unit cell.

In addition to this unit-cell average during the computations, temporal and spatial averages were performed afterwards. Let $\phi=\phi(x,y,z,t)$ represent an arbitrary flow variable. The following averages were used:

[TABLE]

Here $\overline{\phi}$ is a time average over the statistically stationary part of the signal $\phi$ (which starts at $t_{\mathrm{stat}}$ and has duration $T$). An average over one spatial coordinate is denoted by $\left\langle\phi\right\rangle_{x_{i}}$. For instance, $\left\langle\phi\right\rangle_{x}$ is a streamwise average. A similar notation is adopted for an average over two spatial coordinates. For example, $\left\langle\phi\right\rangle_{xy}$ is a streamwise and spanwise average. Finally, $\left\langle\phi\right\rangle_{V}$ represents an average over the entire volume $V=L_{x}L_{y}L_{z}$. For instance, the bulk velocity is defined by $U_{b}=\left\langle u\right\rangle_{V}$.

3 Drag Reduction Quantification

This section describes the quantification of drag reduction. It explains under what conditions the smooth- and textured-wall flows are compared, with special attention for flow generation, the definition of a reference case and the formulation of the benefit of drag reduction.

In general, channel flow is generated by either a constant flow rate (CFR), a constant pressure gradient (CPG) or a constant power input (CPI) [27]. The obtained drag reduction depends slightly on the choice of CFR, CPG or CPI. In this study, the flow was generated by CFR.

The quantification of drag reduction requires the definition of a reference case compared to which drag is reduced or increased. Special attention should be paid to the Reynolds number and channel height, as drag depends heavily on both of them. Conceptually, this study compares two channels with the same fluid, the same flow rate per unit of spanwise width and the same outer dimension $L_{z}^{*}$. With $U_{b}^{*}L_{z}^{*}$ being the volumetric flow rate per unit width, the first two conditions imply that $Re_{b}$ is the same for both flows. The third condition guarantees that both channels are geometrically identical, apart from the texture that is applied to the inside of the channel walls in one case. It implies that the wall location is not adjusted to compensate for the texture volume. The fluid volume for the case with textured walls is thus slightly less than that of the smooth-wall case. This is a conservative choice: drag reduction cannot result from an increased fluid volume or a locally increased channel height [18, 52].

This paper uses the drag change $DC$ as a quantifying parameter. To account for possible differences in domain width and length, $DC$ was computed from the time-averaged drag force per unit volume. As the driving term balances the total drag, the instantaneous drag force per unit volume $f_{d}^{*}=\rho^{*}\left(U_{b}^{*}\right)^{2}\left(L_{z}^{*}\right)^{-1}V^{-1}\int_{V}\gamma_{f}f^{\mathrm{driving}}\mathrm{d}V$. Since the dimensional prefactor is equal for the smooth- and textured-wall channel flows, the drag change is given by

[TABLE]

It measures the increase of the driving force that is required to maintain a given flow rate. The drag reduction $DR=-DC$, so drag is reduced in case $DC$ is negative.

The thus computed drag change is supplemented by a 95% confidence interval. The error in $DC$ can be attributed to the uncertainty in the drag computed for both the textured and smooth wall. For ease of notation, define $d_{t}\equiv d^{\mathrm{tot}}_{\mathrm{textured}}$ and $d_{s}\equiv d^{\mathrm{tot}}_{\mathrm{smooth}}$. Let $u_{\phi}$ for now denote the uncertainty in $\phi$. Given the independence of $d_{t}$ and $d_{s}$, the uncertainty in $DC$ follows from the law of error propagation:

[TABLE]

The uncertainties in $d_{t}$ and $d_{s}$ were computed using the method outlined by Hoyas and Jiménez [39]. It accounts for correlation in the drag time signal. The thus obtained error bar only results from the finite simulation time. Errors of other origins (e.g. discretization errors) were not considered.

4 Validation

The numerical methods were validated with simulations of smooth walls, parallel riblets, and parallel riblets in yaw, as described in the following subsections.

4.1 Smooth Wall

Smooth-wall DNSs were performed at three bulk Reynolds numbers, namely 5500, 11000 and 22000. Table 1 shows the corresponding friction Reynolds number $Re_{\tau}=u_{\tau}\delta/\nu$, where $\delta=L_{z}/2$ is the half-channel height. Compared to recent DNSs reaching $Re_{\tau}=4000$ [6], $Re_{\tau}\approx 4200$ [50] and $Re_{\tau}\approx 5200$ [46], the simulations in this study are considered to be standard. Therefore, this subsection suffices to validate drag and flow statistics.

Drag is validated by comparing $Re_{\tau}$ with the value predicted by an approximate relation $Re_{\tau}^{\mathrm{Pope}}=0.09Re_{b}^{0.88}$ [60]. The deviation of $Re_{\tau}$ from $Re_{\tau}^{\mathrm{Pope}}$ is quantified using

[TABLE]

Table 1 lists $\Delta Re_{\tau}^{\mathrm{Pope}}$ for all smooth-wall simulations. Although the relation for $Re_{\tau}^{\mathrm{Pope}}$ is approximate, good agreement is obtained for all $Re_{b}$.

Flow statistics of $U$, $u_{\mathrm{rms}}$, $v_{\mathrm{rms}}$, $w_{\mathrm{rms}}$ and $p_{\mathrm{rms}}$ were compared with simulations at $Re_{\tau}=180$ [84] and $Re_{\tau}=590$ [85]. Here $U=\left\langle\overline{u}\right\rangle_{xy}$, $u_{\mathrm{rms}}=(\langle\,\overline{u^{2}}\,\rangle_{xy}-\left\langle\overline{u}\right\rangle_{xy}^{2})^{1/2}$ and similarly for the other root-mean-square quantities. Let $Q(z)$ be one of these statistics, then the root-mean-square relative deviation

[TABLE]

is used to quantify the difference between $Q$ (current) and $Q^{\mathrm{ref}}$ (reference). It was computed after piecewise cubic spline interpolation of $Q$ and $Q^{\mathrm{ref}}$ to a uniform grid, $z^{+}(k)=k$ for integers $1\leq k\leq k_{\max}$, with $k_{\max}=175$ at the lowest and $k_{\max}=587$ at the highest $Re_{b}$. Table 1 shows $\delta Q$ for the five flow quantities. The root-mean-square relative deviation is smaller than 1% for the mean velocity and smaller than 2% for the root-mean-square fluctuations. The somewhat larger value for $p_{\mathrm{rms}}$ at the lowest $Re_{b}$ is attributed to the difference between $Re_{\tau}$ and $Re_{\tau}^{\mathrm{ref}}$.

4.2 Parallel Riblets

A parallel riblet texture consists of riblets aligned with the mean flow direction and can reduce turbulent drag up to almost 10% [5]. An overview of previous research can be found elsewhere [19]. A thorough DNS study has been performed quite recently [30, 31, 32].

This paper investigates the blade riblet texture (see Figure 1). The blades have zero thickness, spacing $s$ and height $h$ with $h/s=0.5$. In a small parametric study, mainly $s^{+}$ and $Re_{b}$ were varied (see also Table 3). The grid resolution (specified in wall units) is about the same for all cases.

Figure 2 (left) compares the simulated drag change as function of $s^{+}$ with experiments performed by Bechert et al. [5]. The top axis shows the experimental bulk Reynolds number $Re_{\mathrm{Bechert}}$, which is based on the horizontal channel width and the average velocity between the test plates. In the experiments only the Reynolds number was varied, whereas in the present numerical study both the Reynolds number and the riblet spacing were varied.

The drag change varies slightly with Reynolds number for fixed $s^{+}$: the drag at $Re_{b}=5500$ is higher than at $Re_{b}=11000$ and $22000$. The approximate overlap of the data points at the two highest Reynolds numbers (for $s^{+}=24$) suggests a low-Reynolds-number effect, which is underpinned by the observation that drag reduction data below $Re_{\mathrm{Bechert}}\approx 10000$ deviated more and more from previous high-Reynolds-number data [5]. The deviation of $DC$ at $Re_{b}=5500$ from that at higher Reynolds numbers is also larger at larger $s^{+}$, which might be explained by riblet height increase. For $Re_{b}=5500$, an increase of $s^{+}$ from 10 to 24 is accompanied by a decrease of $\delta/h$ from 35 to 15. The blades protrude farther into the channel, which is presumably detrimental to drag reduction. This is supported by the finding that, for $\delta/h\lesssim 50$, the effect of roughness extends across the boundary layer and the original wall flow dynamics is changed significantly [42].

The simulations capture the experimental trend quite well when $Re_{b}$ is close to $Re_{\mathrm{Bechert}}$. Especially good agreement is obtained around and below the optimum spacing. At $s^{+}=17$, a maximum drag reduction of 9.3% is achieved, which is very close to the 9.9% of the experiments. The difference is slightly larger at $s^{+}=24$: there the total drag is reproduced within 5%. A higher reproduction accuracy probably requires a combined experimental/numerical study, a more extensive grid resolution study, the incorporation of blade thickness and an analytical correction of momentum fluxes near the riblet tips to resolve the high gradients there. However, the accuracy demonstrated here is sufficient for the herringbone riblet simulations, as drag differences for these were found to be much larger than 5%.

4.3 Parallel Riblets in Yaw

Among other factors, the performance of riblets deteriorates in yaw, i.e. when they are not aligned with the mean flow direction. An overview of past research is given by Koeltzsch et al. [45]. The study of Hage et al. [38] shows the significant influence of riblet geometry and spacing $s^{+}$. Simulations of turbulent flow over parallel riblet textures in yaw are rare in the literature. Grüneberger et al. [36] have performed DNSs at $Re_{b}=5750$ ($Re_{\tau}=180$) for trapezoidal grooves. The driving pressure gradient was rotated such that the streamwise direction no longer coincides with the $x$-direction. Good agreement with experimental data was obtained.

This paper investigates the blade riblet texture in yaw (see Figure 1) with $h/s=0.5$ and $s^{+}=17$. The yaw angle $\alpha$ was varied for fixed $Re_{b}=5500$. The horizontal grid resolution was comparable to that used for aligned parallel riblets: ${\Delta x}^{+}\leq 4.1$ and ${\Delta y}^{+}\leq 1.1$ for all cases. As the trick of driving-pressure-gradient rotation cannot be applied for herringbone-riblets simulations, the parallel riblet texture was rotated with respect to the grid. That required a different Immersed Boundary Method (IBM), as the blades are not anymore aligned with the Cartesian directions (see Appendix B for details).

Figure 2 (right) shows $DC$ as function of $\alpha$ (indicated by the numbers in the figure). Clearly, the drag reduction diminishes with increasing yaw angle, as is expected from the literature. Unfortunately, to the best of the authors’ knowledge, no experimental data of blade riblets in yaw is available. For comparison of the simulated results with experimental data for other riblet geometries, the drag reduction loss $DR_{\mathrm{loss}}$ is introduced [38]:

[TABLE]

where $DR=-DC$, and $DR_{\mathrm{max}}$ is the maximum drag reduction that can be obtained with a given geometry: $DR_{\mathrm{max}}\!=DR(\alpha\!=\!0,s^{+}\!=\!s^{+}_{\mathrm{opt}})$. It was assumed that $DR_{\mathrm{max}}\!=DR(\alpha\!=\!0,s^{+}\!=\!17)$ for the simulations.

Figure 2 shows the drag reduction losses as function of yaw angle for different riblet geometries at $s^{+}=17$. The experimental data of Hage et al. [38] (at $Re_{b}\approx 14900$) and the numerical data of Grüneberger et al. [36] (at $Re_{b}=5750$) are included. The values of $h/s$ are geometry dependent, namely $h/s=0.5$ for the trapezoidal and blade, $h/s=0.7$ for the semicircular, and $h/s=1$ for the triangular geometries.

The dependence of $DR_{\mathrm{loss}}$ on $\alpha$ is different for each geometry, which might be explained by differences in $h/s$ and riblet shape [38]. Of the two trapezoidal geometries, the one with the sharpest tip ($30^{\circ}$) is slightly more sensitive to misalignment. The triangular geometry has the largest tip angle of about $54^{\circ}$, but its drag reduction loss is nevertheless larger than for the trapezoidal geometries, presumably because of its larger height ($h/s=1$). Of the experimental data, the semicircular geometry is most sensitive to yaw, likely because of its larger height ($h/s=0.7$) and sharp tip.

The simulation results suggest that the blade riblet geometry is more susceptible to yaw than any of the other geometries. This is likely not a low-Reynolds-number effect, as the simulations at $Re_{b}=5750$ by Grüneberger et al. [36] reproduced well the experimental data (see the figure). Instead, it might be ascribed to the sharper riblet tips (and the associated pressure drag increase) and the broader riblet valleys (and the associated increase of sloshing) [38].

5 Herringbone Riblets

5.1 Texture Description

The feather texture (described in section 1) is modeled as shown in Figure 3, fairly similar to the spatial 3D (s-3D) texture proposed in Chen et al. [11]. It consists of $N_{\mathrm{feather}}$ feathers placed in parallel such that their edges touch. The modeled feather has no physical shaft, although the term ‘shaft’ is employed to denote the symmetry axis of a feather. The barbs are modeled as blade riblets with zero thickness. As opposed to the s-3D texture that consists of sawtooth riblets, the present study uses blades because of their superior drag-reducing performance in the conventional riblet texture. The angle between the positive $x$-direction and the blades is called $\alpha$, which is restricted to $0\leq\alpha<180^{\circ}$. The shortest distance between the blades is $s$. Between two neighboring blades a groove forms, which has spanwise blade spacing $s_{y}=s/\cos(\alpha)$. $N_{\mathrm{groove}}$ is the integer number of spanwise blade spacings $s_{y}$ that fits in one feather half-width (see also Figure 3). The feather width or spanwise texture wavelength $\Lambda_{f}=2N_{\mathrm{groove}}s_{y}$. The riblet height $h_{r}$ decreases linearly with distance to the shaft. Let $h$ represent the riblet height at the shaft and $y_{s}$ the $y$-coordinate of the shaft, then $h_{r}$ is given by

[TABLE]

Note that Koeltzsch et al. [45] and Nugroho et al. [56] used herringbone riblets of constant height. The texture was implemented using the same IBM that was used for blade riblets in yaw; only the texture indicator functions were different.

As Figure 3 shows, a difference is made between forward and backward flow. ‘Forward’ is used for bulk flow in the positive, ‘backward’ indicates flow in the negative $x$-direction. Forward flow over a texture with angle $\alpha$ is the same as backward flow over a texture with angle $\alpha_{\mathrm{backward}}=180^{\circ}-\alpha$. Using this trick, textures with angle $\alpha$ and $\alpha_{\mathrm{backward}}$ can be compared to study the influence of mean flow direction on drag.

In addition to the standard herringbone texture, two texture variations were considered. The first variation is called ‘shifted’. It differs from the standard herringbone texture by a spanwise shift of the top wall texture by half a spanwise texture wavelength as compared to the bottom wall (see Figure 4). This shifted texture was investigated for its ability to generate the drag-reducing secondary flow that is described by Schoppa and Hussain [70], namely one that extends from the bottom to the top wall. The second variation is a riblet texture with $\alpha=0^{\circ}$ (see Figure 4), which results in a parallel blade riblet geometry with blade height variation in the spanwise direction given by Equation 11. This texture does not suffer from yaw and the resulting pressure drag, but it still might give rise to secondary flows.

The herringbone texture and its variations have been subjected to a parametric study (see also Table 3). According to Chen et al. [11], bird-feather parameters are typically $s^{+}\approx 20$, $h/s\approx 0.5$ and $\alpha\approx 30^{\circ}$. Values for $N_{\mathrm{groove}}$ or $\Lambda_{f}$ were not given. Their SEM-pictures show that the feather barbs are comparable to blades with finite thickness and rounded tips. Although the present study was inspired by the bird-feather texture, it did not attempt to exactly reproduce that texture. Instead, texture parameters were chosen to allow a close comparison with the parallel-riblet studies described in subsections 4.2 and 4.3. Specifically, ridge spacing and maximum blade height were fixed at $s^{+}=17$ and $h/s=0.5$. Mainly two parameters were varied, namely $\alpha$ and $N_{\mathrm{groove}}$. The angle $\alpha$ was either $15^{\circ}$ or $165^{\circ}$, which permits the direct comparison with the conventional riblet texture in 15 degrees yaw. Larger angles (e.g. $\alpha=30^{\circ}$ or $\alpha=150^{\circ}$) were not considered, as Figure 2 suggests that such angles will result in a drag increase. The number of grooves $N_{\mathrm{groove}}$ was varied between 1 and 128, which resulted in a feather-width change from $\Lambda_{f}=0.10$ to $12.9$ in outer units or $\Lambda_{f}^{+}=35$ to $4506$ in wall units.

Apart from texture parameters, the other simulation parameters were (almost) the same for all herringbone simulations, which facilitates a fair comparison. Specifically, Table 3 shows that domain size, grid resolution, Reynolds number and simulation time do not change much among the simulations. In addition, these parameters are also close to the ones for parallel-riblet (without/in yaw) simulations. All herringbone simulations were performed at $Re_{b}=5500$. Based on the validation of blade-riblet simulations (see subsection 4.2), a low-Reynolds-number effect can be anticipated. However, as all simulations were performed at the same $Re_{b}$, a proper comparison can still be made. The number of grid cells per groove $N_{cg}=16$, so 16 grid cells span the groove width.

To substantiate the accuracy of the used numerical methods, two issues have to be addressed. The first issue relates to the boundary conditions at the texture surface. Appendix B.4 shows that the no-slip and no-penetration conditions are sufficiently satisfied. The second issue relates to the flow around the blade tip. The exact solution is singular there, which might introduce errors in the numerical solution. The improvement of the numerical accuracy for problems involving singularities is a research in itself [73], but it is not the focus of the current study. For the present purpose, the used IBM is sufficiently accurate. Very similar IBMs have been used in DNSs of flow around other obstacles with sharp corners, such as cubes [7, 57], square and triangular elements [58], and a flat plate normal to the free stream [68, 55]. The penultimate example shows that the IBM can deal with obstacles that are not aligned with the Cartesian grid, whereas the last example demonstrates that the IBM can also accurately capture separating flows at sharp corners.

5.2 Drag Reduction

Trend Figure 5 shows the drag change as function of feather width for the herringbone riblet geometry. Drag changes significantly with the spanwise texture wavelength or feather width. For the texture with $\alpha=165^{\circ}$, the drag increases by 61% for very narrow feathers, while the drag decreases by 2% for very wide feathers. In general, wider feathers experience less drag. Only the texture with $\alpha=15^{\circ}$ and $N_{\mathrm{groove}}=1$ does not comply with this trend. Drag reduction was only obtained when the texture approaches the parallel riblet texture in yaw, i.e. in the limit of very wide feathers. Only 2% reduction was found in that limit, which is indeed comparable to the 0.9% reduction that was obtained for parallel riblets at yaw angle $\alpha=15^{\circ}$. It is significantly less than the 7.6% reduction for parallel riblets aligned with the mean flow.

Difference forward and backward flow There is a clear drag difference between textures with $\alpha=15^{\circ}$ (forward flow) and $\alpha=165^{\circ}$ (backward flow). For example, at $N_{\mathrm{groove}}=16$ drag increases by 34% for $\alpha=15^{\circ}$, and only by 7.5% for $\alpha=165^{\circ}$. This shows that mean-flow reversal has a considerable effect on drag, in agreement with Chen et al. [11]. Here, forward flow experiences more drag than backward flow, except for the $N_{\mathrm{groove}}=1$ texture.

The drag difference between forward and backward flow decreases for increasingly wide feathers, as can be expected. The textures with $\alpha=15^{\circ}$ and $165^{\circ}$ differ only because of riblet convergence or divergence in the feather shaft or feather edge regions. When $\Lambda_{f}\rightarrow\infty$, only a very small portion of the complete texture consists of converging or diverging riblets, so their contribution to the total drag becomes negligible.

Standard texture variations Figure 5 also shows the drag for the shifted and parallel riblet ($\alpha=0^{\circ}$) variations. The drag of textures with and without shift is about the same for all feather widths. Hence, shifting of the top wall texture has almost no effect on drag, although the next subsection will explain that the mean flow is different for some cases.

In contrast with most herringbone textures, the parallel riblet texture with spanwise riblet height variation is able to reduce drag. The drag reduction is the least for $N_{\mathrm{groove}}=1$, namely 0.5%. The texture with $N_{\mathrm{groove}}=1$ and $\alpha=0^{\circ}$ is the conventional parallel-riblet geometry with $s^{+}=17$, but with every second blade removed. It is known that this geometry is not optimal for drag reduction. When $N_{\mathrm{groove}}$ increases, drag reduction increases as well. In the limit of very large $N_{\mathrm{groove}}$, the texture approaches the standard parallel-riblet geometry, apart from a very slight spanwise height variation. Therefore, in that limit one might expect the drag reduction to be close to that for standard riblets.

5.3 Flow Description

Instantaneous streamwise velocity To understand the drag-reduction results presented in the previous subsection, a detailed flow analysis is indispensable. Figure 6 shows three cross sections of the instantaneous streamwise velocity. The V-shaped contours in the horizontal cross section reveal the presence of the herringbone texture. The plumes that appear in the $yz$-plane represent up- or downdrafts that result from the converging/diverging riblets in the texture.

Figure 7 shows part of the mean streamwise-averaged flow fields for textures with $\alpha=165^{\circ}$ and $N_{\mathrm{groove}}=1,4,16$. Because of flow symmetry, the figures show only one feather half. The feather shaft is located at the left side and the feather edge at the right side of the figures.

$\boldsymbol{\alpha=165^{\circ}}$ For textures with $\alpha=165^{\circ}$, the flow near the shaft converges, a local updraft of fluid results, and an approximately square vortex appears. The vortex does not reach the channel centerline, probably because of the small misalignment of the riblets with the mean flow ($\alpha=165^{\circ}$ is relatively close to $180^{\circ}$). For $N_{\mathrm{groove}}\geq 4$, it is limited in vertical extent to approximately $z=0.2$, as can be seen for $N_{\mathrm{groove}}=4$ and $16$ in the figure. Because of spanwise confinement, for small feather widths the vortex decreases in size. For $N_{\mathrm{groove}}=1$, it only reaches $z=0.05$. Secondary flow is thus suppressed for small spanwise roughness spacings. This accords with the experimental finding that secondary flow disappears for decreasing spanwise spacing of roughness elements [80].

While mainly one vortex constitutes the secondary flow for $N_{\mathrm{groove}}=1$ and $4$, tertiary flows appear for larger feather widths. For example, the mean flow for $N_{\mathrm{groove}}=16$ (see Figure 7) shows a counterclockwise-rotating flow that extends to the channel centerline. Less pronounced is the small clockwise-rotating vortex near the feather edge. Several tertiary flows were also observed for even wider feathers ($N_{\mathrm{groove}}=32,128$). This agrees with the experimental finding that tertiary flows appear when the spanwise spacing of roughness elements increases above the boundary layer thickness [80]. These tertiary flows are likely similar to the secondary flows that form over streamwise-aligned roughness strips [3]. The latter are stress-induced (i.e. Prandtl’s secondary flows of the second kind), as opposed to the flow-curvature-induced secondary flows (i.e. Prandtl’s secondary flows of the first kind).

$\boldsymbol{\alpha=15^{\circ}}$ Reversing the flow direction from backward to forward also changes the secondary flow direction. For textures with $\alpha=15^{\circ}$, riblets near the shaft diverge, resulting in a local downdraft of fluid and the appearance of a counterclockwise-rotating secondary flow to the right of the shaft. A tertiary flow again appears for $N_{\mathrm{groove}}\geq 16$. However, this time it is less well-defined, as its rotation direction is also counterclockwise. So, the most dominant tertiary flow does not change rotation direction by flow reversal.

Shifted variant Shifting of the top-wall textures changes the mean-flow patterns only slightly. For all cases, the dominant vortex near the top wall shifts in the spanwise direction. As this vortex does not reach the channel centerline, it has no noticeable effect on the mean flow in the bottom channel half. For $N_{\mathrm{groove}}=1$ and $4$, this spanwise shift is the only change of the mean flow. For $N_{\mathrm{groove}}=16$, however, the pronounced counterclockwise-rotating tertiary flow is also modified. It reaches the centerline for the not-shifted texture (see Figure 7), but it occupies the whole channel for the shifted variant. It is almost square, as for $N_{\mathrm{groove}}=16$ the feather half-width ($\Lambda_{f}/2=0.80$) approximately equals the channel height ($L_{z}=1$). Tertiary flow for $N_{\mathrm{groove}}=32$ also extends from bottom to top wall. However, these modified tertiary flows apparently have not much influence on drag (considering Figure 5).

Spanwise modulation The strong secondary-flow vortex near the shaft causes a spanwise modulation of the boundary layer, as is clear from the streamwise-velocity contours in Figure 7. Regions of updrafts (downdrafts) are characterized by low (high) streamwise velocity. The trend for streamwise velocity fluctuations is opposite, as the top-right subfigure shows: updrafts (downdrafts) are associated with increased (decreased) fluctuations. The same was found in experimental studies of convergent/divergent riblets [45, 56] and in DNSs of turbulent boundary layers with uniform blowing/suction [43].

Instantaneous vortical structures and wall shear Figure 8 shows the instantaneous vortical structures and wall shear for three cases, namely smooth walls, parallel riblets in yaw and herringbone riblets. The vortical structures are iso-surfaces of the second invariant of the velocity gradient tensor (the $Q$-criterion, see e.g. [22]). The contours represent the shear stress $\tau_{w}^{+}$ on the bottom wall.

The smooth-wall plot shows some well-known features, such as low-speed streaks and hairpin-type vortical structures. The streaks are much less apparent in the second plot, which is attributed to the parallel riblets that adjust the turbulence. The streaks are more evident in the flow above the riblets. The vortices are comparable to the ones for the smooth wall. The plot for parallel riblets without yaw is not shown, as it is very similar to the one for parallel riblets in yaw.

The herringbone riblet texture exhibits the largest changes in vortical structures and shear. High shear is associated with diverging riblets, and low shear with converging riblets. The vortices seem to be ordered as well: they are abundant over regions with updrafts, but almost absent over regions with downdrafts. The same trend was again found for blowing/suction: vortices are enhanced by blowing in spite of the reduced wall shear stress, while vortices are suppressed by suction despite the increase of wall shear stress [43].

5.4 Secondary Flow Strength

To quantify the strengths of the secondary flows described in the previous subsection, the secondary flow strength $\Gamma_{s}$ is introduced:

[TABLE]

with $v^{\prime}=v-\overline{v}$ and $w^{\prime}=w-\overline{w}$. The secondary flow strength is decomposed into the mean-secondary-flow strength $\Gamma_{m}$ and the turbulent-secondary-flow strength $\Gamma_{t}$. These variables still depend on the spatial coordinates. A volume- and $xz$-average of $\Gamma_{m}$ and $\Gamma_{t}$ is shown in Figure 9.

Influence of feather width Figure 9 (left part) shows that both mean and turbulent secondary flow become stronger with decreasing feather width. This is attributed to the converging or diverging riblets that trigger updrafts and downdrafts. When the feather width decreases, the spanwise density of these texture-generated secondary flows increases. In other words, the secondary flows that are created near the feather shaft and edges cover a relatively larger portion of the total fluid volume. This yields a stronger volume-averaged secondary flow for smaller feather widths.

Mean and turbulent contribution The mean- and turbulent-secondary-flow strengths follow the same trend for most textures (namely an increase with a decrease of $\Lambda_{f}$). For all $N_{\mathrm{groove}}$, the converging or diverging riblets near the shaft trigger updrafts and downdrafts. As Figure 6 demonstrates, these are not steady flow patterns. Instead, they can be considered as fluctuating ejections and sweeps that are generated by the texture. The resulting instantaneous secondary flow ($v^{2}+w^{2}$) projects both onto $\Gamma_{m}$ and $\Gamma_{t}$ (see Equation 12). A stronger instantaneous secondary flow yields in general an increase of both the mean and turbulent secondary flow. That explains why the mean and turbulent strength follow the same trend for $N_{\mathrm{groove}}\geq 4$. The results for $N_{\mathrm{groove}}=1$, however, deviate in this respect: $\left\langle\Gamma_{m}\right\rangle_{V}$ more than halves as compared to $N_{\mathrm{groove}}=4$. This suppression of mean secondary flow is due to spanwise confinement (see previous subsection). However, the fluctuating updrafts and downdrafts are not suppressed, so the turbulent secondary flow remains strong.

Change with spanwise distance Figure 9 (top right) shows how $\Gamma_{m}$ varies with spanwise distance. The mean secondary flow is clearly strongest near the shaft, which is due to the counter-rotating vortices that form there. Such vortices are also generated near the feather edges, but the riblets have a small height there, which yields only a relatively weak secondary flow. Compared to $\alpha=165^{\circ}$, textures with $\alpha=15^{\circ}$ have a stronger mean secondary flow at the shaft, which is probably due to the downdraft of high-momentum fluid there.

Figure 9 (bottom right) shows how $\Gamma_{t}$ varies with spanwise distance. The relatively uniform turbulence for the textures with $N_{\mathrm{groove}}=1,4$ shows that the flow is well-mixed. In contrast, the turbulence changes with $y$ for the other textures. When approaching the shaft, turbulence reduces when $\alpha=15^{\circ}$ and increases when $\alpha=165^{\circ}$. The same was observed for the streamwise velocity fluctuations in the previous subsection. When compared with the smooth wall, however, the turbulence at the shaft is increased for all textures. Away from the shaft, three textures exhibit local turbulence reductions.

Variants Figure 9 also shows the volume-averaged secondary flow strength for the shifted texture variation. The turbulent contribution is practically the same as for the not-shifted texture. The same applies to the mean contribution for $N_{\mathrm{groove}}=1$ and $4$. This agrees with the observation that, apart from the spanwise shift of the main vortex near the top wall, the mean secondary flow does not change. For $N_{\mathrm{groove}}=16$ and $32$, however, the mean secondary flow is stronger for the shifted texture, which is ascribed to the formation of a tertiary flow that extends from bottom to top wall.

The secondary flow for the texture variation with $\alpha=0^{\circ}$ is much weaker than for the herringbone riblet geometries, which is attributed to alignment of the riblets with the mean flow. Both $\left\langle\Gamma_{m}\right\rangle_{V}$ and $\left\langle\Gamma_{t}\right\rangle_{V}$ are still larger than the values for the parallel riblet geometry with constant blade height, although they seem to approach those values in the limit of large $N_{\mathrm{groove}}$. The figure shows that parallel riblet geometries can reduce turbulence.

5.5 Spanwise Transport of Streamwise Momentum

To confirm that advective transport is most important around the shaft, this subsection considers the streamwise momentum balance as function of spanwise distance. The Navier-Stokes equation for streamwise momentum (see Equation 1) is rewritten as $\partial u/\partial t=\sum f$, where $f$ is one of the terms in that equation. A time, streamwise and wall-normal average (denoted by $\left\langle\overline{f}\right\rangle_{xz}$) is applied to this equation. Assuming statistical stationarity and using the boundary conditions, the balance becomes:

[TABLE]

The terms in this equation only depend on the spanwise coordinate $y$. The third term is called ‘shear’, since it equals the sum of bottom and top-wall shear as a result of the wall-normal integration. All terms are divided by $d^{\mathrm{tot}}_{\mathrm{smooth}}$ and shown in Figure 10 for four herringbone textures. The abscissa represents the spanwise distance to the feather shaft. Most curves exhibit approximately equidistant wiggles that result from staircasing: the riblet height increases in 17 steps of size $\Delta z_{w}$ from 0 at the feather edges to $h$ at the shaft.

Before the differences between $N_{\mathrm{groove}}=4$ and $128$ are highlighted, the general behavior of the different terms is clarified. The driving term is almost constant with $y$, because $\gamma_{f}$ changes only marginally from the feather shaft to the edges. For $N_{\mathrm{groove}}=4$, this term is clearly larger than 1, indicating a significant drag increase.

The shear term is negative, as it tends to decelerate the fluid. Its magnitude quantifies how the flat-wall streamwise shear changes with the spanwise coordinate. When moving from the feather edges towards the shaft, the wall-shear term first decreases (in magnitude), which is attributed to shielding of the flat wall by blades which increase in height. This trend continues for textures with $\alpha=165^{\circ}$. However, the shear magnitude peaks near the shaft for $\alpha=15^{\circ}$, which is ascribed to the local downdraft that transports high momentum towards the wall. For very wide feathers ($N_{\mathrm{groove}}=128$), the shear term equals -1 at the feather edges. There the blade height is zero and the smooth-wall result is recovered.

The IBM term represents the streamwise drag force on the riblets. It is negative, as it is responsible for a velocity decrease (like the shear term). It equals zero at the feather edges, because the riblets have no height there. When approaching the feather shaft from the edges, the IBM-force magnitude first increases, which is due to riblet-height increase. Near the shaft, it has a local minimum for $\alpha=165^{\circ}$. In contrast, it is very large there for $\alpha=15^{\circ}$, likely due to the downdraft of high-speed fluid.

The balance for $N_{\mathrm{groove}}=4$ shows significant contributions from spanwise advective and diffusive transport, especially near the shaft. For $\alpha=15^{\circ}$, advection is on average responsible for an increase of streamwise momentum near the shaft. This is attributed to the secondary-flow vortex that transports low momentum away from the shaft (near the wall) and high momentum towards the shaft (closer to the channel centerline). In contrast, diffusion transports high momentum away from the shaft. These trends are opposite to that of textures with $\alpha=165^{\circ}$, for which spanwise advection causes a streamwise momentum decrease and diffusion an increase near the shaft.

Compared to $N_{\mathrm{groove}}=4$, the relative importance of the terms in Equation 13 is very different for $N_{\mathrm{groove}}=128$. The prominent peaks near the shaft are very narrow. Spanwise transport by advection and diffusion is close to zero for the largest part of the feather. Away from the shaft, the texture behaves as parallel riblets in yaw with a local balance between the driving force on the one hand, and the IBM and wall-shear force on the other hand.

5.6 Wall-normal Transport of Streamwise Momentum

The previous subsections show that drag increase is accompanied by a strong secondary flow, which suggests that enhanced advection is responsible for the drag augmentation. To underpin this suggestion, this subsection considers the streamwise momentum balance as function of the wall-normal coordinate. A time, streamwise, and spanwise average (denoted by $\left\langle\overline{f}\right\rangle_{xy}$) is applied to the Navier-Stokes equation $\partial u/\partial t=\sum f$. Assuming statistical stationarity and using the boundary conditions, the balance reads:

[TABLE]

The terms in this equation only depend on the wall-normal coordinate $z$. All terms are divided by $d^{\mathrm{tot}}_{\mathrm{smooth}}$ and displayed in Figure 11 for four cases. The profiles for parallel riblets without yaw (not shown) are very similar to that for riblets in yaw. In addition, the profiles for herringbone textures with $\alpha=15^{\circ}$ (not shown) are similar to the ones for $\alpha=165^{\circ}$. The behavior of the different terms is clarified below. The driving term is almost constant, as before.

The IBM term represents the drag force on the texture, so it is only present near the wall and it is responsible for a velocity decrease. The drag force is especially large near the blade tips. That explains the peak at blade height for riblets in yaw, for which all blades have the same height. The large IBM force near the blade tips is in Figure 11 not evident for herringbone textures. Instead, the peak is smeared out due to the spanwise blade-height variation between [math] and $h$ (see Equation 11). For the two herringbone textures shown in the figure, the narrower feather clearly experiences a larger IBM force than the wider feather.

Diffusive transport is especially important near the wall. It has a negative tendency for smooth walls, but for textured walls it becomes positive in between the riblets. The latter is associated with an inflection point in the mean streamwise-velocity profile.

Advective transport is significant throughout the whole channel. Near the channel centerline, it balances the driving force (for all cases). Near the wall, turbulent advective transport is responsible for a velocity increase and the associated drag augmentation. Parallel riblets suppress turbulent advection, as is clear from the second sub-figure. Compared to smooth walls, the peak of advection shifts upwards and shrinks. The herringbone case with $N_{\mathrm{groove}}=128$ shows the same trend, although the upward shift is less pronounced. Turbulent transport below $z=h$ is less suppressed as compared to parallel riblets in yaw, which is ascribed to the blade height decrease with spanwise distance to the shaft. Still, weakening and lifting of advective transport is quite apparent. In contrast, advective transport is much stronger for $N_{\mathrm{groove}}=4$ as compared to a smooth wall. This reinforces the suggestion that drag augmentation is caused by enhanced advection, which will be confirmed in the next section with a quantitative analysis.

6 Drag Change Decomposition

6.1 Previous Research

Fukagata et al. [28] derived an equation (called FIK-identity) that decomposes the frictional-drag coefficient into contributions from different dynamical effects. They distinguished four contributions, namely laminar (or bulk), turbulent, inhomogeneous and transient. For homogeneous and steady plane channel flow, the FIK-identity reads (in our notation):

[TABLE]

where it is used that the domain height $L_{z}=1$. Note that the skin-friction coefficient $C_{f}\equiv 2\overline{\tau_{w}^{*}}/\left(\rho^{*}U_{b}^{*2}\right)=\overline{f^{\mathrm{driving}}}$ for smooth-wall channel flow.

This relation has been extended to geometrically more complex surfaces by Peet and Sagaut [59]. They derived analytical relations for streamwise-, spanwise- and quasistreamwise-homogeneous surfaces. To our knowledge, however, these relations cannot be used for the herringbone textures. That asks for a second extension of the FIK-identity.

6.2 Derivation of Extended FIK-identity

The extended FIK-identity is derived from the slab-averaged streamwise momentum balance (Equation 14). For ease of notation, that equation is written as $0=\sum_{l}\left.f^{(l)}\right|_{z}$, where $f^{(l)}$ represents one of the terms in that equation. The subscript $z$ expresses that $f^{(l)}$ only depends on the wall-normal coordinate. In what follows, the global drag balance is needed. It follows from a single integration of Equation 14: $0=\sum_{l}\int_{0}^{1}\!\left.f^{(l)}\right|_{z}\,\mathrm{d}z$ or

[TABLE]

Next, Equation 14 is rewritten as $0=\sum_{l}(1/2)\left(\left.f^{(l)}\right|_{z}+\left.f^{(l)}\right|_{1-z}\right)$ to explicitly account for symmetry in the mean flow. Triple integration is applied to this equation, such that the extended FIK-identity in condensed form reads:

[TABLE]

Using the boundary conditions at the channel walls, the definition of the bulk velocity (i.e. $\int_{0}^{1}\!\left\langle\overline{u}\right\rangle_{xy}\!\,\mathrm{d}z=1$), integration by parts to transform multiple to single integrations, and the global drag balance (i.e. Equation 16), the last equation becomes:

[TABLE]

To arrive at an equation for $d^{\mathrm{tot}}=\langle\gamma_{f}\overline{f^{\mathrm{driving}}}\,\rangle_{V}=\overline{f^{\mathrm{driving}}}\int_{0}^{1}\!\left\langle\gamma_{f}\right\rangle_{xy}\,\mathrm{d}z$, Equation 18 is divided by the prefactor in square brackets and multiplied by $\int_{0}^{1}\!\left\langle\gamma_{f}\right\rangle_{xy}\,\mathrm{d}z$, which yields the final extended FIK-identity:

[TABLE]

Notice that the total advection term is split up into contributions from mean and turbulent flow, using that $\overline{uw}=\overline{u}\,\overline{w}+\overline{u^{\prime}w^{\prime}}$. The symmetry or antisymmetry with respect to the channel centerline at $z=1/2$ is evident from the factors in braces in Equation 18. As the multiplier of $\overline{f^{\mathrm{driving}}}$ in that equation depends on $\gamma_{f}$, the bulk term $d^{\mathrm{bulk}}$ is slightly texture dependent, as was also found by Peet and Sagaut [59]. For smooth walls, $\gamma_{f}=1$ and that multiplier equals 1/12, such that Equation 15 is recovered.

As Equation 19 applies to both smooth and textured walls, the drag change can be decomposed in a similar way:

[TABLE]

The term within braces is abbreviated as $dc^{(l)}$. It represents the change of a certain term for textured walls as compared to that term for smooth walls. This decomposition quantifies which terms contribute to drag reduction or increase.

The identity presented here is somewhat different from the one derived by Peet and Sagaut [59]. First, their decomposition only applies to skin friction, whereas the drag decomposition in Equation 19 also includes the pressure drag. Second, the IBM-term is not present in their identity. They used a body-fitted coordinate system in their derivation. As a result, the shear stress on the texture directly derives from integration of the viscous diffusion term. That approach has the added advantage that the skin-friction coefficient for simple textures in a laminar flow can be computed exactly based on purely geometrical considerations without performing the flow calculations [59]. However, their relation applies to quasi-homogeneous surfaces only and adopts a more complicated integration using a body-fitted grid. In contrast, the Cartesian integration that is employed here is not restricted to certain geometries.

The extended FIK-identity (Equation 19) is not only useful when an Immersed Boundary Method (IBM) is used. In the present work, $f^{\mathrm{IBM}}_{i}$ is a body force that models the shear and pressure forces that the texture exerts on the flow. However, the FIK-identity applies to any body force. Furthermore, the current drag decomposition is also applicable to body-fitted calculations. In that case, the obstacles should be considered as part of the domain, because the identity is based on integration over the entire rectangular channel volume. Three steps are required for a successful use of the decomposition in this situation.

1. The geometry should be translated into a three-dimensional phase-indicator function $\gamma_{f}$.
2. A zero-flow condition should be used for the obstacle volume.
3. The drag force on the obstacle surface should be translated into a three-dimensional body force or IBM force.

In view of the factors between braces in Equation 18, advection and the IBM force contribute differently to the total drag. The weighing factor for advection is largest near the wall, so significant advective transport near the wall contributes most to drag augmentation. In contrast, the weighing factor for the IBM force is largest near the channel centerline. Obstructing the flow there is for two reasons more detrimental than an obstruction near the wall: drag increases due to a larger flow velocity (so a larger IBM force) and a larger weighing factor.

6.3 Smooth Walls

Table 2 shows the decomposition of the total drag into the different contributions for smooth-wall turbulent channel flow. The results for three Reynolds numbers are shown. The total drag is taken as a reference. The rest term equalizes the left-hand and right-hand side of Equation 19. It results from the finite simulation time. It becomes slightly larger for higher $Re_{b}$, which is attributed to a relatively shorter simulation time. Its magnitude is typical for all the other simulations.

As the table shows, only bulk transport and turbulent advection contribute to the total drag in smooth-wall channel flows. The mean advection term is zero to machine precision, and $d^{\mathrm{IBM}}=0$. At $Re_{b}=5500$, the total drag comes for 27% from the bulk and for 73% from the turbulent advection term. These numbers are close to the 26% and 74% reported previously for $Re_{b}=5460$ [59]. As the Reynolds number increases, the relative contribution of turbulence to the total drag increases. This stresses once more the need for reduction of turbulent drag in high-Reynolds-number flows.

6.4 Parallel Riblets

Using Equation 20, the drag change for parallel riblets is decomposed into contributions from the bulk, mean advection, turbulent advection and IBM terms. Figure 12 exhibits this decomposition as function of $s^{+}$ and $\alpha$ in a bar graph. The total term equals the total drag change, which was shown already in Figure 2. The bulk term is in both graphs not visible, as it is very close to zero.

Change with blade spacing The IBM term, which represents the drag force on the blades, is responsible for an increase of the total drag. Apart from being positive, it also increases with $s^{+}$. When the riblet spacing increases, the riblet height increases as well ($h/s=0.5$ is fixed), but the total blade area per unit spanwise width remains constant. Therefore, the increase of the IBM term with $s^{+}$ cannot be due to an increased wetted area, but it is ascribed to blades that further protrude into the flow. Although it remains small as compared to the other terms, the contribution from mean advection also increases with blade spacing. That is attributed to a stronger mean secondary flow for larger $s^{+}$ (not shown here, see also e.g. [14, 34]). Finally, the turbulent drag contribution is negative, which indicates suppression of turbulent streamwise momentum transport. For the three $s^{+}$ values shown here, the maximum turbulent-drag suppression is almost 15% at $s^{+}=17$. The figure demonstrates that the optimum $s^{+}$ is a trade-off between an additional drag force on the blades and reduced turbulent transport because of the blades.

Change with yaw angle Figure 12 shows a similar decomposition for parallel riblets in yaw at fixed $s^{+}=17$. The IBM term increases for increasing yaw angle $\alpha$, which might be due to the additional pressure drag. The change in the mean advection term is very small. The contribution from turbulent advection is increasingly less negative when $\alpha$ increases. The figure thus indicates that deterioration of riblets in yaw is both due to an increased drag force on the blades and reduced suppression of turbulent transport.

6.5 Herringbone Riblets

Drag decomposition Figure 13 shows the decomposition of the total drag change into its contributions from the IBM force, turbulent advection and mean advection. The bulk and rest terms are not presented. They never exceed 0.04% (bulk) and 0.08% (rest) in magnitude. The top left graph is the same as Figure 5 and is shown here again for ease of comparison. Apart from the bulk contribution, the IBM term is the least important term for most textures. In the worst case, it causes an 8.1% increase of drag, which is comparable to results for standard parallel riblets (7.0% for $\alpha=0^{\circ}$, 8.6% for $\alpha=15^{\circ}$). Advection by the mean flow can be a significant contribution to drag increase, up to 28% for $N_{\mathrm{groove}}=4,\alpha=165^{\circ}$. The turbulent advection term exhibits the largest changes: it varies between a 6% reduction and a 50% increase.

The variation of the mean and turbulent drag contributions with feather width very much resembles that of the mean and turbulent secondary flow strengths (compare Figure 13 and Figure 9). It thus appears that advection and secondary flow go hand in hand. The same can be observed from the instantaneous flow field in Figure 14, which shows both contours of advection and vectors of secondary flow. Regions with strong secondary flow (in particular strong $w$) are also regions with strong advection. Both secondary flow and the advective flux are split into mean and turbulent contributions in view of $\Gamma_{s}^{2}=\Gamma_{m}^{2}+\Gamma_{t}^{2}$ and $\overline{uw}=\overline{u}\,\overline{w}+\overline{u^{\prime}w^{\prime}}$, respectively. In summary, the results clearly indicate an intimate connection between increased advective transport and stronger secondary flow.

Given the close correspondence between secondary flow and advection, the change of the secondary flow strength with feather width as clarified in subsection 5.4 also explains the trend of the advective drag contributions. In particular, the increased advective drag for smaller feather widths is due to the higher spanwise density of the converging or diverging riblets that generate the secondary flows. Also, reduction of the mean advective drag for $N_{\mathrm{groove}}=1$ is ascribed to a weaker mean secondary flow due to spanwise confinement. There is, however, one major difference between Figure 13 and Figure 9, namely the effect of shifting of the top-wall texture. For 16 and 32 grooves, the mean secondary flow for the shifted textures is clearly stronger than for the not-shifted textures. In contrast, the drag due to mean advection is about the same, independent of the shift. This is explained by the fact that a stronger secondary flow near the channel centerline does not contribute much to drag because of the factor $\left(1/2-z\right)$ in Equation 18. This demonstrates that mean-secondary-flow strength is a good indicator for the contribution of mean advection to drag, provided that the mean secondary flow near the centerline is weak.

Spanwise decomposition of drag change Figure 13 shows that drag reductions are possible in the limit of large $N_{\mathrm{groove}}$. Like for the conventional riblet texture, these reductions originate from weakened turbulent advective transport. To reveal the origin of this weakening, the spanwise dependence of the advective FIK-terms is investigated with a spanwise decomposition. Let $\left\langle\overline{f}\right\rangle_{xy}$ be one term in Equation 14 and $d^{(l)}$ the corresponding drag contribution. The calculation of this FIK-term can be rewritten as:

[TABLE]

The function $g=g(z)$ results from conversion of a triple to a single integral, and normalization. The function $d_{sp}^{(l)}=d_{sp}^{(l)}(y)$ represents the spanwise decomposition of the FIK-term, as indicated by the subscript $sp$. The drag change contribution $dc^{(l)}$ can be decomposed in a similar way:

[TABLE]

where $dc_{sp}^{(l)}$ quantifies how the drag change depends on the spanwise coordinate. Note that $d_{sp}^{(l)}$ and $dc_{sp}^{(l)}$ have a spanwise dependence, whereas $d_{\mathrm{smooth}}^{(l)}$ and $d^{\mathrm{tot}}_{\mathrm{smooth}}$ have not. The spanwise decomposition of the mean and turbulent advective terms is shown in Figure 15 for $\alpha=15^{\circ}$ and $165^{\circ}$, and four feather widths. Small asymmetries with respect to $y=y_{s}$ are attributed to slow convergence of weak secondary flows. The small oscillations in the curves for $N_{\mathrm{groove}}=128$ probably appear for the same reason.

The effect of convergent/divergent riblets can be isolated from that of parallel riblets in yaw for the textures with $N_{\mathrm{groove}}=128$. The parallel riblets dominate the flow in a region sufficiently far away from the shaft. In that region, the mean advection term fluctuates around zero. The turbulent term is zero at the feather edges. Riblets have no height there and the smooth-wall result is recovered. Away from the edges, the term decreases more or less linearly, which is attributed to riblet-height increase and the associated suppression of turbulent transport. It approaches approximately $-10$%, which belongs to full-height parallel riblets with yaw angle $\alpha=15^{\circ}$ or $165^{\circ}$. These favorable trends of mean and turbulent advective drag stop near the shaft because of the strong secondary flow there.

The converging/diverging riblets dominate the flow in a region around the shaft. That is especially evident from the drag change due to mean advection, which shows a clear signature of the mean flow described before, in particular the strong secondary-flow vortices near the shaft. For $\alpha=165^{\circ}$ and any $N_{\mathrm{groove}}$, the updraft of fluid around the shaft results in a local drag reduction, as is apparent from the dip in $dc_{sp}$ at $y=y_{s}$. The peak next to this dip is associated with that part of the vortex that transports momentum towards the wall. For $N_{\mathrm{groove}}=16$, a second dip appears, which is ascribed to the tertiary flow shown in Figure 7. The central dip is narrower for larger $N_{\mathrm{groove}}$, because the secondary-flow vortex near the shaft is smaller compared to the total feather width. A similar (but opposite) description applies to textures with $\alpha=15^{\circ}$.

Local turbulent drag reduction can be obtained by diverging riblets. For $\alpha=15^{\circ}$ and $N_{\mathrm{groove}}\geq 16$, the strong wall-directed mean flow is accompanied by reduced turbulent transport. The reverse is true for $\alpha=165^{\circ}$, namely that the strong wall-leaving mean flow is accompanied by increased wall-directed turbulent transport, as can also be seen from Figure 14. These observations fully agree with the findings for uniform blowing or suction. In Direct Numerical Simulations, it has been found that uniform blowing reduces mean advective drag and enhances turbulent drag, while uniform suction enhances mean advective drag and reduces turbulent drag [43].

The local contribution of advection to drag might be very different from its global (or volume-averaged) contribution. For instance, mean advection might seem much more important than turbulent advection in view of the scales of Figure 14 and Figure 15. However, the volume-averaged turbulent term is often at least as important as the mean term. As a second example, the texture with $\alpha=165^{\circ}$ and $N_{\mathrm{groove}}=4$ exhibits locally a drag reduction as high as 400% due to a strong wall-leaving flow. However, one should realize that the favorable updraft of low-momentum fluid is compensated by an adverse downdraft of high-momentum fluid. The influence of the whole vortex on the drag should be considered. In the case of $\alpha=165^{\circ}$ and $N_{\mathrm{groove}}=4$, the vortex covers one feather half. The corresponding spanwise-averaged mean advective drag is 28%, which establishes again an unfavorable effect of the vortices near the shaft on the drag. Therefore, one should be careful to judge the performance of this (or any) texture based on a local drag determination. The mean or turbulent advective drag might be reduced locally. However, when the flow is dominated by strong advection (such as near the shaft), the volume-averaged drag generally increases due to an overall increase of both the mean and turbulent advective contributions.

In summary, the present study confirms two effects of the herringbone riblet texture on the turbulent drag. The first effect relates to the texture-generated secondary flows around the shaft. Although turbulent drag might be reduced locally due to a wall-directed mean flow, the overall trend is an increased turbulent drag because of the fluctuating secondary flows that are generated by the converging/diverging riblets near the shaft. The second effect relates to the parallel riblets sufficiently far away from the shaft. Those riblets suppress turbulent transport, which results in turbulent drag reduction.

7 Discussion

Detrimental effect of convergent/divergent riblets The results presented in this paper clarify the influence of the herringbone riblets on the drag. The drag increase for $\Lambda_{f}/L_{z}\lesssim O(1)$ is attributed to an increase of advective transport (Figure 13), which in turn is associated with the secondary-flow vortex near the shaft (Figure 15). Although the vortex might locally be responsible for a drag reduction, as a whole it is responsible for a drag increase. That vortex originates from the herringbone riblets near the shaft (Figures 7, 9, 10). That strongly suggests that these converging/diverging riblets are detrimental to DR. That is confirmed by the finding that $\Lambda_{f}\rightarrow\infty$ is most beneficial for DR (Figure 5). In that limit, the herringbone texture approaches the conventional parallel-riblet texture in yaw. Hence, the presence of convergent/divergent riblets in the texture seems unfavorable for DR.

Research on spanwise-periodic spanwise forcing Research on the drag-reducing spanwise forcing further underpins the detrimental effect of the herringbone riblets on the drag. The spanwise traveling wave of spanwise body force is described by $f_{y}=A(z)\sin\{(2\pi/\Lambda_{f})y-\omega t\}$. Here $f_{y}$ represents the spanwise body force, $A(z)$ the forcing amplitude (only nonzero in the vicinity of the wall), $\Lambda_{f}$ the spanwise wavelength of the forcing, and $\omega$ the angular frequency of the forcing. The herringbone texture is similar to this forcing with respect to its spanwise periodicity, though the texture is static ($\omega=0$). Research shows that the best-performing spanwise traveling wave (in terms of DR or net energy saving) is the one with infinite wavelength, i.e. the spanwise wall oscillation [21, 63]. In other words, $\Lambda_{f}\rightarrow\infty$ is most beneficial for DR, like for the herringbone texture. Hence, the addition of spanwise periodicity seems detrimental to the drag-reducing performance of both the spanwise forcing and the riblet texture in yaw.

Alternative passive riblet textures As the herringbone texture with a spanwise variation of the forcing appears to be unprofitable for DR, a streamwise variation of the forcing might be more advantageous. Indeed, streamwise traveling waves of spanwise forcing are superior to spanwise traveling waves of spanwise forcing, presumably because of the unaltered wall-normal gradient of the wall-normal velocity component at the wall [63]. Both experimentally and numerically, DR has been obtained with a stationary streamwise variation of spanwise forcing [82, 36]. DR by streamwise variation of wall-normal forcing (e.g. suction and blowing) has been reported as well, although the net energy saving has been small [64, 51].

The drag increasing/reducing trends in the current study are similar to what has been found for uniform suction/blowing (Figures 8, 15), but a net drag-reducing effect could not be confirmed. It is difficult (if not impossible) to eliminate the unfavorable effects of downdrafts, as mass conservation dictates that updrafts need to be compensated by downdrafts. Also, textures with converging/diverging riblets are not fully comparable to uniform blowing/suction. For instance, uniform blowing originates from a nonzero mass flux through the wall, in contrast to the texture-generated updrafts. Still, textures with streamwise variation of the forcing seem more promising than the herringbone texture.

Comparison with experiments The conclusion that the herringbone texture seems detrimental to turbulent drag reduction apparently contrasts with the experimental study of Chen et al. [11]. Although it has not been the aim of the present study to reproduce their experiments numerically, a comparison might still be illuminating. Chen et al. [11] obtained drag reduction in forward flow for $\Lambda_{f}/D\lesssim 1$ ($D$ being the pipe diameter), with a maximum of 20%. In contrast, the present study only achieved DR for $\Lambda_{f}/L_{z}>10$, with a maximum of 2% in backward flow. These contrasting results might be ascribed to differences in riblet texture and Reynolds number.

First, the textures were different, particularly the riblet shape and angle. The feather width might have been different too, but it is unclear what value for $N_{\mathrm{groove}}$ was used in the experiments. The experimental texture consisted of sawtooth riblets at an angle of $30^{\circ}$ with the flow direction. The numerical texture was composed of blades at an angle of $15^{\circ}$. Although blades seem to be more sensitive to yaw (see Figure 2), they were studied at a smaller yaw angle. Therefore, there is currently no clear indication that the use of blade riblets in the numerical study contributed to a lesser drag-reducing performance of the herringbone texture.

Second, the Reynolds numbers differed significantly. The 20% DR was obtained at $Re_{b}\approx 2.6\cdot 10^{5}$ (based on pipe diameter and bulk velocity), while the numerical study was performed at $Re_{b}=5500$. For DR techniques that rely on near-wall flow manipulation, the attainable DR is approximately independent of the Reynolds number (although not fully [40, 41, 75, 33]) when near-wall scaling is applied. However, herringbone riblets (in contrast to conventional riblets) cannot be regarded as viscous-region modifiers, because they generate a secondary flow throughout a large part of the channel (Figure 7). Therefore, viscous scaling is probably inappropriate in this case. Because of the unknown scaling and the Reynolds number difference, the texture parameters (e.g. $\Lambda_{f}$) and the drag reduction results of the experiments and the simulations cannot be directly compared.

8 Conclusions and Outlook

The drag-reducing performance of a herringbone riblet texture was studied with Direct Numerical Simulations (DNSs) of turbulent flow in a channel with height $L_{z}$. The FIK-identity for drag decomposition was extended to textured walls and was used to study the underlying drag change mechanisms. For validation, simulations with smooth walls, parallel blade riblets and parallel blade riblets in yaw were performed, which showed good agreement with literature. The parallel-riblet simulations exhibited an expected but small low-Reynolds-number effect. A maximum drag reduction (DR) of 9.3% was obtained, close to the 9.9% that has been found experimentally [5].

The herringbone texture can both increase or reduce the drag, depending on the spanwise texture wavelength $\Lambda_{f}$. For $\Lambda_{f}/L_{z}\lesssim O(1)$ (i.e. narrow feathers), the drag increases with a maximum of 73% for $\Lambda_{f}/L_{z}=0.4$. This increase is ascribed to the convergent/divergent riblets. They generate a fluctuating secondary flow, which on average consists of two counter-rotating vortices centered above the regions of riblet convergence/divergence. The strong secondary flow increases both mean and turbulent advective transport, which in turn results in the significant drag increase.

A slight drag reduction of 2% was found for $\Lambda_{f}/L_{z}\gtrsim O(10)$ (wide feathers). Due to the large feather width, the secondary flow generated by the converging/diverging riblets now influences only a relatively small part of the whole texture. Its drag-increasing contribution is therefore small. The largest part of the texture behaves similarly to a conventional parallel-riblet texture in yaw. Specifically, suppression of turbulent advective transport is responsible for the small DR that was obtained.

As was found by other researchers for spanwise traveling waves of spanwise forcing, the current study confirms that $\Lambda_{f}\rightarrow\infty$ is most beneficial for DR. In that limit, the texture approaches the conventional parallel-riblet texture in yaw. Therefore, the presence of convergent/divergent riblets in the texture seems detrimental to turbulent drag reduction, which apparently contrasts with the experiments of Chen et al. [11]. However, differences in Reynolds number and texture parameters (riblet shape, feather width, angle between riblets and flow direction) hindered a one-to-one comparison between the present simulations and the experiments.

More elaborate experiments and simulations are required to further investigate the drag-reducing potential of the herringbone texture. First, the maximum DR of 20% should be reproduced and the optimum texture parameters ($\Lambda_{f}$, $s$, $h/s$, $\alpha$) should be determined. In numerical simulations, the more realistic sawtooth riblet geometry should be implemented. The Reynolds-number influence and the parameter scaling also need further attention. As the herringbone riblets generate a secondary flow throughout a large part of the channel, the scaling of drag reduction with the texture and flow parameters is presumably non-trivial.

Whether feather riblets have an aerodynamic function remains an open question. Chen et al. [11] claimed that feather riblets greatly impact flight performance because of drag reduction. The DR was ascribed to suppression of turbulent momentum transport, but evidence has been inconclusive so far. One should realize that wings of birds are not flat and operate at a relatively low Reynolds number, so flow-separation delay seems a more plausible aerodynamic function of feather riblets. Indeed, several studies confirm that roughness on a bird wing contributes to separation control [8, 48, 79]. Furthermore, separation delay has been obtained with vortex generators that resemble the herringbone texture [49]. So, future studies might investigate the potential of the herringbone texture for flow-separation control.

Acknowledgements

The research leading to these results has received funding from the European Union Seventh Framework Programme in the SEAFRONT project under grant agreement nr. 614034. The simulations were performed on the Dutch national supercomputer Cartesius at SURFsara, Amsterdam, The Netherlands. This work was sponsored by NWO Physical Sciences for the use of supercomputer facilities. H.O.G.B. would like to thank his colleague Pedro Costa for his supercomputational support.

Appendix A Time Advancement at Fixed Bulk Velocity

For simplicity, Equation 1 is rewritten as:

[TABLE]

where $r_{i}$ contains the advection, diffusion and IBM terms. Let the integer $n$ denote the time steps. The used Runge-Kutta scheme (RK3) employs three sub-steps, which are numbered by the integer $s$. RK3 introduces intermediate velocities $u_{i}^{(s)}$, where $u_{i}^{(0)}=u_{i}^{n}$ and $u_{i}^{(3)}=u_{i}^{n+1}$. Similarly, intermediate pressures $p^{(s)}$ and driving forces $f^{\mathrm{driving}(s)}$ are introduced. Because of the Crank-Nicolson scheme for pressure, $p^{(s)}=p^{(s-1)}+\widetilde{p}^{(s)}$ with correction pressure $\widetilde{p}^{(s)}$. The time advancement is illustrated here for an arbitrary sub-step $s$:

[TABLE]

where $u_{i}^{*({s})}$ is the first and $u_{i}^{**({s})}$ the second prediction velocity. Note that the asterisk here is not used to denote dimensional quantities. The parameters $\alpha_{(s)}$, $\gamma_{(s)}$ and $\zeta_{(s-1)}$ are RK3 parameters (see e.g. [89]). To obtain the driving force, the equation for the streamwise velocity is volume-averaged, which yields:

[TABLE]

The first term represents the intermediate bulk velocity in sub-step $s$, which is set equal to one to obtain the constant bulk flow. The second term disappears because of periodic boundary conditions. The first prediction velocity is known, so the third term can be computed. The resulting equation can be solved for $f^{\mathrm{driving}(s)}$, which yields:

[TABLE]

Next, $u_{i}^{**({s})}$ is computed. With use of the continuity equation, the divergence of Equation 24 yields a Poisson equation for the correction pressure:

[TABLE]

When this is solved, the updated velocity and pressure are computed:

[TABLE]

This procedure guarantees that the bulk velocity equals one in each RK3 sub-step, which results in three values of the driving force. The total forcing (i.e. one per time step) is calculated as follows:

[TABLE]

Appendix B Immersed Boundary Method for Blade Riblet Textures

The IBM adjusts diffusive and advective fluxes around blades with help of indicator functions, as is described below for parallel riblets without and in yaw. The third subsection explains how the $f^{\mathrm{IBM}}_{i}$ term in Equation 1 is computed. The final subsection describes how well the IBM approximates the no-slip and no-penetration conditions at the texture surface.

B.1 Parallel Riblets without Yaw

Figure 16 shows the numerical grid used for parallel riblets without yaw. As a staggered grid is used, the velocity vectors are located at the faces of the grid cells. A few of these vectors are shown in the right figure. The blade coincides with a few grid points of the spanwise velocity $v$. The right figure helps to explain the IBM. The red boxes in that figure indicate three staggered grid cells that belong to the velocity components $u_{\mathrm{left}}$, $v_{(j,k)}$ and $w_{\mathrm{left}}$. The IBM implements an adjustment of the advective and diffusive fluxes at the cell faces marked with a thick solid line, as is explained below.

The streamwise velocity $u$ was only adjusted in grid cells next to the blades (both on the left and right side). For the grid cell of $u_{\mathrm{left}}$, the two spanwise fluxes through its thick solid face were changed: the advective flux $uv=0$ and the diffusive flux $\partial u/\partial y=-2u_{\mathrm{left}}/\Delta y$.

The spanwise velocity $v$ was changed for two reasons, namely to enforce no-penetration at the blades and to adjust the fluxes above the blade tip. To enforce no-penetration, the prediction velocity was set to zero at the grid points that coincide with a blade. Let $v_{\mathrm{blade}}$ represent a spanwise velocity component that coincides with a blade (see Figure 16) and let ∗ represent the first prediction velocity (as in Equation 24), then $v_{\mathrm{blade}}^{*}=0$. The actual velocity $v_{\mathrm{blade}}$ follows from the correction step (Equation 28) and is very close to zero, although not exactly zero (details follow in subsection B.4).

The second change to $v$ comprises the adjustment of advective and diffusive fluxes in the grid cell just above the blade tip (shown in red in Figure 16). The vertical fluxes of $v$ at the bottom face of that cell were adjusted. The diffusive flux $\partial v/\partial z$ was split into two contributions, namely from the left and right side of the blade (indicated by $\phi_{\mathrm{left}}$ and $\phi_{\mathrm{right}}$ in the figure). It accounts for the fact that the thin blade does not inhibit vertical transport. Specifically, $\partial v/\partial z=0.5\phi_{\mathrm{left}}+0.5\phi_{\mathrm{right}}$. The flux $\phi_{\mathrm{right}}$ was computed by linear interpolation of the four velocity components labeled with indices in the figure: $\phi_{\mathrm{right}}=(v_{t}-v_{b})/\Delta z$ with $v_{t}=0.75v_{(j,k)}+0.25v_{(j+1,k)}$ and $v_{b}=0.75v_{(j,k-1)}+0.25v_{(j+1,k-1)}$. The flux $\phi_{\mathrm{left}}$ was computed in a similar way. As advection near the blade tips is likely less important than diffusion, the advective flux was not split but simply set to zero (i.e. $vw=0$ at the thick solid face).

The wall-normal velocity $w$ was only adjusted in grid cells next to the blades, similarly to what was done for $u$. For the grid cell of $w_{\mathrm{left}}$, two spanwise fluxes at the thick solid face were changed: the advective flux $vw=0$ and the diffusive flux $\partial w/\partial y=-2w_{\mathrm{left}}/\Delta y$. These adjustments were also applied for the $w$-cell next to the blade tip, so this diffusive flux was not split into two contributions (although the blade covers only half of the cell face). To justify this choice, the simulation with $s^{+}=24$ at $Re_{b}=5500$ was repeated. The diffusive flux near the blade tip was separated into two contributions, namely from above and below the blade tip. No significant difference in drag was found.

B.2 Parallel Riblets in Yaw

To simulate the turbulent flow over blades in yaw, the riblet texture was rotated with respect to the grid, see Figure 17 (left). That required a different Immersed Boundary Method (IBM), as the surfaces were not anymore aligned with the Cartesian directions. As the left figure shows, the grid was generated in such a way that the obstacle surface still coincides with part of the staggered velocity vectors. The main disadvantage of this approach is that the streamwise and spanwise grid spacings cannot be chosen independently: they depend on $\alpha$. However, the major advantage is a relatively simple IBM. Due to the specific alignment of grid and texture, the number of grid cells per groove $N_{cg}$ is equal for the streamwise and spanwise directions. The grid shown in Figure 17 has $N_{cg}=16$.

The IBM is explained with use of Figure 17 (right). Staggered grid cells for $u_{\mathrm{top}}$ and $w_{\mathrm{left}}$ are shown as red boxes. The IBM implements an adjustment of the advective and diffusive fluxes at the cell faces marked with a thick solid line, as is explained below.

The streamwise and spanwise velocities were adjusted in a similar way, so only the change to $u$ is described here. The changes to $u$ were almost identical to that for $v$ in the previous subsection. To enforce no-penetration, the streamwise prediction velocity was set to zero (i.e. $u_{\mathrm{blade}}^{*}=0$) at the grid points that coincide with a blade. At the bottom cell face of the grid cell for $u_{\mathrm{top}}$, both the advective and diffusive fluxes were modified. The advective flux was set to zero: $uw=0$. The diffusive flux was not split into two contributions: $\partial u/\partial z=2u_{\mathrm{top}}/\Delta z$.

For the wall-normal velocity, two indicator functions were used, namely one for streamwise and another one for spanwise transport of wall-normal momentum. As the treatment of streamwise and spanwise transport of $w$ is analogous, only the former is described here. The wall-normal prediction velocity was not changed directly, as $w$ never coincides with blades. Therefore, only the fluxes in $w$-cells next to the blades were modified. Specifically, the fluxes at the right cell face of the grid cell for $w_{\mathrm{left}}$ were adjusted as follows: $uw=0$ and $\partial w/\partial x=-2w_{\mathrm{left}}/\Delta x$.

B.3 Computation of $f^{\mathrm{IBM}}_{i}$

Although the IBM forcing $f^{\mathrm{IBM}}_{i}$ appears as a separate term in Equation 1 and following equations, it is not implemented as an explicit relation. The two previous subsections have shown instead that the numerical code for smooth walls was adjusted to obtain the zero velocity at the blades. However, the IBM term was needed as a separate term for the postprocessing. Its computation is illustrated with reference to Appendix A. The first prediction velocity $u_{i}^{*({s})}$ that appears there contains all adjustments that are needed to obtain the zero velocity at the blades. Next, all these adjustments are left out to compute the prediction velocity $u_{i}^{*({s})\text{ no IBM}}$ in absence of the IBM forcing. The intermediate IBM forcing $f_{i}^{\mathrm{IBM}(s)}$ then follows from $u_{i}^{*({s})}=u_{i}^{*({s})\text{ no IBM}}+\alpha_{(s)}\Delta tf_{i}^{\mathrm{IBM}(s)}$. Finally, the total IBM forcing $f^{\mathrm{IBM}}_{i}$ is obtained from summing $\alpha_{(s)}f_{i}^{\mathrm{IBM}(s)}$, similar to Equation 29.

B.4 Boundary Conditions at the Texture Surface

The no-slip and no-penetration conditions at the texture surface need to be sufficiently satisfied for the IBM to be accurate. Let a penetration velocity denote the absolute value of a texture-collocated velocity component. As $u_{\mathrm{blade}}$ and $v_{\mathrm{blade}}$ denote two velocity components that coincide with the texture, then the corresponding penetration velocities are $u_{\mathrm{pen}}=\left\lvert u_{\mathrm{blade}}\right\rvert$ and $v_{\mathrm{pen}}=\left\lvert v_{\mathrm{blade}}\right\rvert$. Note that $w_{\mathrm{pen}}$ does not exist, as $w$ does not coincide with the textures studied in the present work. To check whether the no-penetration condition was sufficiently satisfied, each tenth simulation time-step three penetration-related velocities were computed, namely $\left\langle u_{\mathrm{pen}}\right\rangle$, $\left\langle v_{\mathrm{pen}}\right\rangle$ and $\max{(u_{\mathrm{pen}},v_{\mathrm{pen}})}$, where $\left\langle...\right\rangle$ denotes an average. The resulting time series were time averaged, yielding $\overline{\left\langle u_{\mathrm{pen}}\right\rangle}$, $\overline{\left\langle v_{\mathrm{pen}}\right\rangle}$ and $\overline{\max{(u_{\mathrm{pen}},v_{\mathrm{pen}})}}$. The magnitude of these penetration velocities (normalized with the bulk velocity) is similar for all the simulations, namely a mean on the order of $10^{-6}$ to $10^{-5}$ and a maximum on the order of $10^{-4}$ to $10^{-3}$. This shows that the penetration velocities are very small, so the no-slip and no-penetration conditions are sufficiently satisfied.

Figure 18 presents a streamwise velocity profile as function of the spanwise distance for a herringbone riblet texture with 4 grooves per feather half (so $N_{\mathrm{groove}}=4$). The left subfigure shows the texture together with a line parallel to the $y$-axis. The right subfigure shows the time-averaged streamwise velocity profile that was extracted along that line. Clearly, the velocity is zero at the 8 riblet locations. This illustrates that the boundary conditions at the riblet surfaces are satisfied.

Appendix C Simulation Parameters

An overview of all simulations with the corresponding parameters can be found in Table 3.