Drag reduction by herringbone riblet texture in direct numerical simulations of turbulent channel flow

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 Λ_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 Fukagata-Iwamoto-Kasagi identity for drag decomposition was extended to textured walls and was used to study the drag change mechanisms. For Λ_f/L_z ≳ 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 Λ_f/L_z ≲ 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 [1], gas injection [2], wall mass transfer [3], wall cooling or heating [4–6], wall motion [7], wall deformation [8] and electromagnetic forcing [9]. Passive methods include the change of surface chemistry (e.g. superhydrophobicity [10]), elasticity (e.g. compliant walls [11]), shape (e.g. airfoil or ship hull shape) and texture (e.g. riblets [12]). Surface roughness and surface textures have been successfully used for DR by transition delay in laminar flow [13], separation delay in turbulent flow over curved surfaces [14,15] and turbulence modification in turbulent flow over flat surfaces.

Of the investigated flat plate textures – such as sifted sand grains [16], transverse square grooves [17], dimples [18] and V-shaped protrusions [19] – 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 [21]. In the search for even higher values of DR, many variations on the standard riblet geometries have been investigated [12], such as hierarchical or compound riblets [22], riblets on a spanwise travelling surface wave [23], oscillating riblets [24–26], riblets in a wave-like pattern (either in phase [27] or out of phase [28]) and riblets combined with drag-reducing polymers [29]. 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 colouration (e.g. for camouflage or visual signals). The following description of the feather morphology is derived from Chen et al. [30,31,32]. 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 α ≈ 30°. The groove spacing s remains approximately constant with s⁺ ≈ 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. [31] is one of the first to investigate the herringbone riblet texture for turbulent DR. 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 moulding, 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. [19] 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 (∼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 [33] or transition delay [34]. Furthermore, several studies confirm that roughness on a bird wing contributes to separation control [35–37].

Herringbone riblet textures have been studied for several other reasons, such as for heat-transfer enhancement [38,39] and mixing of laminar flows in microchannels [40]. The hot-wire study by Koeltzsch et al. [41] 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. [42] 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 visualised in a flat-plate laminar boundary layer over convergent riblets [43]. So, the motivation for herringbone riblet studies has been DR, heat-transfer enhancement, mixing improvement, flow control and large-scale flow structures manipulation.

The study by Nugroho et al. [42] 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 [44–47]. 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 [48], a regular array of cubes [49], a bed with two parallel lanes of different roughnesses [50], the irregular surface of a replica of a damaged turbine blade [51], streamwise strips of elevated roughness [52] and the aforementioned herringbone riblet texture [41,42]. Interestingly, Schoppa and Hussain [53] obtained 20% DR in direct numerical simulations (DNSs) 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 DNSs. The numerical methods are described in Section 2. Section 3 explains how DR 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, 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 nondimensionalised 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*_bL*_z/ν*, with kinematic viscosity ν^*. The superscript ⁺ is used for nondimensionalisation with ν and , with wall shear stress τ_w and fluid density ρ. Nondimensionalisation 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 × 2.5 × 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 [54–57]. In addition, it is full-span, as opposed to the recently reconsidered minimal-span channels [58]. Finally, it is larger than the moderate box of size (π × π/2 × 1) that is large enough to reproduce the one-point statistics of larger boxes [59].

Texture: Textures are applied to the inside of both channel walls to enforce symmetry in the mean flow [60]. Unless stated otherwise, the top wall texture is the bottom wall texture mirrored in the centreline plane. To simulate flow over a non-smooth surface, two methods can be adopted: coordinate transformation or the immersed boundary method (IBM) [56]. Both methods have been applied to simulate turbulent flow over riblet walls [61,62]. We used an IBM similar to the one employed by Breugem and Boersma [55] and Pourquie et al. [63], which is based on Fadlun et al. [64]. The IBM forcing is direct, i.e. a forcing term is added to the discretised 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 Δx and Δ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 Δ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 centreline. The maximum vertical grid spacing occurs at the centreline and is called Δz_c.

The grid-cell size is important to correctly resolve small-scale fluid motions. The recommendations of Vreman and Kuerten [65] for finite difference codes were followed, namely Δx⁺ = 4.4, Δy⁺ = 2.9, Δz_w⁺ = 0.49 and Δz_c⁺ = 2.2 as maximum grid spacings for smooth-wall turbulent flows. For textured walls, the spanwise grid spacing was reduced to Δy⁺ ≈ 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: (1) (2) 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 ^driving the spatially uniform forcing term to obtain a constant bulk velocity. The phase-indicator function γ_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 δ_i1 guarantees that fluid is driven in the streamwise direction with i = 1. Appendix A describes how f ^driving is calculated. The IBM forcing f^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 discretised using the finite-volume method combined with a pressure-correction scheme [66]. 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 [67]. This discretisation procedure yields a Poisson equation, which is solved using a non-iterative solver based on fast Fourier transforms (FFTs). 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 × 1 × 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 φ = φ(x, y, z, t) represent an arbitrary flow variable. The following averages were used: (3) (4) (5) Here, is a time average over the statistically stationary part of the signal φ (which starts at t_stat and has duration T). An average over one spatial coordinate is denoted by . For instance, ⟨φ⟩_x is a streamwise average. A similar notation is adopted for an average over two spatial coordinates. For example, ⟨φ⟩_xy is a streamwise and spanwise average. Finally, ⟨φ⟩_V represents an average over the entire volume V = L_xL_yL_z. For instance, the bulk velocity is defined by U_b = ⟨u⟩_V.

3. Drag reduction quantification

This section describes the quantification of DR. 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 DR.

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

The quantification of DR 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*_bL_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: DR cannot result from an increased fluid volume or a locally increased channel height [69,70].

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 = ρ*(U_b*)²(L*_z)^{− 1}V^{− 1}∫_Vγ_ff^drivingdV. Since the dimensional prefactor is equal for the smooth- and textured-wall channel flows, the drag change is given by (6) 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 ≡ d^tot_textured and d_s ≡ d^tot_smooth. Let u_φ for now denote the uncertainty in φ. Given the independence of d_t and d_s, the uncertainty in DC follows from the law of error propagation: (7) The uncertainties in d_t and d_s were computed using the method outlined by Hoyas and Jiménez [71]. 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. discretisation 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_τ = u_τδ/ν, where δ = L_z/2 is the half-channel height. Compared to recent DNSs reaching Re_τ = 4000 [72], Re_τ ≈ 4200 [59] and Re_τ ≈ 5200 [73], the simulations in this study are considered to be standard. Therefore, this subsection suffices to validate drag and flow statistics.

Table 1. Validation of drag and flow statistics of smooth-wall simulations. The relative deviation of Re_τ from that predicted by Pope's relation is given by (Equation (8(8) )). Statistics of U, u_rms, v_rms, w_rms and p_rms are compared with simulations by Vreman and Kuerten [65,57] at and 590. The relative difference between current and reference flow statistic Q is measured with the root-mean-square relative deviation δQ (Equation (9(9) )).

Reb	Reτ			δU (%)	δurms (%)	δvrms (%)	δwrms (%)	δprms (%)
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

Drag is validated by comparing Re_τ with the value predicted by an approximate relation [74]. The deviation of Re_τ from is quantified using (8) Table 1 lists for all smooth-wall simulations. Although the relation for is approximate, good agreement is obtained for all Re_b.

Flow statistics of U, u_rms, v_rms, w_rms and p_rms were compared with simulations at Re_τ = 180 [65] and Re_τ = 590 [57]. Here, , and similarly for the other root-mean-square quantities. Let Q(z) be one of these statistics, then the root-mean-square relative deviation (9) is used to quantify the difference between Q (current) and Q^ref (reference). It was computed after piecewise cubic spline interpolation of Q and Q^ref to a uniform grid, z⁺(k) = k for integers 1 ≤ k ≤ k_max, with k_max = 175 at the lowest and k_max = 587 at the highest Re_b. Table 1 shows δ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_rms at the lowest Re_b is attributed to the difference between Re_τ and .

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% [21]. An overview of previous research can be found elsewhere [12]. A thorough DNS study has been performed quite recently [60,76,77].

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 C1). The grid resolution (specified in wall units) is about the same for all cases.

Figure 1. Parallel blade riblet texture. Left: Bird's-eye view of the texture aligned with the mean flow, showing six unit cells in the spanwise direction. Right: Bird's-eye view of the texture in yaw with yaw angle α, showing six unit cells in the spanwise and three in the streamwise direction.

Figure 2 (left) compares the simulated drag change as function of s⁺ with experiments performed by Bechert et al. [21]. The top axis shows the experimental bulk Reynolds number Re_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.

Figure 2. Validation of drag change for the parallel blade riblet texture. Left: Drag change as function of riblet spacing in wall units (no yaw). Numerical results at three different Reynolds numbers are compared with experimental data of Bechert et al. [21]. The experimental bulk Reynolds number is denoted by Re_Bechert. Right: Drag reduction loss (Equation (10(10) )) as function of yaw angle for different riblet geometries at s⁺ = 17. EXP 2001 are experiments by Hage et al. [75] at Re_b ≈ 14900. SIM 2012 are DNSs by Grüneberger et al. [27] at Re_b = 5750. SIM are the current simulations. The numbers in parentheses in the legend represent the maximum drag reduction DR_max in percent. The numbers in the figure represent the drag change DC in percent for the simulated blade riblet geometry.

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 DR data below Re_Bechert ≈ 10000 deviated more and more from previous high-Reynolds-number data [21]. 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 δ/h from 35 to 15. The blades protrude farther into the channel, which is presumably detrimental to DR. This is supported by the finding that, for δ/h ≲ 50, the effect of roughness extends across the boundary layer and the original wall flow dynamics is changed significantly [78].

The simulations capture the experimental trend quite well when Re_b is close to Re_Bechert. Especially good agreement is obtained around and below the optimum spacing. At s⁺ = 17, a maximum DR 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. [41]. The study of Hage et al. [75] 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. [27] have performed DNSs at Re_b = 5750 (Re_τ = 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 α was varied for fixed Re_b = 5500. The horizontal grid resolution was comparable to that used for aligned parallel riblets: Δx⁺ ≤ 4.1 and Δy⁺ ≤ 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 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 α (indicated by the numbers in the figure). Clearly, the DR 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 is introduced [75]: (10) where DR = −DC, and is the maximum DR that can be obtained with a given geometry: . It was assumed that (α = 0, s⁺ = 17) for the simulations.

Figure 2 shows the DR losses as function of yaw angle for different riblet geometries at s⁺ = 17. The experimental data of Hage et al. [75] (at Re_b ≈ 14900) and the numerical data of Grüneberger et al. [27] (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 on α is different for each geometry, which might be explained by differences in h/s and riblet shape [75]. Of the two trapezoidal geometries, the one with the sharpest tip (30°) is slightly more sensitive to misalignment. The triangular geometry has the largest tip angle of about 54°, but its DR 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. [27] 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) [75].

5. Herringbone riblets

5.1. Texture description

The feather texture (described in Section 1) is modelled as shown in Figure 3, fairly similar to the spatial three-dimensional (s-3D) texture proposed in Chen et al. [31]. It consists of N_feather feathers placed in parallel such that their edges touch. The modelled feather has no physical shaft, although the term ‘shaft’ is employed to denote the symmetry axis of a feather. The barbs are modelled 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 α, which is restricted to 0 ≤ α < 180°. The shortest distance between the blades is s. Between two neighbouring blades a groove forms, which has spanwise blade spacing s_y = s/cos (α). N_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 Λ_f = 2N_grooves_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 (11) Note that Koeltzsch et al. [41] and Nugroho et al. [42] 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.

Figure 3. Herringbone blade riblet texture with 4 grooves per feather half (N_groove = 4), and angle between shaft and riblets of 15 degrees (α = 15°). Left: Bird's-eye view of the texture, showing five unit cells in the streamwise and three unit cells (or feathers) in the spanwise direction. The edges and shaft of the middle feather are shown. Backward flow over this texture is equivalent to forward flow over a texture with α = 165°. Right: Top view of the texture, showing two unit cells in the streamwise and one unit cell in the spanwise direction. The feather width is Λ_f. Four grooves on one feather half are numbered.

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 α is the same as backward flow over a texture with angle α_backward = 180° − α. Using this trick, textures with angle α and α_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 [53], namely one that extends from the bottom to the top wall. The second variation is a riblet texture with α = 0° (see Figure 4), which results in a parallel blade riblet geometry with blade height variation in the spanwise direction given by Equation (11(11) ). This texture does not suffer from yaw and the resulting pressure drag, but it still might give rise to secondary flows.

Figure 4. Herringbone riblet texture variations. Left: Bird's-eye view of a shifted texture variation with N_groove = 4 and α = 15°. ‘Shifted’ refers to the spanwise shift of the top wall texture by half a spanwise texture wavelength as compared to the bottom wall. The two vertical planes demarcate the middle unit cell. Note that the distance between bottom and top wall is not to scale. Right: Bird's-eye view of a parallel texture variation with N_groove = 4 and α = 0°, showing three unit cells (or feathers) in the spanwise direction.

The herringbone texture and its variations have been subjected to a parametric study (see also Table C1). According to Chen et al. [31], bird-feather parameters are typically s⁺ ≈ 20, h/s ≈ 0.5 and α ≈ 30°. Values for N_groove or Λ_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 α and N_groove. The angle α was either 15° or 165°, which permits the direct comparison with the conventional riblet texture in 15 degrees yaw. Larger angles (e.g. α = 30° or α = 150°) were not considered, as Figure 2 suggests that such angles will result in a drag increase. The number of grooves N_groove was varied between 1 and 128, which resulted in a feather-width change from Λ_f = 0.10 to 12.9 in outer units or Λ⁺_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 C1 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 [79], 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 [55,80], square and triangular elements [56], and a flat plate normal to the free stream [81,82]. 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 α = 165°, 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 α = 15° and N_groove = 1 does not comply with this trend. DR 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 α = 15°. It is significantly less than the 7.6% reduction for parallel riblets aligned with the mean flow.

Figure 5. Drag change as function of feather width for the herringbone riblet geometry, including results for the parallel (α = 0°) and shifted variants (see Figure 4). The shown values of Λ_f and Λ⁺_f apply to the textures with α = 15° or 165°, and are approximate for α = 0° textures. The data points on the right vertical axis belong to conventional parallel riblets with yaw angle α = 0° and α = 15°.

Difference forward and backward flow: There is a clear drag difference between textures with α = 15° (forward flow) and α = 165° (backward flow). For example, at N_groove = 16 drag increases by 34% for α = 15°, and only by 7.5% for α = 165°. This shows that mean-flow reversal has a considerable effect on drag, in agreement with Chen et al. [31]. Here, forward flow experiences more drag than backward flow, except for the N_groove = 1 texture.

The drag difference between forward and backward flow decreases for increasingly wide feathers, as can be expected. The textures with α = 15° and 165° differ only because of riblet convergence or divergence in the feather shaft or feather edge regions. When Λ_f → ∞, 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 (α = 0°) 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 DR is the least for N_groove = 1, namely 0.5%. The texture with N_groove = 1 and α = 0° 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 DR. When N_groove increases, DR increases as well. In the limit of very large N_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 DR to be close to that for standard riblets.

5.3. Flow description

Instantaneous streamwise velocity: To understand the DR 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 6. Three cross sections of the instantaneous streamwise velocity in the simulation domain that belongs to the herringbone texture with N_groove = 4, α = 165°. The horizontal plane is located at z = 0.016 or z/h = 0.68.

Figure 7 shows part of the mean streamwise-averaged flow fields for textures with α = 165° and N_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.

Figure 7. Streamwise-averaged mean flow in a plane perpendicular to the streamwise direction for herringbone textures with α = 165° and N_groove = 1 (top left), 4 (top centre and right), 16 (bottom). Vectors exhibit in-plane secondary flow. Contours represent streamwise velocity (top left, top centre, bottom) or streamwise velocity fluctuations (top right, ).

: For textures with α = 165°, 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 centreline, probably because of the small misalignment of the riblets with the mean flow (α = 165° is relatively close to 180°). For N_groove ≥ 4, it is limited in vertical extent to approximately z = 0.2, as can be seen for N_groove = 4 and 16 in the figure. Because of spanwise confinement, for small feather widths the vortex decreases in size. For N_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 [52].

While mainly one vortex constitutes the secondary flow for N_groove = 1 and 4, tertiary flows appear for larger feather widths. For example, the mean flow for N_groove = 16 (see Figure 7) shows a counterclockwise-rotating flow that extends to the channel centreline. Less pronounced is the small clockwise-rotating vortex near the feather edge. Several tertiary flows were also observed for even wider feathers (N_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 [52]. These tertiary flows are likely similar to the secondary flows that form over streamwise-aligned roughness strips [83]. 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).

: Reversing the flow direction from backward to forward also changes the secondary flow direction. For textures with α = 15°, 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_groove ≥ 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 centreline, it has no noticeable effect on the mean flow in the bottom channel half. For N_groove = 1 and 4, this spanwise shift is the only change of the mean flow. For N_groove = 16, however, the pronounced counterclockwise-rotating tertiary flow is also modified. It reaches the centreline 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_groove = 16 the feather half-width (Λ_f/2 = 0.80) approximately equals the channel height (L_z = 1). Tertiary flow for N_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 characterised 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 [41,42] and in DNSs of turbulent boundary layers with uniform blowing/suction [3].

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. [84]). The contours represent the shear stress τ⁺_w on the bottom wall.

Figure 8. Instantaneous vortical structures and wall shear in a part of the bottom channel half. Vortical structures are iso-surfaces of the second invariant of the velocity gradient tensor with Q⁺ = 0.03. The coloured contours represent the wall shear stress τ⁺_w on the bottom wall. Three cases are shown (from top to bottom): smooth walls; parallel riblets in yaw of α = 15°; herringbone riblets with N_groove = 4, α = 165°. Re_b = 5500 for all cases and s⁺ = 17 for all textures.

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 [3].

5.4. Secondary flow strength

To quantify the strengths of the secondary flows described in the previous subsection, the secondary flow strength Γ_s is introduced: (12) with and . The secondary flow strength is decomposed into the mean-secondary-flow strength Γ_m and the turbulent-secondary-flow strength Γ_t. These variables still depend on the spatial coordinates. A volume- and xz-average of Γ_m and Γ_t is shown in Figure 9.

Figure 9. Strength of mean and turbulent secondary flow for the herringbone riblet geometry. Top left: Volume-averaged mean-secondary-flow-strength as function of feather width (specified by N_groove). The data points on the right vertical axis represent the conventional parallel blade riblets with yaw angle α = 0° and α = 15° (not visible with ⟨Γ_m⟩_V = 0.027). Top right: Streamwise- and wall-normal-averaged mean-secondary-flow-strength as function of spanwise distance for textures with α = 15° (left part) and α = 165° (right part). Bottom left: turbulent-secondary-flow-strength plotted as in top-left subfigure. Bottom right: turbulent-secondary-flow-strength plotted as in top-right subfigure. In the bottom figures, the turbulent-secondary-flow-strength of the smooth wall is subtracted.

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 Λ_f). For all N_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² + w²) projects both onto Γ_m and Γ_t (see Equation (12(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_groove ≥ 4. The results for N_groove = 1, however, deviate in this respect: ⟨Γ_m⟩_V more than halves as compared to N_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 Γ_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 α = 165°, textures with α = 15° 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 Γ_t varies with spanwise distance. The relatively uniform turbulence for the textures with N_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 α = 15° and increases when α = 165°. 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_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_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 α = 0° is much weaker than for the herringbone riblet geometries, which is attributed to alignment of the riblets with the mean flow. Both ⟨Γ_m⟩_V and ⟨Γ_t⟩_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_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(1) )) is rewritten as ∂u/∂t = ∑f, where f is one of the terms in that equation. A time, streamwise and wall-normal average (denoted by ) is applied to this equation. Assuming statistical stationarity and using the boundary conditions, the balance becomes: (13) 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^tot_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 Δz_w from 0 at the feather edges to h at the shaft.

Figure 10. Streamwise momentum balance as function of the spanwise coordinate (Equation (13(13) )). Six textures are shown: narrow herringbone feathers (N_groove = 4) for α = 15° (top left) and α = 165° (bottom left); wide herringbone feathers (N_groove = 128) for α = 15° (top centre) and α = 165° (bottom centre); parallel riblets aligned with the mean flow (top right) and in yaw of α = 15° (bottom right). The blade location is denoted by y_b. For all cases, Re_b = 5500 and s⁺ = 17.

Before the differences between N_groove = 4 and 128 are highlighted, the general behaviour of the different terms is clarified. The driving term is almost constant with y, because γ_f changes only marginally from the feather shaft to the edges. For N_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 α = 165°. However, the shear magnitude peaks near the shaft for α = 15°, which is ascribed to the local downdraft that transports high momentum towards the wall. For very wide feathers (N_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 α = 165°. In contrast, it is very large there for α = 15°, likely due to the downdraft of high-speed fluid.

The balance for N_groove = 4 shows significant contributions from spanwise advective and diffusive transport, especially near the shaft. For α = 15°, 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 centreline). In contrast, diffusion transports high momentum away from the shaft. These trends are opposite to that of textures with α = 165°, for which spanwise advection causes a streamwise momentum decrease and diffusion an increase near the shaft.

Compared to N_groove = 4, the relative importance of the terms in Equation (13(13) ) is very different for N_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 ) is applied to the Navier–Stokes equation ∂u/∂t = ∑f. Assuming statistical stationarity and using the boundary conditions, the balance reads: (14) The terms in this equation only depend on the wall-normal coordinate z. All terms are divided by d^tot_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 α = 15° (not shown) are similar to the ones for α = 165°. The behaviour of the different terms is clarified below. The driving term is almost constant, as before.

Figure 11. Streamwise momentum balance as function of the wall-normal coordinate (Equation (14(14) )). Four cases are shown (from left to right): smooth walls; parallel riblets in yaw of α = 15°; narrow herringbone feathers (N_groove = 4) for α = 165°; wide herringbone feathers (N_groove = 128) for α = 165°. Re_b = 5500 for all cases and s⁺ = 17 for all textures.

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 0 and h (see Equation (11(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 centreline, 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_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_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. [85] derived an equation (called Fukagata-Iwamoto-Kasagi (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): (15) where it is used that the domain height L_z = 1. Note that the skin-friction coefficient for smooth-wall channel flow.

This relation has been extended to geometrically more complex surfaces by Peet and Sagaut [86]. 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(14) )). For ease of notation, that equation is written as 0 = ∑_lf^(l)|_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(14) ): 0 = ∑_l∫¹₀.f^(l)|_z dz or (16) Next, Equation (14(14) ) is rewritten as 0 = ∑_l(1/2)(f^(l)|_z + f^(l)|_{1 − z}) 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: (17) Using the boundary conditions at the channel walls, the definition of the bulk velocity (i.e. ), integration by parts to transform multiple to single integrations and the global drag balance (i.e. Equation (16(16) )), the last equation becomes: (18) To arrive at an equation for , Equation (18(18) ) is divided by the prefactor in square brackets and multiplied by ∫¹₀⟨γ_f⟩_xy dz, which yields the final extended FIK-identity: (19) Notice that the total advection term is split up into contributions from mean and turbulent flow, using that . The symmetry or antisymmetry with respect to the channel centreline at z = 1/2 is evident from the factors in braces in Equation (18(18) ). As the multiplier of in that equation depends on γ_f, the bulk term d^bulk is slightly texture dependent, as was also found by Peet and Sagaut [86]. For smooth walls, γ_f = 1 and that multiplier equals 1/12, such that Equation (15(15) ) is recovered.

As Equation (19(19) ) applies to both smooth and textured walls, the drag change can be decomposed in a similar way: (20) 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 DR or increase.

The identity presented here is somewhat different from the one derived by Peet and Sagaut [86]. First, their decomposition only applies to skin friction, whereas the drag decomposition in Equation (19(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 [86]. 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(19) )) is not only useful when an IBM is used. In the present work, f^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 3D phase-indicator function γ_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 3D body force or IBM force.

In view of the factors between braces in Equation (18(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 centreline. 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 equalises the left-hand and right-hand side of Equation (19(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.

Table 2. Contributions to the total skin friction for smooth-wall turbulent channel flow. The terms are as given in Equation (19(19) ) and divided by d^tot. The rest term equalises the left-hand and right-hand side of that equation.

Reb	d^tot (%)	d^bulk (%)	d^mean advection (%)	d^turbulent advection (%)	d^IBM (%)	d^rest (%)
5500	100	26.9	4.3 · 10^−17	73.1	0	−0.013
11000	100	16.0	1.7 · 10^−17	83.9	0	0.086
22000	100	9.6	1.4· 10^−16	90.3	0	0.10

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^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 [86]. 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(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 α 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.

Figure 12. Different contributions to the total drag change (Equation (20(20) )) for parallel blade riblets. Note that five bars belong to only one abscissa. Left: Decomposition as function of riblet spacing in wall units for riblets aligned with the mean flow (α = 0°). Right: Decomposition as function of yaw angle for fixed riblet spacing s⁺ = 17.

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. [61,87]). 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 α, 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 α 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 α = 0°, 8.6% for α = 15°). Advection by the mean flow can be a significant contribution to drag increase, up to 28% for N_groove = 4, α = 165°. The turbulent advection term exhibits the largest changes: it varies between a 6% reduction and a 50% increase.

Figure 13. Decomposition of the total drag change (top left) for herringbone riblets in contributions from the IBM-force (top right), turbulent advection (bottom left) and mean advection (bottom right) according to Equation (20(20) ). All figures show the contribution change relative to a smooth wall as function of feather width. The data points on the right vertical axes represent the conventional parallel blade riblets with yaw angle α = 0° and α = 15°.

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 Γ²_s = Γ_m² + Γ²_t and , respectively. In summary, the results clearly indicate an intimate connection between increased advective transport and stronger secondary flow.

Figure 14. Advective transport in a plane perpendicular to the streamwise direction for the herringbone texture with N_groove = 4, α = 165°. Contours represent the instantaneous advection uw (left), the averaged mean advection (centre) or the averaged turbulent advection (right). Vectors depict the instantaneous (left) or averaged (centre, right) in-plane flow velocity.

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_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 centreline does not contribute much to drag because of the factor (1/2 − z) in Equation (18(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 centreline is weak.

Spanwise decomposition of drag change: Figure 13 shows that DRs are possible in the limit of large N_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 be one term in Equation (14(14) ) and d^(l) the corresponding drag contribution. The calculation of this FIK-term can be rewritten as: (21) The function g = g(z) results from conversion of a triple to a single integral, and normalisation. The function d^(l)_sp = 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: (22) where dc^(l)_sp quantifies how the drag change depends on the spanwise coordinate. Note that d^(l)_sp and dc^(l)_sp have a spanwise dependence, whereas d_smooth^(l) and d^tot_smooth have not. The spanwise decomposition of the mean and turbulent advective terms is shown in Figure 15 for α = 15° and 165°, 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_groove = 128 probably appear for the same reason.

Figure 15. Spanwise decomposition of advective contributions to the total drag change (Equation (22(22) )). The total advection term is split into a mean (top) and turbulent (bottom) part. Results for herringbone textures with α = 15° (left) and α = 165° (right) are shown. The number in the legend quantifies dc, the spanwise average of each curve (Equation (22(22) )), which was shown already in Figure 13.

The effect of convergent/divergent riblets can be isolated from that of parallel riblets in yaw for the textures with N_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 α = 15° or 165°. These favourable 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 α = 165° and any N_groove, the updraft of fluid around the shaft results in a local DR, 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_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_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 α = 15°.

Local turbulent DR can be obtained by diverging riblets. For α = 15° and N_groove ≥ 16, the strong wall-directed mean flow is accompanied by reduced turbulent transport. The reverse is true for α = 165°, 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 DNSs, 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 [3].

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 Figures 14 and 15. However, the volume-averaged turbulent term is often at least as important as the mean term. As a second example, the texture with α = 165° and N_groove = 4 exhibits locally a DR as high as 400% due to a strong wall-leaving flow. However, one should realise that the favourable 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 α = 165° and N_groove = 4, the vortex covers one feather half. The corresponding spanwise-averaged mean advective drag is 28%, which establishes again an unfavourable 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 DR.

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 Λ_f/L_z ≲ 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 DR, 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 Λ_f → ∞ 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 unfavourable 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 travelling wave of spanwise body force is described by f_y = A(z)sin {(2π/Λ_f)y − ωt}. Here, f_y represents the spanwise body force, A(z) the forcing amplitude (only nonzero in the vicinity of the wall), Λ_f the spanwise wavelength of the forcing, and ω 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 (ω = 0). Research shows that the best-performing spanwise travelling wave (in terms of DR or net energy saving) is the one with infinite wavelength, i.e. the spanwise wall oscillation [88,89]. In other words, Λ_f → ∞ 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 travelling waves of spanwise forcing are superior to spanwise travelling waves of spanwise forcing, presumably because of the unaltered wall-normal gradient of the wall-normal velocity component at the wall [89]. Both experimentally and numerically, DR has been obtained with a stationary streamwise variation of spanwise forcing [27,90]. 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 [91,92]. 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 unfavourable 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 DR apparently contrasts with the experimental study of Chen et al. [31]. 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. [31] obtained DR in forward flow for Λ_f/D ≲ 1 (D being the pipe diameter), with a maximum of 20%. In contrast, the present study only achieved DR for Λ_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_groove was used in the experiments. The experimental texture consisted of sawtooth riblets at an angle of 30° with the flow direction. The numerical texture was composed of blades at an angle of 15°. 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 ≈ 2.6 · 10⁵ (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 [44–47]) 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. Λ_f) and the DR 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 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 DR of 9.3% was obtained, close to the 9.9% that has been found experimentally [21].

The herringbone texture can both increase or reduce the drag, depending on the spanwise texture wavelength Λ_f. For Λ_f/L_z ≲ O(1) (i.e. narrow feathers), the drag increases with a maximum of 73% for Λ_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 centred 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 DR of 2% was found for Λ_f/L_z ≳ 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 travelling waves of spanwise forcing, the current study confirms that Λ_f → ∞ 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 DR, which apparently contrasts with the experiments of Chen et al. [31]. 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 (Λ_f, s, h/s, α) 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 DR with the texture and flow parameters is presumably non-trivial.

Whether feather riblets have an aerodynamic function remains an open question. Chen et al. [31] claimed that feather riblets greatly impact flight performance because of DR. The DR was ascribed to suppression of turbulent momentum transport, but evidence has been inconclusive so far. One should realise 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 [35–37]. Furthermore, separation delay has been obtained with vortex generators that resemble the herringbone texture [33]. So, future studies might investigate the potential of the herringbone texture for flow-separation control.

Acknowledgments

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.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was sponsored by the European Union Seventh Framework Programme in the SEAFRONT project [grant agreement 614034]; and the Netherlands Organisation for Scientific Research (NWO) [project number SH-023-15].

Appendix A. Time advancement at fixed bulk velocity

For simplicity, Equation (1(1) ) is rewritten as: (A1) 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^(s)_i, where u⁽⁰⁾_i = u_iⁿ and u⁽³⁾_i = u_i^{n + 1}. Similarly, intermediate pressures p^(s) and driving forces f^driving(s) are introduced. Because of the Crank–Nicolson scheme for pressure, with correction pressure . The time advancement is illustrated here for an arbitrary sub-step s: (A2) where u*^(s)_i is the first and u*^*(s)_i the second prediction velocity. Note that the asterisk here is not used to denote dimensional quantities. The parameters α_(s), γ_(s) and ζ_{(s − 1)} are RK3 parameters (see e.g. [67]). To obtain the driving force, the equation for the streamwise velocity is volume-averaged, which yields: (A3) 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^driving(s), which yields: (A4) Next, u*^*(s)_i is computed. With use of the continuity equation, the divergence of Equation (A2(A2) ) yields a Poisson equation for the correction pressure: (A5) When this is solved, the updated velocity and pressure are computed: (A6) 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: (A7)

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^IBM_i term in Equation (1(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 B1 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 dash-dotted boxes in that figure indicate three staggered grid cells that belong to the velocity components u_left, v_{(j, k)} and w_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.

Figure B1. Numerical grid and IBM used for parallel riblets without yaw. Left: Part of the numerical grid used for s⁺ = 17 simulations at Re_b = 5500, showing two unit cells in the spanwise direction. Both for blade spacing and blade height, 17 grid cells were used. Right: Small part of the numerical grid (8 grid cells) around the blade tip. The text explains the IBM with use of this subfigure. The three dash-dotted boxes represent staggered grid cells that belong to the velocity components shown in their centres. Each of these cells has one face indicated with a thick and solid line. At that cell face, the IBM adjusts the advective and diffusive fluxes.

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_left, the two spanwise fluxes through its thick solid face were changed: the advective flux uv = 0 and the diffusive flux ∂u/∂y = −2u_left/Δ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_blade represent a spanwise velocity component that coincides with a blade (see Figure B1) and let ^* represent the first prediction velocity (as in Equation (A2(A2) )), then . The actual velocity v_blade follows from the correction step (Equation (A6(A6) )) 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 as a dash-dotted box in Figure B1). The vertical fluxes of v at the bottom face of that cell were adjusted. The diffusive flux ∂v/∂z was split into two contributions, namely from the left and right side of the blade (indicated by φ_left and φ_right in the figure). It accounts for the fact that the thin blade does not inhibit vertical transport. Specifically, ∂v/∂z = 0.5φ_left + 0.5φ_right. The flux φ_right was computed by linear interpolation of the four velocity components labelled with indices in the figure: φ_right = (v_t − v_b)/Δ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 φ_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_left, two spanwise fluxes at the thick solid face were changed: the advective flux vw = 0 and the diffusive flux ∂w/∂y = −2w_left/Δ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 B2 (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 α. 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 B2 has N_cg = 16.

Figure B2. Numerical grid and IBM used for parallel riblets in yaw. Left: Top view of the texture with α = 20°, showing one unit cell in the streamwise and two in the spanwise direction. The numerical grid consists of N_cg = 16 grid cells per groove. The markers on one blade indicate that the blades intersect the grid cells at the locations of the staggered velocity vectors u (circle) and v (square). Right: Small part of the numerical grid (8 grid cells) around the blade tip. The text explains the IBM with use of this subfigure, see also Figure B1.

The IBM is explained with the use of Figure B2 (right). Staggered grid cells for u_top and w_left are shown as dash-dotted 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. ) at the grid points that coincide with a blade. At the bottom cell face of the grid cell for u_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: ∂u/∂z = 2u_top/Δ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_left were adjusted as follows: uw = 0 and ∂w/∂x = −2w_left/Δx.

B.3. Computation of f^IBM_i

Although the IBM forcing f^IBM_i appears as a separate term in Equation (1(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*^(s)_i 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 in absence of the IBM forcing. The intermediate IBM forcing f^IBM_i(s) then follows from . Finally, the total IBM forcing f^IBM_i is obtained from summing α_(s)f^IBM_i(s), similar to Equation (A7(A7) ).

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_blade and v_blade denote two velocity components that coincide with the texture, then the corresponding penetration velocities are and . Note that w_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 ⟨u_pen⟩, ⟨v_pen⟩ and max (u_pen, v_pen), where ⟨...⟩ denotes an average. The resulting time series were time averaged, yielding , and . The magnitude of these penetration velocities (normalised with the bulk velocity) is similar for all the simulations, namely a mean on the order of 10⁻⁶ to 10⁻⁵ and a maximum on the order of 10⁻⁴ to 10⁻³. This shows that the penetration velocities are very small, so the no-slip and no-penetration conditions are sufficiently satisfied.

Figure B3 presents a streamwise velocity profile as function of the spanwise distance for a herringbone riblet texture with four grooves per feather half (so N_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 eight riblet locations. This illustrates that the boundary conditions at the riblet surfaces are satisfied.

Figure B3. Spanwise profile of streamwise velocity. Left: Bird's-eye view of one unit cell of a herringbone texture with N_groove = 4, α = 165°. The thick line parallel to the y-axis cuts through the centre of the first grid cell above the bottom wall. Right: Time-averaged streamwise velocity as function of spanwise distance along the thick line in the left subfigure.

Appendix C. Simulation parameters

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

Table C1. 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 refers to the conventional parallel-riblet texture. The averaging time T is normalised with δ/u_τ derived from smooth-wall flow at the same bulk Reynolds number.

Short name	s^+	α (°)	h/s	Ngroove	Lx	Ly	Nx	Ny	Nz	Δx^+	Δy^+	Δz^+w	Δz^+c	Ncg	Reb	Tuτ/δ	d^tot · 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 ± 0.02	0.0 ± 0.4
Re11000	-	-	-	-	4.9	2.5	800	800	512	4.0	2.0	0.50	2.0	-	11000	66	6.81 ± 0.02	0.0 ± 0.5
Re22000	-	-	-	-	5.4	2.7	1600	1600	1024	3.9	2.0	0.49	1.8	-	22000	14	5.70 ± 0.02	0.0 ± 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 ± 0.04	−6.9 ± 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 ± 0.05	−7.6 ± 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 ± 0.04	3.0 ± 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 ± 0.04	3.4 ± 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 ± 0.03	−9.3 ± 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 ± 0.04	−0.5 ± 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 ± 0.04	−1.6 ± 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 ± 0.05	−4.0 ± 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 ± 0.03	−0.9 ± 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 ± 0.04	−1.2 ± 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 ± 0.04	3.7 ± 0.6
Herringbone riblets:
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 ± 0.06	53.8 ± 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 ± 0.1	73.4 ± 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 ± 0.09	75.7 ± 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 ± 0.05	34.1 ± 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 ± 0.02	2.3 ± 0.4
Herringbone riblets:
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 ± 0.1	61.2 ± 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 ± 0.07	47.9 ± 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 ± 0.04	7.5 ± 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 ± 0.04	2.1 ± 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 ± 0.02	−2.0 ± 0.4
Herringbone riblets: , 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 ± 0.1	60.3 ± 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 ± 0.1	47.2 ± 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 ± 0.03	8.5 ± 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 ± 0.05	2.1 ± 0.7
Herringbone riblets: , 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 ± 0.03	−7.6 ± 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 ± 0.03	−0.5 ± 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 ± 0.04	−4.1 ± 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 ± 0.08	−5.7 ± 1.0