Riblet Drag Reduction Modeling and Simulation

: One of the most interesting passive drag reduction techniques is based on the use of riblets or streamwise grooved surfaces. Detailed ﬂow features inside the grooves can be numerically detected only by Direct Numerical Simulations (DNS), still unfeasible for high Reynolds numbers and complex ﬂows. Many papers report the DNS of ﬂows on microgrooved surfaces providing fundamental details on the drag reduction devices, but all are limited to plate or channel ﬂows far from engineering Reynolds numbers. The numerical simulation of riblets and other drag reduction devices at very high Reynolds numbers is difﬁcult to perform due to the riblet dimensions (microns in aeronautical applications). To overcome these difﬁculties, some models for riblet simulation have been developed in recent years, due to the data provided by DNS, experiments, and theoretical analyses. In all these models, the drag reduction is modeled rather than effectively captured; however, the analysis of some nonlocal effects on practical aeronautical conﬁgurations with riblets, requires their adoption. In this paper, the capabilities of these models in predicting riblets’ performance and some interesting features of the riblets’ effect on form drag and shock waves are shown. Two models are discussed and compared showing their respective advantages and limitations and providing possible enhancements. A comparison between the two models in terms of accuracy and convergence is discussed, and two new formulae are proposed to improve one of these models. Finally, a review of the results obtained by the two models is provided showing their capabilities in the analysis of the riblet effect on complex conﬁgurations.


Introduction
The requirements of pollution emission reduction is driving aviation and other industries towards a growing interest in drag reduction technologies. One of the most interesting passive drag reduction techniques is based on the use of riblets or streamwise grooved surfaces. These are very likely the only aircraft profile drag reduction that have the possibility of being applied in the next generation of aircraft [1]. Detailed flow features inside the grooves can be numerically detected only by Direct Numerical Simulations (DNS), still unfeasible for high Reynolds numbers and complex flows. There are several papers reporting the DNS of flows on microgrooved surfaces, providing fundamental details on the drag reduction devices, but all of these are limited to plate or channel flows far from engineering Reynolds numbers. Large Eddy Simulations (LES) with riblets installed were performed by Martin and Bhushan [2] on a flat plate and by Zhang et al. [3], who employed an implicit LES to predict riblet performance on an airfoil. However, the resolved LES still require great computational efforts at a high Reynolds number, and it is still unfeasible on aircraft configurations. The numerical simulation of riblets and other drag reduction devices at very high Reynolds numbers is difficult to perform due to the riblet dimensions (microns in aeronautical applications). To overcome these difficulties, some models for riblet simulation have been developed in recent years, due to the data provided by DNS [4,5], experiments [6][7][8][9][10], and theoretical analyses [11]. Aupoix et al. [12] modified the Spalart-Allmaras and k − ω turbulence models for riblet simulations; the model was adapted by Koepplin et al. [13] for turbomachinery applications and the misalignment of riblets. Mele and Tognaccini [14] modeled riblets as a singular roughness problem in the transitional regime [15], adopting a wall boundary condition for the turbulence model. Results in agreement with the experiments in the case of a flat plate, airfoil, and practical aeronautical configurations in a wide range of Reynolds numbers were presented in [16][17][18]. Jiahe et al. [19] adopted a similar method to compute the riblet effect on a missile surface. All these models are linked to turbulence modeling with the related limitations. Different methods have been recently introduced in [20,21]. Ran et al. [20] proposed a free-simulation approach based on volume penalization in the Navier-Stokes equations with the help of eddy-viscosity models for the turbulence effect on the mean quantities. Wang et al. [21] adopted an Artificial Neural Network (ANN)-based approach in a Lattice Boltzmann method to extract the micro flow characteristics in the near-wall region.
Following the physical mechanism described in Luchini et al. [11] and in Bechert et al. [7] where the drag reduction was characterized in terms of the difference between the two so-called protrusion heights in the longitudinal and cross-flow directions, Mele and Tognaccini [22] proposed a slip boundary condition for riblet simulations. The two protrusion heights, defined as the distance between the tips of riblets and the average of the velocity profiles, are equivalent to the slip lengths. A new law for protrusion height difference (i.e., slip length) as a function of riblet shapes valid also in the nonlinear region of drag reduction curve was derived. The model was validated in flat plate and airfoil flows showing a good agreement with the experiments. The capabilities of the model in predicting the riblet effect were confirmed in a more recent paper [17] employing a resolved LES.
The link between the slip flow and drag reduction mechanism has been widely discussed in the literature. Min and Kim [23] implemented a slip boundary condition in a DNS solver to compute the effect of a hydrophobic surface on skin-friction drag. DNS by Rastegari and Akhavan [24] showed that super-hydrophobic microgrooves and riblets share a common mechanism of drag reduction related to a slip flow. Longitudinal and transverse protrusion heights (i.e., slip lengths) in superhydrophobic surface flows to be adopted in DNS were computed by Alinovi and Bottaro [25]; the homogenization model was recently enhanced employing a high-order approach in [26]. The effective slip length in the case of a super-hydrophobic surface was calculated by Chang et al. [27] finding that the slip was anisotropic in the case studied. More recently Zhang et al. [28] proposed a slip model for riblets, performing LES on different riblet shapes, showing that the drag reduction by riblets was well reproduced by the slip model. Luchini [29] adopted the slip length concept to propose a linearized boundary condition in DNS solvers for the simulation of rough walls, allowing calculations without resolving the surface roughness in the grid.
In the present paper, the models introduced in [16,22] are discussed, analyzing their limitations and advantages. A comparison between the two models is performed for the first time providing useful details for their adoption in CFD codes. Two new formulae are also proposed for the slip boundary condition that introduce a substantial improvement to the original formulation. Finally, a review of the main results obtained in the numerical simulation of the riblets on complex aeronautical configurations are shown. The results show the capabilities of these models in predicting riblet performance and in analyzing some interesting features of riblet effects.

Governing Equations
The results of the numerical simulation of riblets presented in the next sections were obtained by a RANS (Reynolds Averaged Navier Stokes) and an LES solver. In the case of the RANS computations, a CFD (Computational Fluid Dynamics) code that solves the compressible steady or unsteady RANS equations was employed. The continuity, momentum, and energy equations are reported here: where H = E + p/ρ is the total enthalpy and q the thermal diffusive flux vector, τ is the dissipative stress tensor, and the bold symbol indicates a vector. The fluid dynamic variables are mean quantities obtained by a Favre averaging [30] of the Navier-Stokes equations. The thermodynamic model is the ideal gas model, the dissipative stress tensor is τ = (µ + µ T ) ∇V + (∇V ) T , where µ T is the turbulence viscosity computed by the turbulence model equations, and µ is the molecular viscosity µ computed by the Sutherland law. The k − ω SST [31] (Shear Stress Transport) turbulence model has been adopted in all calculations returning the best performance for the turbulence resolution in external flows. All results refer to the external compressible flow calculations; in this case, the only boundary condition at the boundaries far from the body is a "far-field" boundary condition based on the characteristic curves theory (see [32], for instance). The boundary conditions on the bodies will be discussed in the next sections. The flow field is always initialized with free-stream conditions. Both the experiments and theory highlighted that the physical mechanism responsible for the drag reduction by the riblets is local [7,11], i.e., it only depends on the local Reynolds number, and it can be summarized in a modification of the constant ∆U + in the log-law of the turbulent velocity profile: where the superscript + specifies the nondimensional quantities obtained by using wall variables, k = 0.41 is the Kármán constant, and B = 5 is the value of the constant adopted for smooth surfaces. In the next sections, two methods for riblet simulations are presented.

Turbulence Model Based Boundary Condition
∆U + can be introduced in different ways. Mele and Tognaccini [14,16] introduced the idea to model riblets as an ordered roughness in the transitional regime. They proposed a boundary condition for ω in the k − ω turbulence models written as a function of where A + g is the riblet nondimensional cross-section area. In fact, following Mayoral and Jiménez [33], l + g provides a better characterization of riblet performance than the nondimensional riblet spacing s + and height h + . The boundary condition for ω in the k − ω turbulence models can be written as [34]: where ρ is the density, u τ = τ w /ρ is the friction velocity (τ w is the wall shear stress), and µ is the dynamic viscosity. In the case of riblets, this law for S R has been proposed: The coefficients were obtained by numerical experiments matching the experimental data reported in [33]: C 1 = 2.5 × 10 8 ; C 2 = 10.5; C 3 = 1.0 × 10 −3 ; n = 3. C 2 is equal to the value of l + g corresponding to the maximum value of S R , while C 1 and C 3 are related to the maximum value of S R .

Slip Length Based Boundary Condition
The literature discussed in the introduction showed that drag reduction devices such as riblets, streamwise-traveling waves of spanwise velocity or a super-hydrophobic surface, and also roughness induce a slip flow that can be correctly modeled by a slip boundary condition. The slip boundary condition relates the components of the velocity tangent to the surface to the shear rate at the surface through the so-called slip length λ: u w = λ ∂u ∂y w . By increasing λ, the shift in the log-law ∆U + increases. The relation between λ and l + g derived by Mele and Tognaccini [22] is: where , and Re λ is the Reynolds number based on λ.
Two new equations are proposed. The first one is a slight modification of Equation (7) with different coefficients and the introduction of a new term for a better characterization of riblet performance in the drag increase regime: where The second equation links λ/s to l + g , where s is the riblet spacing: where C 1 = 8, C 2 = 90, and C 3 = 0.0023, while C 4 and λ 0 are the same as in Equations (7) and (8).
In the next section, the results obtained applying the modified equations will show the improvement in the characterization of the riblet performance.

Performance Comparison
There is an interesting symmetry between the two models that is clearer when considering their numerical implementation. In the case of the slip length model: where u 1 is the velocity of the first inner point near the wall, and d is its distance from the wall. The finite difference expression of Equation (6) is: showing a clear symmetry with the equivalent slip condition Equation (10), which indeed only differs for a scale factor (u w is a velocity, and ω w is a reciprocal of time). It is worth remarking that the two models, as they have been built, are alternative, i.e., each model must be applied alone. The first evidence when comparing the two models is that the ω-based model is isotropic, while the slip length model can take into account the riblet alignment with the flow direction.
The performance of the two models in terms of drag reduction prediction are quite similar. Flat plate computations were performed on a grid with 640 cells in the streamwise direction and 256 cells in the normalwise direction, and 512 cells were distributed along the flat plate. In Figures 1 and 2, the computed drag reduction applying Equations (6)- (9) is compared with the experiments and a DNS by García-Mayoral and Jiménez [33]. In particular, in Figure 2, the improvement obtained by adopting the new formulae is evident. Equation (9) can be easily adapted for the drag reduction prediction of a specific riblet shape. In Figure 3, the computed drag reduction in the case of blade riblets is compared with the experiments of Bechert et al. [7] showing a very good agreement. Even if both models reproduce the drag reduction curve well, the ω-based model is simpler to apply to the drag increase regime, while the slip length model should have some kind of adaptation to obtain a drag increase. On the other hand, the slip length model is not linked to a turbulence model and can be adopted in an LES or DNS solver. A limitation of the ω-based model is that the values of ω needed to achieve a drag reduction of 6-8% are very large, and this may lead to lack of accuracy. Indeed, a different solver found a strong dependency of the drag reduction on grid y + until y + ≈ 0.25 was reached. This is not surprising because the large values of ω at the wall imply strong gradients near the wall; thus, if the number of grid points inside the boundary layer is not sufficient, the gradient computation may be inaccurate. In Figure 4, a grid dependency study in the case of the flat plate flow is shown. The finest grid had 640 cells in the streamwise direction and 256 cells in the normalwise direction, 512 cells were distributed along the flat plate; the coarser grids were obtained by halving the number of grid cells. The coarser grid had y + ≈ 1; it is evident that the slip length model was almost grid independent, while the ω-based model converged at the finest mesh size. However, it has been found that this effect is less evident in 3D simulation cases. It is also clear that this aspect limits the amount of drag reduction that can be predicted by the ω-based model, differently from the slip length model that can theoretically predict any percentage of drag reduction.  It is well known that the drag reduction due to the riblets depends on the Reynolds number; in [35], a theoretical analysis showed that there was an evident decay of riblet performance with the Reynolds number; however, the reduction in performance decreased with the Reynolds number. For instance, the decrease in drag reduction from the Reynolds number 3 × 10 6 to 10 7 was less than 1%. Equation (8) has an explicit dependence on the Reynolds number, while in Equation (9), the dependence on the Reynolds number is implicit through riblet spacing s. The results of the flat plate computations reported in Figure 5 show that Equation (9) better reproduced the theoretical Reynolds number dependence.

Effect of Riblets on Form Drag and Shock Wave
It has been shown by experiments and theoretical analysis that the physical mechanism responsible for the drag reduction by riblets is local, and that the riblets act on friction drag. However, some experiments investigating the flow over a flat plate under an adverse pressure gradient showed an increased effectiveness of riblets [36], a result also confirmed by Large Eddy Simulations [37]. The increasing efficiency of riblets in airfoil flows has been measured in other experiments, where an additional drag reduction while increasing the angle of attack was found, implying an effect of the pressure distribution on riblet performance. A variation in the displacement thickness in the DNS computation of two different drag reduction techniques was reported in Stroh et al. [38]. The slip length concept was also applied in a DNS by Banchetti et al. [39] who analyzed the drag reduction over a curved wall, applied by streamwise-traveling waves of spanwise velocity. In this case, the authors also found a pressure drag reduction that improved the drag reduction performance. Very recently, Mollicone et al. [40] showed how superhydrophobic surfaces reduced form drag in bluff bodies. In [22], this effect was analyzed in detail, applying the slip length model previously described in a RANS solver. In [35], a more reliable analysis was carried out employing a resolved LES, which confirmed the results. The analyses provided a possible explanation for the increased performance of riblets in pressure gradient flows. It has been shown that riblets induce small but significant modifications of the pressure distribution, which tends towards the inviscid; this effect leads to a form drag reduction in addition to the well known friction drag reduction. Indeed, the log-law is not influenced by the pressure gradient; however, the boundary layer developed with the modified ∆U + has a secondary influence on the outer inviscid flow and on the pressure distribution in particular, as is well known by Prandtl's boundary layer theory. In Figure 6, the LES computation results show the effect of riblets on the pressure coefficient distribution along three different airfoils; it is evident that the pressure recovery at the trailing edge and the expansion peak are influenced by the riblets. A quantitative analysis with the help of the classical matched asymptotic expansion technique shows that the riblets reduce the displacement thickness of the boundary layer. The reduced thickness of the equivalent body leads to the form drag reduction. In Figure 7, the effect of the riblets on the displacement thickness is shown. The significant influence of riblets on the location and strength of the shock waves was first observed in [16], where the ω-based model was applied to a transonic wing body configuration in the NASA Common Research Model (CRM). In Figure 8, the effect of the riblets on the shock wave of a wing section is shown. The effect of the riblets on the shock wave position and strength is evident. A deeper analysis of a typical transonic test case, the RAE 2822 case 9, showed that the riblets may induce a separation. In Figure 9, the pressure coefficient and skin friction coefficient show how the shock wave was stronger and moved downstream with the riblets installed also causing a shock-induced separation. A very recent DNS computation [41] confirmed the effect of the reduced skin friction on the shock wave strength and position due to a drag reduction device. In the next sections, the effect of riblets on the form drag and shock wave of a full aircraft configuration will be shown.

Modeling the Riblet Effect on Aircraft Configuration
The boundary condition (6) was applied to full aircraft configurations in [16][17][18], showing how this model was able to analyze in detail riblet performance on practical aeronautical configurations. In particular, a transonic wing body configuration (NASA CRM) subject of the fifth AIAA CFD drag prediction workshop (DPW5) and two different wing body configurations designed to have a wide Natural Laminar Flow (NLF) along the wing were analyzed. In the case of the NASA CRM transonic configuration, the effect of the partial riblet installation was analyzed. The medium grid of DPW5 was used providing y + < 1 over the aircraft surface. In Figure 10, the lift drag polar curves with riblets installed on the wing, on the fuselage, and on both the wing and fuselage are shown. The analysis helps to evaluate the cost-benefit of riblet installation on single parts of the aircraft.
In the case of the NLF wing body, the effect of the riblets together with the NLF technology was considered. In Figure 11, the wing body configuration is shown with the pressure coefficient distribution on the body and the module of x-vorticity in the wake.
A structured mesh with 40 millions cells at the finest level was used. A grid convergence analysis while reducing the mesh size was first performed; in Figure 12, the computed polar curves for the three grid levels (h = 1 is the finest) shows that the medium grid (h = 2), with 20 millions cells, returned good accuracy, and it was used for the following calculations.   In Figure 13, the lift drag polar curves in the case of an NLF wing body with and without riblets are shown. The NLF conditions were obtained by imposing the transition line on the wing using data provided by the designers; the production terms of the turbulence model were set to zero in the laminar zone [42]. It is interesting to note that also in the case of NLF, riblet installation in the turbulent zones returned a significant additional drag reduction; indeed, in terms of drag reduction in cruise condition (C L = 0.5), 29 drag counts were computed with NLF only against 40 drag counts with NLF and riblets. It must be remarked that the riblets were installed both on the wing and fuselage with constant lg + . With this choice, the model could evaluate, as a result of the numerical simulation, the optimal riblet height distribution over the aircraft surface ( Figure 14). A drag breakdown by far-field methods [43,44] allowed the analysis of the effects of the riblets on the form drag and on the shock wave. In Figures 15 and 16, the viscous drag breakdown showed that the reduction in form drag contributed to the total drag reduction, and that the decrease in the total drag reduction with the angle of attack was certainly due to the increase in induced drag on which the riblets did not act but also to the decrease in form drag reduction.
The strong increase in the form drag at a 4 degree of angle of attack was due also to the presence of a shock wave not present at the previous angles of attack. The analysis of the shock wave with and without riblets provides interesting details. In Figure 17, the shock wave zone is shown (the supersonic zone is in red). It can be noted that the supersonic zone in the case of NLF control and without NLF control with riblets was slightly wider with respect to the case without NLF and riblets (no flow control). The computed wave drag in the case of NLF control (21 drag counts) was substantially the same as that computed in the case of the fully turbulent flow with riblets (22 drag counts). This result shows that the effect of riblets on the shock wave was similar to the effect due to the NLF control and must be considered in the aircraft design phase.

Conclusions
Two models for riblet simulation were analyzed, and two new formulae for the slip length model were proposed. The new formulae returned a better characterization of riblet performance; in particular, the equation that linked the slip length to the riblet spacing ratio to l + g better accounted for the Reynolds number dependence of the riblets' performance. The comparison between the ω-based and slip length models showed a greater adaptability of the slip length model. The comparison with the experiments confirmed that both models reproduce the riblets' performance in term of friction drag reduction in the whole range of l + g well. The application of these models to complex flows allows a deep analysis of additional riblet effects, which had not been considered in previous experiments and theoretical analysis, such as the effect on the form drag and shock wave. It was shown that the application of riblets on the fuselage could be useful, and that riblets provide additional drag reduction when applied to the NLF configuration. The present analysis shows that, even if the drag reduction is modeled rather than effectively captured, these methods are very useful in providing insights into some unexpected effects of riblets.