Numerical investigation of the drag reduction effect in turbulent channel flow by superhydrophobic grooved surfaces

Abstract

Superhydrophobic surfaces (SHSs) are considered to be a promising technology for achieving skin-friction drag reduction. Development of more efficient techniques for simulating the turbulent boundary layer on SHSs continues to be a subject of interest. In this study, numerical simulations were carried out to capture near-wall behaviours due to the effect of the SHS on wall-bounded flows. To achieve this, high- to intermediate-fidelity turbulence models including Reynolds-averaged Navier–Stokes, detached eddy simulation and large eddy simulation were utilized. With regard to slip conditions, the well-known Navier slip velocity method was used over the SHS. For validating the numerical solutions, the slip velocity and skin friction over the SHS were compared with the experimental output. Results showed that the velocity profile and Reynolds stresses on the SHS were comparable to the reported results. Then, the developed models were further extended to investigate the drag reduction effect of SHSs with rectangular grooves. The subsequent results showed that the combination of superhydrophobicity and rectangular grooves led to a better performance with a maximum drag reduction of 46.1%. This is due to the surface slip caused by the SHS and the secondary vortex effect created by the grooves. Our results revealed that Reynolds stresses of the slippery grooved surface were higher than those of the case in which a shear-free condition was employed for the grooved surface. More importantly, the numerical results indicate the previous assumption of the shear-free condition is inaccurate for the geometrically simplified grooved SHSs. Therefore, geometry modifications rather than an overly simplified shear-free boundary condition should be applied in computational fluid dynamics simulations for SHSs with grooves or other complex structures.

Impact Statement

Computational fluid dynamics for drag reduction applications has been a subject of interest in engineering and industries. Superhydrophobic surfaces have been previously investigated as a passive technique to reduce drag. The paper performs high- to intermediate-fidelity turbulence models to predict the flow field around superhydrophobic surfaces to support further understanding of their drag reduction efficiency. Our paper addresses the challenge of integrating superhydrophobicity and grooved surfaces in simulating turbulent flow using the well-known Navier slip boundary condition. Our simulations were implemented to observe the near-wall behaviours in flat and grooved surfaces and reveal further insight into the effects of the superhydrophobic surfaces on the drag. The adoption of these numerical models for superhydrophobic surface design by industry and engineers will enable more reliable estimations of the optimal drag reduction that minimizes the energy consumption in vehicles, contributing to the development of this technology for energy saving.

1. Introduction

Skin friction accounts for a great portion of the total hydrodynamic resistance and fuel consumption on vehicle surfaces. Nearly 60% of the vessel's fuel is expended to overcome the frictional drag on the wetted surface (Mäkiharju, Perlin, & Ceccio, 2012). Techniques to reduce skin friction can therefore produce substantial savings in fuel consumption and operating costs through improvements in marine vehicles’ speed and efficiency. For decades it has been well established that a stable air layer over a solid surface can reduce the skin-friction drag in liquid fluid flows (Ceccio, 2010). Superhydrophobic surfaces (SHSs) have been previously investigated as a passive technique to reduce drag by forming an air layer at the solid–liquid interface (Rothstein, 2010). Whereas the usual no-slip condition exists in the regions of solid–liquid contact, the liquid–air interfaces bridging the roughness of the texture act effectively as slippery boundaries on SHSs. Consequently, the flow partially slips over the texture, and the net shear stress on the wall is decreased because of the reduced contact between the solid substrate and the fluid flow.

A great deal of experimental research has been carried out to explore drag reduction using SHSs. Representative experimental works include different flow geometries such as towing tanks (Aljallis et al., 2013; Taghvaei, Moosavi, Nouri-Borujerdi, Daeian, & Vafaeinejad, 2017), water channels (Gose et al., 2018; Ling et al., 2016), pipe flows (Pakzad, Liravi, Moosavi, Nouri-Borujerdi, & khani, 2020; Rad, Moosavi, Nouri-Boroujerdi, Najafkhani, & Najafpour, 2021) and Taylor–Couette apparatus (Rajappan et al., 2019; Srinivasan et al., 2015). For instance, Rowin et al. investigated the effect of slip boundary on the near-wall statistics of a fully developed turbulent channel flow over an SHS. Experimental results achieved a 37%–42% drag reduction in Reynolds numbers ranging from 6200 to 9400 (Rowin & Ghaemi, 2020). They also measured the slip velocity and near-wall Reynolds stresses using optical techniques.

In a parallel effort, numerical studies have been developed to explore the mechanism of turbulent drag reduction on SHSs. The majority of the work employed direct numerical simulation (DNS). Turbulent structures near the SHS have been regarded as a central challenge due to their dynamic and complex nature. A DNS simulation study showed that the number of spanwise vortical structures in the outer region of a SHS was reduced. The origin of the mean secondary flow in turbulent flows over the SHS revealed that the secondary flow was driven and maintained by the gradient of Reynolds stresses. As a result, SHS was considered as one of the key features in reducing drag in internal turbulent flows (Im & Lee, 2017).

Although DNS is regarded as the gold standard of numerical simulations (Hajisharifi, Marchioli, & Soldati, 2022), DNS remains computationally expensive for studying SHS in a turbulent flow. The DNS studies have been performed in just one or two Reynolds numbers and mostly concentrated on relatively low flow speed due to the extensive computational cost. Most of the computational effort in DNS is expended on resolving the small dissipative motions, whereas the energy and anisotropy are contained mostly in the larger scales of motion.

Hence, other turbulent methods have to be employed that can cover a wide range of Reynolds numbers with less computational expenses. Large eddy simulation (LES) is a method that does have a balanced computational cost and has been broadly used for numerical studies (Pope, 2000; Safari, Saffar-Avval, & Amani, 2018; Sagaut, Deck, & Terracol, 2013; Kadivar, Tormey, & McGranaghana, 2021). In LES, the dynamics of the larger-scale motions are resolved precisely, while the influence of the smaller scales is modelled. However, a few research studies have considered the LES method in the simulation of the turbulent boundary flow on SHSs. For example, the effects of micro-features on near-wall behaviours were investigated by using LES in the case of a fully slip condition on the SHS (Saadat-Bakhsh, Nouri, & Norouzi, 2017). Others have investigated the effects of SHS curvature at the liquid–gas interface by utilizing LES turbulence methods (Yao & Teo, 2020). However, the utilization of LES methods generally comes with certain complications, including the choice of filtering, closure modelling and near-wall treatment (Goc, Lehmkuhl, Park, Bose, & Moin, 2021).

The Reynolds-averaged Navier–Stokes (RANS) method is another approach to account for turbulent phenomena (Hanjalić, 2004; Wilcox, 1998) with much lower computational cost with respect to DNS and LES. It is acknowledged that, unlike Scale-resolving simulation (SRS) models, the use of a RANS solver cannot resolve the spectrum of turbulent length scales. However, the well-established k–ω model has been utilized widely in wall-bounded flows. The k–ω model is a two-equation eddy-viscosity model which represents turbulent flows by solving for turbulent kinetic energy and dissipation rate (Menter, Kuntz, & Langtry, 2003). For instance, the RANS results of previous research have shown a good agreement with DNS results in a microchannel with SHSs (Jeffs, Maynes, & Web, 2010). Moreover, the RANS method has also been used to explore the performance improvement of the tidal turbine by taking the advantage of slip effects of the SHS (Sun & Huang, 2020). However, RANS turbulence models are averaged in time while the transient turbulent structures can only be determined with the SRS methods. Therefore, the RANS models cannot accurately reproduce an entire flow field and it is not capable of resolving the important unsteady flow structures.

There are a limited number of studies that employed two-phase flow approaches to model the liquid–air interfaces in turbulent flows. By implementing the volume of fluid method, the air–liquid interface was simulated on a SHS. The effects of interface and roughness heights were observed and it was shown that the slip tended to shift the velocity profiles towards the SHS whereas roughness pushed it away from the wall. They concluded that roughness increased the negative shear stress and momentum mixing, while the interface prevented the development of shear stress. (Alamé & Mahesh, 2019).

In the literature, a shear-free (fully slip) condition, which implies a full air layer above the boundary, or Navier's slip (partially slip) condition, which corresponds to partial coverage of the SHS with an air layer, have been applied to SHSs. These conditions are applied without considering the liquid–air interfaces or employing multiphase-flow simulations. However, experiments have shown that it is nearly impossible to keep air entrapped as a fully air layer in turbulent flows (Heo, Choi, & Lee, 2021). Therefore, the shear-free boundary condition may not be accurate due to the failure of the air layer. Hence, a portion of the boundary may be in direct contact with water, and consequently, Navier's slip conditions can be established.

In DNS simulations, a simple smooth SHS can be modelled as a flat wall consisting of solid–water and air–water interfaces with streamwise and/or spanwise slips (Min & Kim, 2004). Furthermore, the effects of slip boundary condition in both the streamwise and spanwise directions on turbulent channel flow have been investigated parametrically. Results showed that a high spanwise slip length disrupted streamlines adjacent to the wall, while a high streamwise slip tended to straighter and more regular streamlines (Busse & Sandham, 2012). Textured SHSs with patterns can be modelled through boundary conditions by alternating regions that are either no slip or shear-free on the boundary. In the research by Park et al., the SHS was designed as an array of micro-grooved surfaces placed parallel to the flow direction (Park, Park, & Kim, 2013). The effect of pitch and air fraction (shear-free region) of the micro-grooved surfaces on turbulence structure alteration as well as drag reduction was studied (Park, Park, & Kim, 2013). It should be noted that geometrical complexity such as a textured surface or grooving may add complications in terms of accuracy and stability of the numerical simulation (Falcucci et al., 2021)

To the best of our knowledge, the previous studies did not aim to evaluate the suitability of various turbulence models as well as appropriate slip boundary condition over grooved SHSs. Accordingly, in this study, we assessed low- to high-fidelity turbulence methods including RANS, detached eddy simulation (DES) and LES, for the wall-bounded flow over a SHS. The validation study for the flat SHS proved the accuracy of Navier's slip condition in accounting for the flow characteristics such as shear stress and slip velocity. Moreover, in the case of a grooved SHS, a comparison between two common slip boundary conditions, i.e. shear free and Navier slip, has been made as another objective of this study.

The methodology and details of numerical simulations are described in § 2. Thereafter, in § 3, the validation results of the flat SHS are discussed by comparing the velocity, shear stress and Reynolds stresses. § 4 presents the flow structures on grooved SHSs and their drag reduction enhancement. The summary and concluding remarks are presented in § 5.

2. Computational methodology

2.1 Problem description and mesh generation

We set up our simulation for a three-dimensional rectangular channel flow that has been experimentally measured by Rowin and Ghaemi (2020). The geometrical details as well as channel dimensions are shown in figure 1(a). The channel length is L = 1.2 m, with a rectangular cross-section of $40\;\textrm{mm} \times 6\;\textrm{mm}$. To investigate the SHS effect in the fully developed flow regime, the channel surface is made superhydrophobic at 720 mm downstream of the channel entrance, with a length and width of L = 240 mm and W = 40 mm, respectively.

Figure 1. (a) Schematic view of the numerical set-up and main dimensions (Rowin & Ghaemi, 2020). (b) Generated mesh in the channel.

The utilized grids for each turbulence model are purely structured meshes, which were generated using the well-known feature-based blocking approach in the ANSYS ICEM package. As shown in figure 1(b), the domain is decomposed into three regions namely the inlet, SHS and outlet zones, in order to employ a sufficiently fine and uniform grid in the SHS zone. It is worth noting that the SHS zone not only includes the superhydrophobic boundary but also is extended by around 20 mm (twice the hydraulic diameter) on both ends to ensure fully developed turbulent structures adjacent to the SHS. Furthermore, the mesh size is gradually increased towards the inlet and outlet boundaries.

With regards to the RANS simulation, the whole domain including all three zones was modelled. Different mesh sizes were considered for each zone. As depicted in figure 1(b), the element sizes of the ‘Inlet zone’ and ‘Outlet zone’ decrease towards the ‘SHS zone’ with a rate of 0.9, while maintaining uniform element sizes at the SHS zone. This corresponded to approximately $6 \times {10^6}$ elements for the three zones.

However, implementation of the SRS requires a uniformly fine grid for the whole domain and thus would be computationally expensive. Hence, as for the LES and DES simulations, only the SHS zone was modelled. Accordingly, the inlet and outlet boundary conditions, i.e. velocity and pressure profiles, were extracted on those sections from the RANS simulation and were utilized for the SRSs. Finally, the DES and LES models consisted of $4 \times {10^6}$ and $3 \times {10^7}$ elements, respectively.

Table 1 shows the number of nodes in each direction for different scenarios. A fine grid should be considered near the SHS to resolve the near-wall phenomena. Accordingly, the dimensionless wall distance index, ${y^ + } = y{u_\tau }/\upsilon $, was kept less than unity for all cases to resolve the viscous sublayer phenomena. Moreover, a smooth grid without intense cell volume changes was also generated in the other walls.

Table 1. Domain size and mesh resolution for turbulent channel flow at Re = 9400.

^a As mentioned in figure 1, the RANS simulation was performed in the whole domain of the channel. For other cases, grid points were reported for the SHS zone only.

^b These two cases will be discussed in § 3.3 in detail.

2.2 Numerical methodology and governing equations

In the following, three approaches for modelling the turbulent flow over a SHS constituting RANS, DES and LES are briefly described.

To perform the numerical simulation of the SHS in channel flow, the filtered and averaged form of the Navier–Stokes equation is numerically solved using the well-known commercial computational fluid dynamics (CFD) code ANSYS Fluent. The unsteady RANS (URANS) formulation is adopted for the governing equations of the incompressible flow as follows:

(2.1)\begin{equation}\frac{{\partial \rho \overline {{U_i}} }}{{\partial t}} + \frac{{\partial \rho \overline {{U_i}}\, \overline {{U_j}} }}{{\partial {x_j}}} = {-} \frac{{\partial P}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {\mu \left( {\frac{{\partial \overline {{U_i}} }}{{\partial {x_j}}} + \frac{{\partial \overline {{U_j}} }}{{\partial {x_i}}}} \right) - \mathrm{\rho }\overline {{{u^{\prime}}_i}{{u^{\prime}}_j}} } \right],\end{equation}

where $\rho$ is fluid viscosity and P is pressure. Also, $\overline {{U_i}} $ and $\mathrm{\rho }\overline {{{u^{\prime}}_i}{{u^{\prime}}_j}} $ are the averaged velocity and Reynolds stress tensor, respectively. The latter can be written as

(2.2)\begin{equation}{-}\mathrm{\rho }\overline {{{u^{\prime}}_i}{{u^{\prime}}_j}} = 2{\mu _t}\overline {{S_{ij}}} - \frac{2}{3}k\rho \overline {{\delta _{ij}}} ,\end{equation}

where the mean strain rate tensor is

(2.3)\begin{equation}\overline {{S_{ij}}} = \frac{1}{2}\left( {\frac{{\partial \overline {{U_i}} }}{{\partial {x_j}}} + \frac{{\partial \overline {{U_j}} }}{{\partial {x_i}}}} \right).\end{equation}

The turbulent viscosity term, ${\mu _t}$, can be addressed with common one- or two-equation eddy viscosity models to establish turbulence closure.

As for the LES model, the filtered Navier–Stokes equation is

(2.4)\begin{equation}\frac{{\partial \rho \widetilde {{U_i}}}}{{\partial t}} + \frac{{\partial \rho \widetilde {{U_i}}\widetilde {{U_j}}}}{{\partial {x_j}}} = {-} \frac{{\partial P}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_j}}}\left[ {\mu \left( {\frac{{\partial \widetilde {{U_i}}}}{{\partial {x_j}}} + \frac{{\partial \widetilde {{U_j}}}}{{\partial {x_i}}}} \right)} \right] - \frac{{\partial {\tau _{ij}}}}{{\partial {x_j}}},\end{equation}

where the last term, ${\tau _{ij}}$, is the sub-grid-scale (SGS) stress and is defined as follows:

(2.5)\begin{equation}{\tau _{ij}} = \rho \widetilde {{u_i}{u_j}} - \rho \widetilde {{u_i}}\widetilde {{u_j}}.\end{equation}

There exist a couple of closure models to address SGS in terms of the local resolved flow. The selection of the appropriate turbulence closure model for RANS and LES is of great importance. The k–ω shear stress transport (SST) was used which yields the most accurate results in terms of RANS models for wall-bounded problems (Menter, Kuntz, & Langtry, 2003). In addition, all the spatial and temporal terms in (2.1) were discretized by the second-order method, and the semi-implicit pressure-linked equations algorithm was employed to solve the set of equations. Furthermore, regarding the LES equations, the wall-adapting local eddy-viscosity model was implemented to capture the SGS structures as it is the most suitable modelling technique for the LES computation of channel flows. In addition, as for the LES simulations, the convection term in (2.4) was discretized by the central differencing method.

In the DES approach, which is regarded as a hybrid turbulence model, the URANS and LES methods were implemented in the near-wall and detached regions, respectively (Krastev, Ilio, Falcucci, & Bella, 2018; Krastev, Silvestri, & Falcucci, 2017). The LES region is usually associated with the core turbulent area where large turbulence scales play a dominant role. In the near-wall region, a corresponding RANS model was considered, by taking advantage of the k–ω SST model. The DES requires a lower computational cost with respect to LES and it provides more accurate results than RANS, therefore, it is currently being more widely applied for numerical studies of the turbulent boundary layer (Spalart, 2021).

2.3 Numerical procedure

2.3.1 Boundary conditions

Setting appropriate boundary conditions is essential in obtaining a real physical model, especially in the case of SRSs. Regarding the RANS simulation, the inlet velocity and pressure outlet boundary conditions were considered for the start and end of the domain, respectively. The constant velocity corresponding to Reynolds numbers ranging from 6200 to 9400 was prescribed, while zero gauge pressure was considered for the outlet boundary. As mentioned in § 2.1, RANS simulations were performed first. Afterwards, the velocity, pressure and turbulence values at the start and end sections of the SHS zone were extracted. Then, these data were used as the boundary conditions of the following SRSs.

The Navier slip and no-slip conditions were applied to SHS and other wall boundaries, respectively. The Navier slip model, (2.6), is a conventional boundary condition for slip velocity at the wall, ${U_S}$, in numerical simulation of SHSs and is defined as follows:

(2.6)\begin{equation}{U_S} = {\left. {{b_x}\frac{{{\tau_w}}}{\mu }} \right|_{y = 0}},\end{equation}

where ${\tau _w}$ is the streamwise shear stress over SHS and ${b_x}$ is the slip length, which is a constant value that is evaluated from experimental data (Rowin & Ghaemi, 2020). Furthermore, this model depends on the local gradient of the velocity on the wall. Subsequently, we wrote a user-defined function to apply Navier's slip condition on the SHS boundary.

2.3.2 Iterative procedures with the Navier slip condition

The implementation procedure for the slip velocity is shown in figure 2. First, the no-slip condition was applied at the start of the simulation. By solving the Navier–Stokes equations once, the wall shear stress can be calculated afterward. Then, the slip velocity was calculated by applying the slip length relationship. Finally, the modified boundary condition was inserted into the second step to solve the Navier–Stokes equation. An under-relaxation method of 0.1 was utilized for the value of slip velocity as well to guarantee the convergence and stability of this procedure.

Figure 2. Flow chart of the simulation procedures for implementing slip velocity.

Simulations were conducted for five different inlet velocities for the RANS simulations, including 1.11, 1.25, 1.35 and 1.46 m s⁻¹. Afterward, the simulation was continued with the DES and LES methods separately to obtain the transient flow simulation. A Courant number of 0.3 was considered for all simulations to determine the time step value. Simulations continued until the stationary condition was achieved (typically eight flow cycles). The convergence criterion of ${10^{ - 6}}$ was considered for steady simulations, i.e. RANS model. Whereas for the DES/LES simulations, the solution was considered to be converged in each time step if the residuals of flow equations reach ${10^{ - 3}}$.

3. Results and discussion

In this section, the results for the flat and grooved SHSs are discussed. First, mesh verification is evaluated to ensure the appropriate element size. Then, the reliability of the utilized Navier slip condition will be examined by validating numerical results to the experimental data. Finally, the slip condition will be applied to the grooved SHS and the underlying physics of the turbulent flow discussed.

3.1 Mesh independency

The mesh independence analysis was first investigated to verify the performance and choice of numerical parameters. For the RANS simulations, to investigate the sensitivity of the numerical solution to the mesh size, three different meshes consisting of 4.5, 6 and 6.5 million cells were generated, and the variation of the wall shear stress and slip velocity was investigated. It was concluded that the difference in the average wall shear stress and slip velocity between the two latter meshes was no more than two per cent. We then chose 6 million cells to perform further simulations.

Moreover, the LES index of quality (LESIQ) proposed by Celik et al. was utilized to assess the LES mesh which is defined as $\textrm{LESI}{\textrm{Q}_v} = 1/(1 + 0.05{[(\mu + {\mu _{sgs}})/\mu ]^{0.53}})$, in which $\mu $ and ${\mu _{sgs}}$ are fluid and SGS turbulent viscosities, respectively (Celik, Cehreli, & Yavuz, 2005). The purpose of this parameter is to evaluate the proportion of total kinetic energy resolved by LES on each element. Accordingly, an LES computation can be judged to be well resolved when 80% of the turbulent kinetic energy is resolved, i.e. $\textrm{LESI}{\textrm{Q}_v} > 0.8$ (Gousseau, Blocken, & Heijst, 2013). As shown in figure 3, the value of $\textrm{LESI}{\textrm{Q}_v}$ ranges between 0.897 and 0.952, which is acceptable for the LES simulation (Celik, Cehreli, & Yavuz, 2005). Moreover, the value of dimensionless wall distance, ${y^ + }$, was calculated across the wall boundaries and is presented in table 1.

Figure 3. The LES index of quality contour for simulations in the middle plane along the channel flow, where SHS is located at the bottom.

3.2 Flat SHSs

Figure 4 shows the velocity contour along the channel, where the SHS is located at the bottom and the flow is in the positive x-direction. Different values of the slip velocity confirm that it is dependent on the local gradients of the flow field. Slight increases in the non-uniformity of the streamwise velocity fluctuation and turbulent kinetic energy were detected along the SHS surface. As illustrated in figure 4, Navier's slip boundary condition leads to inhomogeneous velocity fields on SHSs, enhancing turbulent structures. Enhanced vortical structures near the SHS can be observed from the comparison of figure 6(c,d) as well.

Figure 4. The LES result of the instantaneous velocity contour in channel flow where the SHS is located at the bottom. Four sections across the channel were created to visualize the velocity contour. The inhomogeneous velocity field on the SHS is due to the Navier slip boundary condition (Re = 9400)

3.2.1 Shear stress and slip velocity

The numerical results were compared with an experimental study by previous research (Rowin & Ghaemi, 2020) in which a two-dimensional particle image velocimetry (PIV) was used to measure the slip velocity and slip length, and near-wall Reynolds stresses were measured by using the three-dimensional particle track velocimetry method. Prior to demonstrating the near-wall behaviour of the SHS, the numerical results were validated using experimental data in terms of slip velocity, and smooth and superhydrophobic wall shear stresses. As shown in figure 5, in the presence of a slip velocity, the shear stress of the SHS was reduced compared with a no-slip surface. Overall, the DES and LES results correlated well with the experimental data. Similarly, the RANS results showed an acceptable accuracy, and therefore, all the simulation results were found to be within an acceptable error range. However, as mentioned by previous authors, RANS generally comes with primary shortcomings in terms of turbulence modelling. For instance, secondary flows induced by anisotropic turbulence cannot be precisely captured with this model (Cauwenberge, Schietekat, Floré, Geem, & Marin, 2015). Moreover, the RANS model was unable to accurately account for the fluctuations of the flow field (Salim, Ong, & Cheah, 2011) which led to underpredict turbulent mixing (Verma et al., 2019) and, as a consequence, higher velocities and shear rates near the wall were derived. Hence, as expected, the RANS model demonstrated higher amounts of error with respect to the other turbulence models.

Figure 5. Comparison of numerical results with experimental data for wall shear stress of (a) a no-slip (smooth) surface and (b) an SHS. (c) Slip velocity for different inlet flow velocities.

3.2.2 Mean velocity profiles

Figure 6(a) demonstrates the mean velocity magnitude profiles of the LES models over the smooth surface and SHS for a Reynolds number of 9400. The mass flow rate was kept constant through the simulations, therefore the velocity across the channel slightly changed due to the slip velocity of the SHS. In figure 6(a), the alteration of the mean velocities by the SHS was clear, in contrast to the results of the smooth (no-slip) surface. The near-wall velocities over SHS shifted upward due to streamwise slip, which eventually affected the flow velocity pattern of the entire channel.

Figure 6. (a) Mean velocity magnitude profile at Re = 9400 for smooth (no-slip) and SHSs. An SHS is located at y = 0. (b) A comparison between vorticity magnitudes of ‘Smooth’ and ‘Superhydrophobic’ surface cases in the vicinity of the SHS zone. Contour and vector of velocity at the plane x = 0.9 m for (c) the smooth (no-slip) surface. (d) An SHS.

To be more clear, the vorticity magnitude near the wall is plotted in figure 6(b). The SHS caused more turbulent vorticity and fluctuation near the surface and resulted in enhancing the inner layer mixing in the vicinity of the SHS. Hence, the vorticity magnitude of the SHS is higher throughout the normal direction with respect to the no-slip surface. Figure 6(c,d) shows the contour and vector of velocity for the smooth and SHS at Reynolds number of 9400. The velocity magnitude incorporated a higher value adjacent to the SHS than that of the smooth surface.

Semi-logarithmic profiles of ${u^ + }$ over the SHS for Reynolds numbers 6200 and 9400 are presented in figure 7. The DES and LES results are compared with the PIV measurements as well as the smooth, i.e. no-slip, surface profile (in logarithmic law). There was an upward shift of ${u^ + }$ for the SHS when compared with the logarithmic law, due to the slip velocity. The normalized mean streamwise velocity is, thus, defined as ${u^ + } = u/{u_\tau }$, where the frictional velocity $({u_\tau })$ is calculated from $\sqrt {{\tau _w}/\rho } $. Here, ${\tau _w}$ is calculated by averaging the shear stress value over the SHS. As shown in figure 7, our simulation code could successfully capture the velocity profile near the SHS. The velocity profile by DES can be regarded to be in the uncertainty range of the experimental results.

Figure 7. Semi-logarithmic profile of the mean streamwise velocity over the SHS at (a) Re = 6200, (b) Re = 9400.

3.2.3 High-order flow quantities

The high-order flow quantities such as statistical Reynolds stresses and the velocity fluctuations have been explored to further evaluate the influence of slip boundary conditions of SHS on the channel flow dynamics. The non-zero components of the Reynolds stress tensor, $\langle {u_i}{u_j}\rangle $, adjacent to the SHS in slip boundary cases were calculated for the DES and LES simulations. The profiles were normalized by the lowest value of the inlet velocity throughout all the cases, i.e. 0.95 m s⁻¹, and compared with experimental data (Rowin & Ghaemi, 2020). Subsequently, the LES results were in better correlation with experimental data in comparison with DES. Figure 8 compares the results of these two numerical approaches for a Reynolds number of 9400.

Figure 8. A comparison between LES and DES in terms of (a) streamwise, (b) wall-normal, (c) spanwise and (d) shear Reynolds stresses over the SHS.

Figure 9 compares the LES results of the components of the Reynolds stress tensor over the SHS for three Reynolds numbers of 6200, 8000 and 9400, to those of the experimental data. As expected, the Reynolds stresses increased with the increase of Reynolds number. With regards to the streamwise Reynolds stress, i.e. $\langle {u^{\prime 2}}\rangle $, the maximum value shifted toward the wall with the increase of Reynolds number because of a thinner inner layer at higher velocities. Moreover, due to the small wall-normal and spanwise gradients of the velocity components the Reynolds stresses of these directions were one order of magnitude smaller than that of the streamwise Reynolds stress.

Figure 9. The LES results of (a) streamwise, (b) wall-normal, (c) spanwise and (d) shear Reynolds stresses over the SHS.

By increasing the Reynolds number, the near-wall velocities shifted upwards because of the streamwise slip velocity, as shown in figure 7. Similar behaviours were also observed on Reynolds stresses along the SHS, as illustrated in figure 9. Moreover, the streamwise Reynolds stresses in the outer layer were reduced continuously due to the decrease in the production of the streamwise turbulence. On the other hand, the wall-normal and spanwise Reynolds stresses and Reynolds shear stress increased due to the increased vortex transport term, which was correlated to the production of the near-wall vortices.

3.3 Superhydrophobic rectangular grooved surfaces

After we validated our simulation code on a flat SHS in comparison with the experimental results (Rowin & Ghaemi, 2020), we describe how we further extended our simulation code to the grooved surfaces with superhydrophobicity, which is expected to be more effective in drag reduction (Ermagan & Rafee, 2018) (Fuaad & Prakash, 2019). In conjunction with the SHS, the grooved surface achieves better drag reduction owing to the secondary vortex effect. Streamwise rectangular grooves were located across the SHS with the same geometry configuration for simulation. The width and depth of these grooves were 500 and 250 μm, respectively, corresponding to the width-to-depth ratio of 0.5, which has been reported to be the optimal value (Soleimani & Eckelsb, 2021).

We applied two boundary conditions for the grooved surface configuration, as shown in figure 10. For the first case, the grooved partitions were assumed to be filled with air and therefore the shear-free boundary condition was applied to them, and the rest were considered to have a slip velocity. For the second case, the slip velocity (Navier slip) was applied to the entire surface of the SHS. The simulations for the grooved surface were done at Reynolds numbers of 6200, 8000 and 9400. To achieve good precision in the grooved regions, fine grids were generated (figure 10). The hexahedron structured meshes were utilized with the aid of the blocking technique using the ICEM CFD package.

Figure 10. The boundary conditions for: (a) air-filled SHS (shear-free boundary condition), (b) the grooved SHS (Navier's slip boundary condition). All dimensions are in mm.

3.3.1 Shear stress and slip velocity

A comparison of simulation results for different surface configurations is presented in table 2. The drag reduction was calculated based on the shear stress defined as $DR = \; (\tau _{wall}^0 - {\tau _{wall}})/\tau _{wall}^0$, where $\tau _{wall}^0$ and ${\tau _{wall}}$ are shear stresses measured on the smooth and superhydrophobic walls, respectively. The air-filled case experienced a higher drag reduction than the other scenarios. However, the assumption of air-filled grooves is revealed to be exaggerated, whereas the partially slipped (grooved SHS) case can be regarded as a more realistic representation of the model which experienced a drag reduction of 46% (12% more than a flat SHS).

Table 2. Drag reduction comparison between different surface morphologies for Reynolds number of 9400.

Moreover, the effect of changing the Reynolds number is presented in figure 11. The air-filled case incorporated higher slip velocities with respect to other scenarios even for the grooved SHS. This is due to the effect of the riblets’ sidewalls in decreasing velocity. However, the grooved SHS case demonstrated lower shear stress corresponding to the change of the vortical structures in the boundary layer. There will be further discussion on this topic in § 3.3.3.

Figure 11. Reynolds number dependence for the grooved SHS case of (a) slip velocity and (b) wall shear stress.

3.3.2 Mean velocity profiles

The streamwise velocity profiles of different surface configurations in the semi-logarithmic scale for Re = 9400 are shown in figure 12. In the viscous and buffer layers of the boundary layer, the velocity magnitude shifted upward for both air-filled and grooved SHSs in comparison with the flat SHS. Additionally, due to the shear-free condition in the air-filled model, the boundary experienced higher values of the slip velocity with respect to other cases.

Figure 12. Semi-logarithmic profile of the mean streamwise velocity over different SHS configurations for a Reynolds number of 9400.

3.3.3 Velocity fluctuations

In the following section, the high-order flow quantities including Reynolds stresses as well as turbulent structures are discussed in order to evaluate the effect of the surface grooves on the flow dynamics.

Reynolds stresses: the non-zero components of the Reynolds stress tensor for different surface configurations were compared to assess their contribution to turbulent structures. The grooved SHS demonstrated higher values of Reynolds stress with respect to other scenarios, as depicted in figure 13. As mentioned in § 3.2, on the flat SHS, the slip condition reduces the viscous drag force and enhances turbulence. Moreover, the riblets’ shape controls the vortical structures in the turbulent boundary layer, leading to less momentum transfer and shear stress. This grooved surface morphology enhances turbulent structures near the surface that contribute to the reduction of drag. Therefore, this effect causes more drag reduction in the grooved case over the flat case. As for the air-filled (shear-free) case, the underestimated shear stress may be attributed to the geometrical simplification. Hence, neglecting the effects of the grooves may lead to notable errors. The grooved surfaces significantly increased the intensity of the streamwise and spanwise turbulent fluctuations. Moreover, the peak of Reynolds stresses decreased and became closer to the SHS surface as for the grooved case.

Figure 13. Reynolds stress comparison for the flat, air-filled and grooved SHSs in terms of (a) streamwise, (b) wall-normal, (c) spanwise and (d) shear Reynolds stresses over the superhydrophobic region.

Turbulent structures: as emphasized by previous work (Soleimani & Eckelsb, 2021), surfaces with riblets alter the vortex structures in channel flows. The momentum of turbulent structures in the near-wall region is affected due to the existence of vortices in the riblets, and consequently, the flow field is altered in the spanwise direction. Principally, the visualization of coherent vortical structures can provide valuable insight into the flow field. Hence, to provide a better observation of the relationship between turbulence and drag reduction over the SHS, the representation of the vortical structures based on the line integral convolution (LIC) concept is depicted in figure 14. The LIC technique produces streaking patterns that follow the velocity vectors tangentially. More details of the implementation of this technique can be found in previous research on this (Cabral & Leedom, 1993). Accordingly, velocity contours and streamlines for different cases are shown in figure 14. Both the grooved and air-filled models exhibit higher velocity magnitudes and vortices near the SHS boundary (at the bottom). The grooved SHS generated even more turbulent eddies near the SHS resulting in better flow interaction with respect to the flat SHS. The increased wall-normal and spanwise Reynolds stresses for the grooved SHS are attributed to the generation of near-wall vortices as well. Moreover, in the grooved SHS case, some tiny vortical structures were observed near the SHS, in particular, inside the grooves. This proves that more vortical structures exist in the grooved SHS, which cannot be truly captured by assuming an air-filled model.

Figure 14. The LIC illustration at the plane of x = 0.9 m for (a) flat SHS, (b) air-filled grooved surface and (c) grooved SHS.

The three-dimensionality of the flow dynamics was evaluated in terms of the Q-criterion, which is directly derived based on the second invariant of the velocity gradient tensor. The abundant multiscale vortices of the SHS were determined by the Q-criterion approach to present turbulent structures (Krastev, Amati, Succi, & Falcucci, 2018; Zhan, Li, Wai, & Hu, 2019). In this method, the vortices were defined as areas where the vorticity magnitude is greater than the magnitude of the rate of strain with the following expression:

(3.1) \[Q = {1/2}(\|\varOmega\|^2 -\|S\|^2), \]

where Q > 0 represents the existence of a vortex. In (3.1), Ω is the vorticity tensor and S is the rate of strain.

Figure 15 shows the iso-surfaces of three-dimensional vortical structures for the flat SHS, air-filled and grooved SHS based on a normalized Q-criterion of 0.2. The vortical structures are elongated in the streamwise direction with a reduction in size near the wall boundaries (Im & Lee, 2017). The density of vortical structures in the wall region is dependent on the hydrodynamics of the wall. As indicated in the top and side views of the flat SHS and air-filled cases, the density of the vortical structures increases for the latter case. As mentioned before, this corresponds to the reduction of shear stress and an increase in turbulence. Moreover, another turbulence initiative for the air-filled case is the alternative velocity boundary condition along the spanwise direction on the SHS boundary, which may contribute to the density of the vortical structures.

Figure 15. Contour of the instantaneous three-dimensional vortical structures for (a) flat SHS, (b) air-filled and (c) grooved SHS cases. The cross-section in the x–y direction shows the side view of vortical structure from the bottom of the channel, where the SHS is placed, to the top wall. The cross-section in the x–z direction shows the contour in the top view from the sidewall to the centreline.

By comparing the air-filled and grooved cases, the latter constitutes a higher density of vortical structures. The vortical structures in the air-filled case are mainly dependent on the shear-free boundary condition. However, the geometrical manipulation, as exists in the grooved SHS, generates more vortical structures than the simple shear-free boundary condition.

4. Conclusions

The SHS alters the symmetry, peak magnitude and location of Reynolds stresses, which implies the presence of a significant slip velocity and consequently a high amount of the shear stress reduction on the SHS. Previous research studies have shown the great potential of SHSs for drag reduction applications and have examined the near-wall behaviour in turbulent flows as well. The main purpose of this study was to shed light on the underlying mechanisms of turbulent structure alternations by considering the SHS. High- to intermediate-fidelity turbulent models including RANS, DES and LES were evaluated to explore the superhydrophobicity for wall-related turbulent statistics by utilizing Navier's slip boundary condition. The boundary condition was then developed and validated in terms of the shear stress, velocity and Reynolds stresses.

The focus of this research differs from previous investigations by addressing the modelling of the turbulent channel flow in three important aspects. First, the Navier slip boundary condition was considered for the walls of the SHS. The velocity magnitude on the SHS was then calculated in each time step resulting in a better understanding of near-wall behaviours. Second, simulations were performed in a broad range of Reynolds numbers and the effects of turbulent models and SHS configurations were assessed. Third, the present work compared the grooved SHS case with frequent modelling approaches in the literature including shear-free and slip boundary conditions.

From our simulations, we found that the RANS solver cannot capture flow structure near the SHS as accurately as the LES approach. In contrast, the DES and LES methods captured boundary layer phenomena with excellent agreement with the experimental data. In addition, the shear stress derived from LES was comparable to experiments. However, DES was not able to resolve Reynolds stresses precisely. Moreover, the slip condition causes the turbulent structures to move toward the SHS. The utilization of the LES model can truly predict some of the physics in turbulent wall-bounded flow on the SHS without the need for performing DNS.

Finally, the simulation was implemented to observe the near-wall behaviour in grooved surfaces to achieve further insight into the effects of the SHS on the drag. Flow patterns of the grooved SHS were compared with that of the flat air-filled assumption at a relatively high Reynolds number. The grooved SHS results showed a drag reduction of up to 48.9% over SHS. Furthermore, Reynolds stresses components were considerably higher for grooved SHS than that of the flat, leading to a more efficient reduction of the shear stress. The peak of Reynolds stresses increased and relocated closer to the SHS in the slip model with respect to the air-filled case. For future work, more comprehensive studies of other geometrical parameters are required to investigate the dependence of drag reduction on other factors such as different texture shapes and a range of pitch dimensions.