Drag reduction is of great importance because of the environmental and economic benefits from the reduced fuel consumption [1]. As one of the oldest and most investigated drag reduction methods, riblets, small protruding surfaces along the direction of the flow, had been known to be capable of reducing friction drag up to 8% with appropriate height and spacing in spite of the extra surface area. Early experiments conducted by Walsh et al. [2, 3], Bacher et al. [4], and Bechert et al. [5] tested different shapes including triangular, rectangular (or blade), trapezoidal, semi-circular, and other shape grooves with various sizes to get an optimal drag reduction rate. Experiments on concave and convex riblet shapes showed that drag reduction increases as the peak curvature increases, or as the radius of valley curvature increases, which means that an optimum riblet shape should have sharp peaks and curved valleys [2]. However, later investigations put more emphasis on the riblet tip rather than the valley, arguing that the tip sharpness plays an important role in damping the spanwise velocity fluctuations and thus limiting the momentum transport and turbulence intensity.

Two drag reduction regimes of riblets had been recognized [5], i.e., the “viscous regime,” W + < 10 − 15 , where the drag reduction rate increases almost linearly with the riblet width; and the “breakdown regime,” W + > 15 − 20 , where the drag reduction rate reaches a maximum and then decreases even to the point of drag increase. Here, W + represent the riblet width in wall units and the protruding surfaces would be referred to as roughness when entering the drag increase state. The physical mechanisms for drag reduction by riblets have been extensively studied, for example Walsh [6], Choi et al. [7], Choi [8], and Rstegari et al. [9]. Researchers often owe the credit to the interaction between the riblet surfaces and the longitudinal vortices of the turbulent boundary layer. One widely accepted explanation for riblets of the viscous regime, or for infinitesimal riblets, relies heavily on the concept of protrusion height [10, 11]. Intuitively, the protrusion height is understood as the distance between the riblet tip plane and the origin of the velocity profile, which lies between the riblet tip and riblet valley. The idea behind the emphasis of the protrusion height concept is that the influence of riblets is mainly constrained to within the inner layer and is thus viewed as an inner layer control strategy, and what the outer layer can “feel” is just the shear between the inner and outer layer. Therefore, the protrusion height denotes the location of a hypothetical flat plane that the outer layer perceives. The true protrusion height definition was given by Luchini et al. [11] From the intuitive definition above, we can easily define a longitudinal protrusion height h ∥ and a perpendicular protrusion height h ⊥ based on average streamwise and spanwise velocity respectively. Then the difference ∆ h = h ∥ − h ⊥ is referred to as the protrusion height. The basic idea is that compared with streamwise flow, ∆ h represents the extent at which the cross-flow fluctuations are hampered, which had been viewed as critical for the turbulence generation cycle. In the viscous regime, the drag reduction rate is linearly dependent on this protrusion height ∆ h , a fact that has been verified by experiments and numerical simulations for various riblet shapes. However, the mechanism for the deterioration of drag reduction, or the “breakdown regime” remains controversial. The generation of secondary streamwise vortices by unsteady crossflow separation, which brings high speed flow towards the riblet surfaces, was suggested by Goldstein and Tuan [12], but secondary streamwise vortices do not necessarily lead to drag increase. For example, spanwise oscillation of the wall can also decrease drag [13, 14]. The second group of theories emphasized the scale interaction between riblet surfaces and turbulent coherent structures, i.e., near-wall streaks and streamwise vortices [7, 15]. Generally, the optimal lateral scale of riblets is an order of magnitude smaller than the spacing between low-speed streaks (about 100 in wall units). It was argued that for riblets larger than W + ≈ 20 − 30 , the near wall streamwise vortices could move freely and lodge inside the riblet valleys, which leads to an increase of wall surface area exposed to the sweep motion induced by those vortices. On the other hand, for smaller riblets, the streamwise vortices would stay above the riblet tips, and only the tip region could be subjected to the induced sweeps. Alternatively, Garcia-Mayoral and Jimenez [16] attributed the breakdown to the appearance of large-scale spanwise rollers. The mechanisms suggest plausible reasons based on experimental and numerical observations, but lack substantial support of improved design of riblets based on these theories.

As a bio-inspired technique mimicking the denticles of fast swim sharks, riblets are possibly the only turbulent drag reduction strategy that have been tested in application. Szodruch [17] reported that covering of riblets on 70% of the surface on an Airbus 320 commercial airplane leads to a 2% reduction in oil consumption. Instances of using riblets in sports, for example, on the hulls of boats and the surfaces of racing swimsuits, have been successful [18]. In contrast to the canonical arrangement, novel riblet concepts have been introduced. Nugroho et al. [19] explored the possibility of using ordered and directional surfaces to redirect the near wall flow. Benhamza et al. numerically investigated variable spacing riblets of rectangular shape in turbulent channel flows. Riblets with a sinusoidal variation along the streamwise direction [20, 21] had also been devised to mimic the spanwise oscillation strategy to reduce drag, however, these gained very little further drag reduction. Boomsma [22] numerically studied denticles resembles sharkskin in turbulent boundary layer, however, obtained drag increase instead of drag reduction. These expeditions generally failed to surpass the optimum drag reduction rate of conventional riblets, indicating a lack of full understanding of the drag reduction mechanism.

The plausible mechanisms of breakdown imply a requirement of tip sharpness of riblet for high drag reduction performance. However, sharp tips pose challenges for manufacturing and maintenance. Actually, Walsh [2] experimentally found that the optimum riblet appears to have sharp peaks and significant valley curvature. Launder and Li [23] also reported that the U-form (or scalloped) riblets can achieve a superior performance to the corresponding V-shaped riblets with the same height and width. Based on the protrusion height concept and numerical simulation of riblet controlled boundary layer transition, Wang et al. [24] showed that scalloped riblets could have better performance than corresponding triangular riblets with sharper tips. The scalloped riblet shape we considered is constructed by smoothly connecting two third-order polynomials. y γ W = 2 a z W 3 − 2 A z W 2 + 1 , 0 ≤ z ≤ z * 2 b z W − 1 2 3 + 2 B z W − 1 2 2 , z * < z ≤ W 2 (1) where W is the width of the riblet, and the shape of the other half riblet ( W / 2 < z ≤ W ) can be obtained by symmetry about z = W / 2 . Note, we have selected a coordinate system consistent with the following three-dimensional simulations. x , y , and z directions are streamwise, wall-normal, and spanwise directions respectively. Being aligned in the streamwise direction, x is not present in the riblet shape definition. γ is the height-width ratio while A and B are chosen parameters to determine the curvature at the tip and valley of the riblets. Note that an extra factor of 2 is introduced to make A and B the coefficients of second order terms when the width of the riblet is S = 1 and height-width ratio is γ = 0.5 . a , b , and z * are obtained by requiring the two polynomials smoothly connected, i.e., the functions value, first and second derivatives are continuous at z * . Despite the non-linearity in z * , the equations can be solved analytically z * W = 6 − B 2 A − B , a = 2 A + B z * / W + A − B 3 z * / W , b = 2 A + B z * / W − 2 B 3 z * / W − 1 / 2 (2)

To guarantee that z * / W is in the range of 0,0.5 , we have to ensure that the given A and B satisfies 0 ≤ B < 6 < A , or   0 ≤ A < 6 < B (3)

An extreme exists with A = B = 6 , when the two polynomials are identical and the curvatures at the tip and valley are perfectly balanced. Thus, Eq. 1 represents a system of scalloped riblets, the tip and valley curvatures of which can be easily defined. Wang et al. [24] calculated the protrusion heights for scalloped riblets with various A , B , and γ values with a boundary element algorithm proposed by Luchini et al. [11], and found that for a larger height-width ratio γ , the protrusion height is mainly determined by A , i.e., the curvature at riblet tip, while for a smaller height-width ratio γ , the protrusion height is mostly dependent on B , i.e., the curvature at the riblet valley. However, for riblets in practical applications, the height-width ratio is generally in the range of 0.5 ≤ γ ≤ 1.2 , the protrusion height is influenced equally by A and B . Therefore, the shape of the riblet valley cannot be overlooked, and should be designed more carefully. It is also illustrated by direct numerical simulations of a boundary layer transition that a scalloped riblet, which is not as sharp in the tip as a corresponding triangular riblet with the same height-width ratio, nevertheless has a larger drag reduction rate, which is in accordance with protrusion height results. Another advantage of the scalloped riblet over the triangular riblet is its stronger resilience to riblet tip erosion, which means better manufacturing and maintenance performances. For the current study, we will focus on the same scalloped riblet as in Wang et al. [24] with A = 25 , B = 1.6 and γ = 0.5 . It was obtained, to four decimal digits, that a = 100.6970 , b = − 1.4793 and z * / W = 0.0940 . Note that similar groove shapes had been studied by Klumpp et al. [25] for boundary layer transition but without detailed information of the groove shape. In addition, there are two methodologies to resolve the riblet shapes in numerical simulations, which involve the utilization of either body-fitted meshes or an immersed boundary method. With body-fitted meshes, the riblet shapes are better represented, but lack the versability of representing different shapes and sharp corners, especially when high-order numerical methods are adopted. On the other hand, the immerse boundary method can easily be used to simulate different riblet shapes even with the same mesh, which, however, has been criticized for low accuracy near the wall surfaces. For the current study, we choose the immersed boundary method for its versability and adopt Lagrange interpolation [26] to improve the accuracy near the wall surface.

One important aspect in understanding the flow dynamics is the identification of vortices. Noted by Küchemann [27], vortices are the “sinews and muscles of the flow.” However, no consensus on the definition of vortices has been reached. Popular methods including Q , Δ , λ 2 , λ c i , Ω [28–32] and other criteria have been introduced to visualize vortices in various flows. However, these methods mainly suffer from two issues: 1) These methods are different with different physical reasoning and dimension units among each other; 2) as scalar-valued methods, no directional information about the rotational motion is provided. To overcome these issues, a new vortex vector, Liutex vector [33, 34], was introduced with its magnitude as twice the angular velocity of the rigid rotation part of the fluid motion and its direction as the local rotational axis. An explicit expression for the Liutex vector R in terms of vorticity ω , eigenvalues and eigenvectors of the velocity gradient tensor is proposed by Wang et al. [35] which leads to efficiency improvement and physical intuitive comprehension as R = R r = ω ∙ r − ω ∙ r 2 − λ c i 2 r (4) where λ c i is the imaginary part of the complex conjugate eigenvalues and r is the real eigenvector of velocity gradient tensor. Thereafter, multiple vortex identification methodologies have been developed including the Liutex iso-surfaces, Liutex- Ω method [36, 37], Liutex core line method [38, 39], objective Liutex method [40], etc. In addition, a particular feature of the Liutex vector is that it is free from shear contamination. Kolář and Šístek [41] found that only the Liutex vector is not contaminated by stretching and shear in that an arbitrary adding or subtracting of stretching and shear would not alter the resulting vortex identification. This distinctive feature has been utilized by Ding et al. [42] to develop a new Liutex-based sub-grid model for large eddy simulation, which has been proved to outperform the famous Smagorinsky model in homogeneous isotropic turbulence and turbulent channel flows. In the current study, the Liutex vortex identification method will be adopted to analyze the flow field of riblet controlled turbulent channels.

The paper is organized as follows. Section Numerical Methods and Case Setup introduced the numerical methods adopted, especially the customized immersed boundary method that we use to model riblet surfaces. Section Numerical Results presents the numerical results, discussing the drag reduction rate, mean flow and turbulent statistics, premultiplied power spectrum density, and instantaneous flow field, with particular attention paid to the Liutex field. Finally, conclusions are drawn in Section Conclusion.

The canonical case of a turbulent channel flow at Reynolds number R e τ = 180 was considered and used to study the drag reduction capacity of the scalloped riblets. The open-source finite difference software “Incompact3d” [43, 44] was adopted to simulate the flat plate and riblet channels, with which the incompressible Navier-Stokes equations were numerically solved. Periodic conditions were naturally employed in the streamwise ( x ) and spanwise ( z ) direction, whilst no-slip conditions were applied in the upper and lower walls in the wall normal direction ( y ). For the riblet case, the no-slip condition was implemented with the immersed boundary method (IBM) which will be elaborated in the following. Sixth order compact schemes were adopted for spatial derivatives and a low-storage third-order Runge-Kutta scheme was used for time advancement. The pressure Poisson equation was solved with a pseudo-spectral method to ensure divergence-free condition up to machinery accuracy. The lengths for the computational domain were selected to be L x = 4 h , L y = 2 h and L z = 2 h with h as the half channel height. The domain is generally smaller compared to other direct numerical simulations of turbulent channel at R e τ = 180 , but still sufficient to capture all the length scales and meanwhile limits the overall mesh point number to accelerate the simulations. A schematic of the flow with riblets is shown in Figure 1A. Both meshes for the two cases were identical and include 256 × 129 × 360 grid points in the streamwise, wall normal, and spanwise directions respectively. Uniform grids were used in the streamwise and spanwise directions, whilst a stretching grid was used in the normal direction with points clustered towards the wall surfaces. In wall units, the intervals in the streamwise and spanwise directions were ∆ x + ≈ 2.81 and ∆ z + ≈ 1.00 while the first interval adjacent to the wall was ∆ y + ≈ 0.98 . Normally, this resolution is more than sufficient to resolve all the scales in channel flow, especially in the spanwise direction. With 18 riblets arranged spanwise, 20 mesh points were used to resolve one riblet, and the riblet width in wall units was W + = 20 , which was around the optimum width. The distribution of mesh points near one riblet surface is shown in Figure 1B. The flow rate was kept constant every time step and the bulk velocity u b was thus kept as 2 / 3 with a bulk Reynolds number of 2800. Note that the bulk Reynolds number was exact while the friction Reynolds number was only approximately 180 and the wall unit lengths in the current study were all obtained based on R e τ = 180 unless otherwise stated.

The riblet surfaces were modeled with a customized immersed boundary method in “Incompact3d” based on an alternating direction forcing to ensure a no-slip boundary condition at the wall of the solid body [26]. A particular treatment was that no-zero velocities were set inside the solid body when computing derivatives to avoid discontinuities of the velocity field, which, of course, would not be used when outputting numerical results. Because the technique is important for both correctly resolving the riblet wall surfaces and precisely calculating the friction drag, which is certainly critical especially when aiming for about 5∼10 percent drag reduction, the procedure is described in more detail in the following.

For the adopted compact finite difference scheme, the derivatives along a directional line, which could cover both fluid and solid regions, were calculated simultaneously in a so-called implicit manner by solving systems of linear algebraic equations. If the condition of velocity component u i = 0 was to be set everywhere inside the solid domain, then the spatial derivative at the interfaces between fluid and solid would not be ensured to be continuous, which leads to deterioration of the solution. The problem only has a minor impact for low-order numerical scheme. When combined with high-order schemes, however, spurious oscillations would be generated at the interfaces. To deal with this issue, non-zero values inside the solid domain were set based on Lagrange or cubic spline interpolation from points on and adjacent to the interfaces. This described 1D procedure was dealt with sequentially for the three dimensions, and hence the name “alternating direction.” For the riblet case, we selected Lagrange interpolation and two points from each fluid side of the interfaces to do interpolation. Note that the first adjacent fluid nodes near the interfaces were not included for stability reasons.

To calculate wall shear stresses on a riblet surface, velocity gradients need to be calculated on these surfaces, which are not necessarily on mesh nodes. To obtain velocity gradients, we carried out interpolation the same way as described above on one dimensional basis. Then we used a simple finite difference based on the interface point and an interpolated point about the local mesh size away. Once velocity gradients are obtained, the wall shear stress τ w can be calculated by τ w = μ ∂ V ∥ ∂ n | w (5) where V ∥ is the velocity parallel to the wall surface, and n is the normal direction to the wall surface. It is not obvious when calculating wall shear stresses on riblet surfaces than flat surfaces. For riblet surfaces, the normal direction is 0 n 2 n 3 T , the partial derivative ∂ V ∥ / ∂ n can be computed by ∂ V ∥ ∂ n = ∂ u ∂ y n 2 + ∂ u ∂ z n 3 (6)

Then, the drag can be expressed as a line integral (actually a surface integral, however we can simply take average in streamwise direction because of the homogeneity) D = ∫ 0 L μ ∂ V ∥ ∂ n d s = μ ∫ 0 L ∂ u ∂ y n 2 d s + ∫ 0 L ∂ u ∂ z n 3 d s (7) where L is the length of the curved surface in the spanwise direction. Transforming the curve integral into coordinate integral we get D = μ ∫ 0 L 3 ∂ u ∂ y d z + ∫ 0 L 2 ∂ u ∂ z d y (8) where L 2 and L 3 are the normal and spanwise extent of the riblet surfaces. Note before using Eq. 8 to integrate, riblet-wise average should be carried out and the symmetry of the riblet surface used first to make the integrand single-valued in the y direction. Finally, to be comparable to the flat plate case, the skin friction coefficient c f is defined as c f = D / L 3 1 2 ρ u b 2 = 2 D ρ u b 2 (9) where u b is the bulk velocity.

We considered the kind of scalloped riblets constructed by smoothly connecting two third-order polynomials, and selected a scalloped riblet with shape parameters A = 25 and B = 1.6 . The turbulent channel flows with and without riblet control were simulated with W + = 20 , and height-width ratio γ = 0.5 . The drag reduction rate was 5.77%, which is generally larger than corresponding triangular riblets with sharper tips. The flow fields were then investigated carefully. Mean flow fields, Reynolds stress, correlations of vorticity and Liutex components, pre-multiplied spectra of streamwise velocity and Liutex component, and instantaneous flow fields were presented. It was found that streamwise vortices just above the riblet tips, which have a length scale of 200–300 in wall units, play a significant role in the controlling the flow field. Further investigations, especially the causal relations between those streamwise vortices and the drag reduction, should be analyzed. Moreover, it is imperative to consider the influence of the Reynolds number in the context of scalloped riblets, particularly in high Reynolds number flows, where the application of riblets for drag reduction is targeted.