Can large-scale oblique undulations on a solid wall reduce the turbulent drag?

In the context of greener and more cost-effective aviation, industrial and academic researchers have proposed and studied a wide range of drag-reducing control methods mainly over the past three decades. Unfortunately, hardly any offers realistic prospects of being implemented in practice. This applies, in particular, to active control schemes, despite some successful implementations of active devices in laboratory tests.

Among passive concepts, the use of riblets was originally inspired by the narrow grooves observed on sharks’ placoid scales. Although the effectiveness of the dermal denticles of the shark has been questioned by Boomsma and Sotiropoulos,⁸ the use of optimally chosen longitudinal grooves, aligned with the main flow direction, has been shown to be capable of reducing the turbulent skin-friction drag by levels of order 5%–10%.^5,6,11,14 However, a practical, cost-effective implementation has yet to be achieved, mostly hindered by the small optimal spacing required (about 15 $μ$m in cruise-speed conditions) and stringent tolerances on the sharpness of the crests. More complex variants, such as sinusoidal riblets, were studied by Bannier,⁴ Kramer et al.,²¹ and Peet et al.,²⁷ but despite attempts to optimise the geometry, Bannier⁴ showed that conventional (straight) riblets appear to be as effective within the uncertainty margins.

On the active-control front, based upon the work of Jung et al.²⁰ on the drag-reducing properties of transverse wall oscillations, it has been established, computationally, that imparting streamwise-modulated, spanwise in-plane wall motions of the form of a travelling wave $ww(x,t)=Asin(2πx/λx−ωt)$⁠, giving rise to a “Generalised Stokes Layer” (GSL), results in gross friction-drag-reduction levels of up to about 45%^29,30 at low Reynolds numbers, the effectiveness being observed to reduce at higher Reynolds numbers.^16,19,39 An experimental confirmation of the numerical results for the travelling wave was undertaken by Auteri et al.³ in pipe flows, for which drag-reduction levels of up to 33% were achieved, while Bird et al.⁷ designed a Kagome-lattice-based actuator for boundary layers, achieving a drag-reduction level of around 20%. However, it is very challenging to implement the latter method in a physical laboratory, let alone in practice.

A particular case of the GSL is the Spatial Stokes Layer (SSL), consisting of a standing wave $ww=ASSLsin(2πx/λx)$⁠, as shown in Fig. 1. This method was studied by Viotti et al.⁴⁰ by means of direct numerical simulation (DNS) for various forcing amplitudes $ASSL$ and wavelengths $λx$⁠. The maximum net-energy savings of 23%, achieved as a result of Viotti et al.’s exploration at $Reτ=O(200)$⁠, was for a forcing wavelength of $λx+=O(1000)$ that is about two orders of magnitude larger than the optimal dimension of riblets. Thus, while the resulting control method is still active, it is steady, and based on geometric dimensions that, for an entirely passive device, would be compatible with a practical implementation.

In an effort to address the need for practical control methods, the present research examines a passive means of emulating spanwise in-plane wall motions, as proposed in Ref. 9. The geometry considered is a solid wavy wall, with the troughs and crests skewed with respect to the main flow direction. The presence of the skewed wavy wall creates a spanwise pressure gradient that forces the flow in the spanwise direction, thus generating an alternating spanwise motion, as shown in Fig. 2. In contrast to the SSL, where the wall is actuated, the velocity has to vanish at the solid wall so that the spanwise forcing can obviously not be faithfully emulated. However, the premise is that the wavy geometry will generate a spanwise shear strain, somewhat away from the wall, that will weaken turbulence in a similar manner to that affected by the SSL. Such a passive device would benefit from the favourable actuation characteristics of the SSL (large wavelength), resulting in a practical solution, from a manufacturing and maintenance standpoint.

The present study will focus on selected direct numerical simulations of turbulent wavy-channel flows with the aim of examining the degree of Stokes-layer emulation achieved, and the degree to which the drag is reduced relative to the plane channel. As part of this study, some major similarities and differences between the flow arising from in-plane wall motions and that past a wavy wall, as shown in Fig. 2, will be brought to light, including the impact on the near-wall turbulence.

An obvious problem is posed by the absence of any guidance on which combination of geometric parameters offers the promise of maximum drag reduction. An exploration of the three-dimensional parameter space (wave height, wavelength, and flow angle) by a “carpet-bombing” strategy, or classical formal optimisation strategy, is not tenable on cost grounds, especially because of the tight resolution requirements needed for an accurate prediction of the drag increase/decrease margin. For this reason, a preliminary low-order study was undertaken by Chernyshenko⁹ to narrow down the exploration range within which the drag reduction might be maximised. More details on the assumptions and results of this analysis are given in  Appendix A, leading to the estimate for the streamwise-projected wavelength of the wave $λx+≈1500$ and the flow angle $θ≈52°$⁠. However, the analysis does not provide an estimate for the height of the wave. Rather, the condition given is that the forcing amplitude of the emulated SSL should be $ASSL+=2$⁠. In the simulation presented of this particular configuration, the wave height was adjusted so as to approximately satisfy the above-mentioned condition.

The present simulation programme was initiated with the flow configuration arising from the preliminary analysis. Other configurations, a selection of which will be presented below, were later simulated, and the exploration across a reasonable parameter range was mainly undertaken by trial an error, guided by observed trends of emerging results.

Direct numerical simulations were performed using an in-house code, which features collocated-variable storage, second-order spatial approximations implemented within a finite-volume, body-fitted mesh. The equations are explicitly integrated in time by a third-order gear-like scheme, described by Fishpool and Leschziner.¹³ In the fractional-step procedure, the non-solenoidal intermediate velocity field is projected onto the solenoidal space by solving a pressure–Poisson equation. The latter is solved by successive line over-relaxation (Hirsch,¹⁸ p. 510), the convergence of which is accelerated using a multigrid algorithm by Lien and Leschziner.²² The multigrid iterations within any time step are terminated when a convergence criterion, based on the RMS of the mass residuals, made non-dimensional using the fluid density, bulk velocity, and channel half-height, decreased to 10⁻⁹ in all the simulations presented. Stable velocity–pressure coupling is ensured by the use of the Rhie-and-Chow interpolation,³¹ preventing odd-even oscillations.

Simulations were initiated with a laminar Poiseuille profile to which random perturbations were added. The viscosity was then temporarily decreased to allow the growth of instabilities and to prevent relaminarisation. Alternatively, when available, the simulation was restarted from a set of instantaneous snapshots taken from a previous simulation. Throughout the simulation, statistics are obtained over a fixed time interval and saved periodically. This approach allowed a judgement to be made when the transient period ended, and the statistics were discarded, prior to the assembly of definitive statistical data.

The code has been thoroughly verified and validated. A verification of the spatial accuracy of the code was undertaken via the method of manufactured solutions,^33–36,38 which indicated a second-order spatial accuracy for the velocity and pressure fields. The manufactured solution was implemented in a channel with lower and upper walls being wavy and flat, respectively. Validation was performed by independently reproducing the results of flow solutions documented in existing databases. Thus, results were obtained for a turbulent flow past a wavy wall and compared to the experimental and DNS data by Ref. 24, provided in the ERCOFTAC Classic Database (case 77). Very good agreement was found for all the statistical quantities available, including the velocity, Reynolds stresses, and pressure field.

The simulations were performed using surface-conformal meshing. The wavy geometry was created by adding an increment $hw(x,y,z)$ to the wall-normal cell coordinates of a plane channel of half-height h, with the walls located at $y=±h$⁠.

A number of grid configurations have been simulated, with the characteristics of the plane-channel mesh $(hw=0)$ listed in Table I, where all quantities are scaled by reference to the target friction Reynolds number. The labels G1–G6 will be used later to identify these cases.

As will transpire, the changes in drag relative to the plane channel are small, pushing the requirements for the spatial resolution to much more stringent levels than for regular DNS. Therefore, particular emphasis is placed on a few simulations performed at the highest tractable resolution within the resources available. These key simulations correspond to the finest grid G6, at a bulk Reynolds number of $Reb=Ubh/ν=6200$ (⁠$Reτ≈360$⁠). Quantitative evidence for the level of refinement is given in Fig. 3 by scaling the grid dimensions by the Kolmogorov length scale in a plane-channel DNS for the mesh G6, showing that the ratio $Δ/η$⁠, where $Δ=ΔxΔyΔz3$⁠, remains lower than unity throughout the channel.

The wavy mesh is generated by adding an increment to the wall-normal coordinate of each and every cell of the plane channel. Two types of wavy geometry are considered herein: one with both walls wavy and the other with one wall flat, as shown in Fig. 4. In the former, shown in Fig. 4(a), both walls are in phase, i.e., yielding a constant passage height of 2h along the entire channel. In the latter, shown in Fig. 4(b), the local height varies from 2h − $Aw$ to $2h+Aw,$ where 2h is the mean channel height, and $Aw$ is the amplitude of the sinusoidal wall undulations.

The skewed wavy channel can be simulated in two ways: the grid can be aligned with the wave or with the main flow direction. Although the two are physically equivalent, keeping the wavy boundary aligned with the numerical box and skewing the flow at an angle to the grid allow greater flexibility. Specifically, the main advantages of this approach are as follows:

1. if the crests are aligned with the z-direction, this direction becomes statistically homogeneous;

2. the post-processing is significantly eased; and

3. this option allows continuous variations of the domain extent in the z-direction without affecting the periodicity boundary conditions applied to the z-direction boundaries.

However, a disadvantage of this option is that it increases the problem size since there is no longer an alignment between the x–z coordinates and the streamwise and spanwise directions, respectively (cf. Fig. 5), necessitating both the domain sizes and resolution to be increased in the two wall-parallel (x–z) directions. This is in contrast to the usual practice of increasing the spanwise resolution whilst decreasing the spanwise domain width. Nevertheless, with the flow inclined at an angle to the grid, longer structures can be captured, as the flow traverses diagonally, thus mitigating the requirements for large domain sizes. Additionally, the difference between the two orientation strategies is shown to be small in Sec. III C.

An approximately constant flow rate across the channel was maintained by iteratively updating the two orthogonal pressure gradients, P_x and P_z, implemented as explicit body forces in the momentum equations so that the bulk velocity is close to unity in the streamwise (⁠$x̌$⁠) direction and zero in the spanwise (⁠$ž$⁠) direction. The main reasons for using a constant flow rate, rather than a constant pressure gradient to drive the flow, is that tests revealed that the duration of the transient is thereby reduced, and that the angle between the mesh and the main flow direction is then prescribed unambiguously. The data show that the target streamwise bulk velocity is maintained within an error lower than 0.001%. Throughout the discussion to follow, only the incremental part of the pressure is reported, relative to a pressure-reference value located at one of the corners of the computational box.

Quantities expressed in the frame of reference of the flow will be denoted with an overlaying inverted circumflex accent (⁠$⋅̌$⁠), as in Fig. 5, although this notation will be omitted later when not needed. Unless stated otherwise, all physical interpretations are given in the frame of reference aligned with the flow.

All statistical quantities can be averaged in the homogeneous direction parallel to the wave crests and troughs. This groove-wise-averaging procedure is significantly eased by the choice made in Sec. II C 2 of forcing the flow at an angle $θ$ to the numerical grid. Phase-averaging is also performed when multiple waves are included within the domain of solution. Furthermore, in the case of a wavy channel with constant wall separation, both boundaries are statistically equivalent, which allows doubling the amount of data included in the phase-averaging by shifting one of the walls by half a period and then taking advantage of the symmetry to average over both walls. Thus, any time- and phase-averaged quantity $q¯$ only depends upon the phase location $x/λ$ and the wall-normal location y, reducing the dependence to $q¯(x/λ,y)$⁠.

Depending on the objective of the analysis, two types of statistical decomposition of the mean turbulence properties can be considered. The first decomposition is relevant to studying how the flow properties vary in phase

where $q¯$ is any time- and phase-averaged quantity, $Q(y)=∫01q¯(x/λ,y+yw)d(x/λ)$⁠, $yw=Awsin(2πx/λ)$ is the wall-normal location of the wall, and $q̃$ is the phase-varying part of the mean field. The second approach lays emphasis on the action of the wavy boundaries relative to the plane-channel flow,

where $Q0$ is the baseline plane-channel-flow value, and $q^$ is the difference to the plane-channel-flow solution. These decompositions will be referred to as “phase-integrated” (⁠$Q$⁠), “phase-varying” (⁠$q̃$⁠), and “difference” (⁠$q^$⁠).

A drag coefficient is defined for each contribution to the drag as

where $⋆$ identifies the contribution (e.g., “f” for friction), $F¯⋆$ is the mean force exerted on the walls opposing the flow direction, and $Ub$ is the bulk velocity.

As mentioned in Sec. II C 2, in all the cases presented, the flow is driven at approximately constant flow rate so that $∥Ub∥≈1$⁠, via the imposition of a spatially constant (vectorial) pressure gradient in order to balance the total drag force. Given a unit bulk velocity, the correct Reynolds number is set by prescribing the appropriate viscosity. The resulting pressure force driving the flow is

where $Px̌$ is the projection of the pressure gradient onto the flow direction: $Px̌=Pxcosθ+Pzsinθ$⁠.

For the wavy channel, the total drag force is composed of two contributions: friction and pressure drag. The friction and pressure forces were integrated on the wavy surface and projected onto the flow direction, yielding the drag coefficients D_f and D_p, respectively.

Since only pressure and friction forces act on the walls, the total drag is

Simulations were performed at $Reτ≈360$ for various configurations, grid resolutions, and domain sizes. This Reynolds number was chosen so that the cost of the simulation remained tractable, and that the ratio between the height of the wave and the channel height $Aw/h$ was kept relatively modest. The corresponding parameters are given in Table II, along with the drag levels calculated, as explained in Sec. II E. The net drag variation is evaluated in two ways: from the imposed pressure gradient and from the sum of the surface-integrated pressure force and friction-drag force. As expressed by Eq. (4), the two should be identical. However, they differ very slightly in the simulations, and this is expressed by the column “imbal.” in Table II. In all simulations, the imbalance of the forces is regarded as negligible. By way of contrast with previous DNS studies reporting this value, the force imbalance by Yoon et al.⁴¹ was of about 3%–4%.

Previous studies on wavy channels often had only a single wavy wall or featured varicose wall undulations.^{10,17,23,25,26,43} A distinct feature of the present configuration is that the passage height is constant at any phase location. Thus, comparing various wall configurations is interesting in order to ensure that the distance between the two walls is large enough so that the upper wavy wall does not influence significantly the flow above the lower one. This is addressed by comparing simulations for channels with one or both walls wavy. The relevant entries in Tables II and III are G2W1 and G2W1f, and the corresponding plane-channel baseline G2P1. The results in Table III show that the friction on the lower (wavy) boundary of G2W1f has a lower value than for G2W1, whilst the friction on the upper (flat) surface of G2W1f is slightly increased with respect to the baseline drag. The latter is a consequence of the control strategy of keeping the flow rate constant. The drag coefficients reported in Table III show that the drag relative to the baseline level (in this case, drag increase) is about half of that of the case with both walls wavy: the drag increase per wavy wall in G2W1 is 1.58%, compared to 1.54% for the wavy-flat case G2W1f. This supports the assertion that the two configurations are close to each other in terms of the processes effective at each wavy wall separately.

In Table II, TDR represents the drag relative to that of the plane channel. Physically, the latter is unique regardless of the angle between the flow and the mesh. However, this is not so computationally, owing to variations in solution domain and numerical errors, including the finite time of averaging, grid resolution, and the finite size and orientation of the periodic domain which does not allow very long structures to be captured. As the angle between the main flow direction and the grid varies, the length of the longest structure allowed to exist within the periodic boundaries changes. At the same time, the grid resolution is also altered in the flow-oriented directions. This results, in a plane channel, in a slight dependence of the drag on the flow direction. By way of an example of the influence of the domain size on the drag, for a flow aligned with the mesh, Ricco and Wu³² observed, at $Reτ=200$⁠, that changing the streamwise extent from $4πh$ to 21h whilst keeping roughly the same spatial resolution resulted in a drag increase of 0.6%. This difference in drag is already large enough to be of the same order of magnitude as the variations of the total drag level in the cases considered herein. Therefore, a careful assessment of key computational parameters is required, as detailed below.

First, the influence of the domain size is considered for the plane-channel flow. This was investigated by means of highly resolved simulations, with constant resolution, but varying domain sizes. For this purpose, grids G3 and G4 are chosen to have the same spatial resolution $Δx+=Δz+=2.5$ and $Δy+$ ranging from 0.4 at the wall up to 8.5 in the centre, but the domain size of G4 is one half of that of G3 in the x-z directions (cf. Table I). The largest domain considered is L_x = L_z = 15h, which represents about 5400 wall units—i.e., longer than the commonly chosen value of $8πh+=O(4500)$ at $Reτ=O(180)$⁠. Table II shows that the relative difference in drag between the larger domain (grid G3) and the smaller (grid G4) is lower than 0.05%, both at an angle of $0°$ (G3P1 and G4P1) and $45°$ (G3P2 and G4P2).

Next, the impact of resolution is quantified, still for the case of a plane channel. To this end, simulations were run with the domain size kept constant, as in G3, but with a mesh coarser by a factor of 2 in the wall-parallel and wall-normal directions independently. Then, the drag levels for $θ=0°$ and $45°$ configurations were compared. As expected for a consistent discretisation, the difference between the two physically equivalent configurations (⁠$θ=0°,45°$⁠) decreases as the resolution is increased. However, this difference remains significant at about 0.7% for grid G3, despite the mesh being fine with respect to usual DNS standards. An important observation is that increasing the number of cells in the wall-normal direction has little impact on the total drag, whereas increasing the resolution in the wall-parallel directions reduces significantly the difference in friction between $θ=0°$ and $θ=45°$⁠. Drag differences with respect to the angle of the flow were observed to vary monotonically, the minimum drag coefficient being found at $θ=0°$ and the maximum at $θ=45°$⁠.

A possible approach to reducing the error in the predicted drag-reduction level is to evaluate it by reference to a plane-channel flow simulation at exactly the same spatial resolution, domain size, and flow angle. This is the approach preferred here, in light of the work of Gatti and Quadrio,^15,16 who observed some cancellation of the systematic bias associated with the domain size. Thus, for example, the baseline drag for G2W1 is G2P1 (both at $θ=52°$⁠), whereas for G2W1bis (⁠$θ=70°$⁠), the baseline is taken to be G2P2 (also at $θ=70°$⁠). However, even this elaborate approach may not suffice, in the face of the small drag-reduction margin, to remove uncertainties associated with the numerical aspects that contribute to the total error. A factor that may also be influential is the distortion of the cells fitted to the wavy boundary. As shown in Table IV, changes in the total drag from grid G2 to the finest grid G6 are not independent of the wave geometry.

The variation of the total-drag reduction with grid refinement was studied in greater detail for the most promising case, namely, $(Aw+,θ,λ+)=(18,70°,918)$⁠. Since the drag was found to be mainly sensitive to the wall-parallel resolutions, only $Δx$ and $Δz$ are used as indicators of the grid refinement. More precisely, for such fine resolutions, and in light of the results shown in Sec. III C 1, it is reasonable to assume the accuracy to be mostly sensitive to the spacing in the spanwise direction. Since the flow is not aligned with the grid, the measure of the grid spacing used is the spanwise-projected cell spacing, defined as $Δž=min{Δz/cosθ,Δx/sinθ}$⁠. The main outcome of the study, shown in Fig. 6, is that the total-drag reduction predicted decreases as the mesh is refined, with a quadratic dependence on the spanwise mesh spacing, consistent with the order of accuracy of the spatial discretisation. The drag-reduction value appears to tend, asymptotically, to a positive value of 0.6%. Additionally, the ratios of the friction- and pressure-drag coefficients (D_f and D_p, respectively), relative to the total D_tot, remain constant for grids G2 and G6; thus, the value of D_p is 1.314% that of D_tot for the coarsest mesh G2W1A2, compared with 1.315% for the finest mesh G6W1A2. Therefore, whilst the share between the pressure and friction drag is grid-independent, the absolute value of the total drag is a very sensitive quantity that requires substantial computational efforts to be accurately predicted.

One additional simulation, physically equivalent to G6W1A2, was undertaken, but in a computational box aligned with the main flow direction, the waves being at an angle to the grid. Two wavelengths were included in the streamwise and spanwise directions, and the grid resolution is similar to that of the grid G6. The result, shown in Fig. 6 by way of the red cross, is in line with the simulations of the skewed configuration.

Although Fig. 6 may give the impression that there is a large difference between the coarsest- and finest-grid simulations, it needs to be emphasised that the grid-dependence of the absolute drag level and other flow properties is very modest. To substantiate this claim, results for velocity and Reynolds stresses obtained with the coarsest and finest grid are compared in Figs. 7 and 8 for the same flow configuration already studied in Fig. 6 (⁠$Aw+=18$⁠, $θ=70°$⁠, and $λ+=918$⁠). Figure 7(a) shows that the streamwise-velocity profiles are virtually indistinguishable. Differences in the phase-integrated Reynolds stresses, shown in Fig. 7(b), are slightly larger, but also very small. Finally, Fig. 8 demonstrates that the corresponding profiles of phase-averaged transverse (Stokes) strain are virtually indistinguishable.

As the flow travels past the skewed wavy wall, it accelerates on the windward side and then decelerates on the leeward side, a behaviour linked to the pressure being a minimum above the crest and a maximum in the trough region, as shown in Fig. 9. Because the wave is at an angle to the main flow direction, this pressure variation in phase gives rise to a pressure gradient both in the streamwise and spanwise directions. The latter gradient generates a spanwise motion, shown in Fig. 10, which is asymmetric in phase and penetrates quite far into the boundary layer above $y+≈100$⁠. However, as previously mentioned, it is not the velocity but the shear strain which is important with respect to the emulation of the Stokes layer, and which dictates the orientation of the near-wall streaks. In the work by Touber and Leschziner,³⁹ this orientation was observed to vary in time, with strong reduction in turbulence during the reorientation phase. For the wavy wall, the phase-modulation of the shear strain occurs in space and also results in a reorientation following the shear-strain field, as well as in a weakening of the streaks, as shown in Fig. 11 at a constant distance from the wall around $y+≈10$⁠. However, the effect is far less pronounced than observed by Touber and Leschziner³⁹ because the forcing amplitude is much smaller in the present case.

In the work by Viotti et al.,⁴⁰ the reduction in skin friction was observed to increase linearly with forcing amplitudes up to $ASSL+≈5$ and then to increase at a slower rate. Such an observation is not consistent with the expected symmetry of the problem around $ASSL+=0$⁠, which would imply a zero value of the first derivative of the drag reduction as a function of the amplitude. For the wavy wall, the dependence of the friction-drag reduction (FDR) on the wave slope, shown in Fig. 12(a), is compatible with a symmetry condition, thus indicating that this is a local effect taking place below $2Aw/λ≈0.04$⁠: for very small actuation amplitudes, the FDR does seem to exhibit a quadratic behaviour, which then becomes close to linear. The amplitude of the spanwise shear strain at $2Aw/λ≈0.04$ corresponds to that of a SSL with a forcing amplitude of about $ASSL+≈1.1$⁠, thereby corroborating the observations of Viotti et al.,⁴⁰ since the smallest amplitude they considered was $ASSL+=1$⁠, which occurs at the intersection between the quadratic and linear behaviours for the present key simulations.

The decrease in $λx$ from G6W1* (⁠$λ+=918$ and $λx+≈2700$⁠) to G6W2* (⁠$λ+=612$ and $λx+≈1800$⁠) results in a reduced effectiveness of the wavy wall at the same wave slope—i.e., the same equivalent forcing amplitude $ASSL+$⁠. This trend contrasts with the drag-reduction trend of the SSL, which features an optimum around $λx+≈1250$⁠. Therefore, there exists some unfavourable mechanism limiting the drag reduction achievable by the wavy wall. Such considerations will be discussed in Sec. IV C.

The present approach of using a skewed wavy wall is based on the assumption that similar longitudinal patterns of the transverse shear strain, whether created by the SSL or the wavy wall, will lead to the same effects on turbulence and hence lead to some turbulent-drag reduction. In the case of a temporal Stokes layer, Touber and Leschziner³⁹ showed that the wall oscillations led to a bimodal partial decay, reformation, and reorientation of the streaks, a behaviour dictated by the unsteady Stokes strain in the buffer layer. If this mechanism is indeed intimately linked to the friction-drag reduction, the relevant forcing quantity is the resulting shear strain, rather than the forcing velocity itself. It is this key element that allows for a passive surface, with no slip at the solid wall, to emulate the actuation by in-plane wall oscillations, through the action of a transverse pressure gradient that generates an equivalent shear-strain field. Importantly, however, any differences between the shape of the Stokes-strain profiles generated by the SSL and wavy wall may be very influential to the equivalence between the two actuation scenarios. In particular, the wall-normal penetration of the Stokes strain is a critical factor.^2,12 In sub-optimal TSL or SSL actuation (e.g., at too low actuation frequency and/or excessive wave-length), the strain is observed to penetrate into the buffer layer, thus enhancing turbulence generation and counteracting the favourable effects of the Stokes strain in the viscous sublayer. Yakeno et al.⁴² show that, for an oscillating wall, the spanwise shear strain close to the wall counteracts the quasi-streamwise vortical motion, which in turn results in a weakening of the Quadrant 2 events, associated with ejections. Further away from the wall, the presence of spanwise shear strain enhances Quadrant 4 events at sub-optimal actuation, thus increasing turbulent mixing.

In contrast to the Stokes layer, the wavy wall also induces wall-normal forcing that contributes to the shear strain via $∂v¯/∂z$⁠, but this additional effect was observed to be small compared to $∂w¯/∂y$⁠, especially in the near-wall region where the shear is maximum. It follows that there is no significant parasitic effect of the wall-normal velocity on the spanwise forcing, justifying a direct comparison of $∂w¯/∂y$ between wavy-wall and Stokes-layer configurations.

A slight difference between the SSL simulations presented in this section, and those of Viotti et al.,⁴⁰ has been introduced deliberately, in order to maximise the correspondence between the SSL and the wavy-wavy channel configuration. This difference is that the wall forcing on the upper wall is shifted by half a period relative to the lower wall. However, this change does not result in noticeable differences since the thickness of the Stokes layer is much smaller than the channel half-height. The reference SSL considered is for a forcing amplitude of $ASSL+=2$ (based on unactuated friction velocity) and a wavelength close to the optimum at this forcing amplitude, $λx+≈1260$⁠, subject to the assumption that the Reynolds-number change from $Reτ=200$ in the work by Viotti et al.⁴⁰ to the present value of $Reτ≈360$ does not have a significant impact on the optimal wavelength. The simulation was run on a domain $Lx+≈2500$⁠, $Lz+≈1100$⁠, with a grid resolution $Δx+=9.8$⁠, $Δz+=5.8$⁠, and $0.7<Δy+<7.3$⁠. Such a resolution may not be sufficient for an accurate comparison of the drag levels but is acceptable for the comparison of the shear-strain profiles, as shown in Sec. III C 3.

Results for some wavy channels are shown in Fig. 13. The strain profiles demonstrate that the wavy wall emulates reasonably well the shear layer of the SSL. This observation leads to the expectation that the reduction in turbulent skin friction achieved by the wavy wall would be similar and arises from the same physical mechanism as in the SSL. In addition, the amplitude of the phase-wise variations in the shear strain—and hence the corresponding SSL forcing amplitude $ASSL+$—appears to be mostly dictated by the streamwise wave slope $Aw/λx$⁠, as suggested by rescaling the amplitude of the transverse shear strain by this ratio, so as to compensate for the different forcing amplitudes. This implies that for $Aw$ and $λ$ kept constant (same wave shape), increasing $θ$ results in a decrease in the forcing amplitude, at least for angles in the range $θ∈[50°,70°]$⁠.

The similarity between the SSL and the wavy wall, in terms of the friction-reducing transverse strain, is given further weight by a supplementary simulation of the wall-actuated case at the conditions $Reτ≈360$⁠, $ASSL+=1.25$⁠, and $λx+=2700$⁠, approximately corresponding to the flow configuration that yields the maximum drag-reduction level, namely, G6W1A2 (⁠$Aw+=18$⁠, $θ=70°$⁠, and $λ+=918$⁠). The resolution used was $Δx+=4$⁠, $Δy+∈[0.6,4.5]$⁠, and $Δz+=2$⁠, and the domain size was L_x = 7.5h and L_z = 3.2h. The drag-reduction level measured for the SSL simulation is $1.9±0.3%$⁠, which compares favourably with the friction-reduction level in case G6W1A2 of $2.0±0.2%$⁠, as shown in Table II. However, it should be understood that significant differences exist, in terms of the detailed physics, between the wall-actuated flow and the wavy-channel flow, and it is important to point out that the present comparison is limited to the particular flow configuration considered. In other flow configurations, the equivalence shown in Fig. 13 may not be as close, and the previously noted observations of Yakeno et al.⁴² on the importance of the shape of the Stokes-strain profile are pertinent.

As is observed with Stokes layers and, more generally, with most control strategies yielding turbulent-drag reduction, an upward shift in the log-law appears relative to the uncontrolled flow when actual scaling is used (i.e., with the actual friction velocity). This shift is shown in Fig. 14(a) for some key simulations. The corresponding flow configurations are characterised by two height-to-wavelength ratios and two wavelengths. Although there is a noticeable difference between the two profiles for the two wavelengths, it is observed that similar wave slopes yield approximately similar shifts, the configurations at $λ+=918$ featuring a slightly greater shift that implies a higher level of skin-friction reduction (cf. Sec. III E). As the wave slope is linked to the forcing amplitude, the higher this ratio is, the higher the friction-drag reduction is. This close resemblance of the log shift at a given wave slope, but different wavelengths, indicates a limited sensitivity of the friction-drag reduction for the range of wavelengths tested, especially at such a small forcing amplitude.

As far as the Reynolds stresses are concerned, there is a substantial decline in the peak of the streamwise normal stresses—evidence of the weakening of the streaks. The behaviour of the Reynolds stresses is also similar for a given height-to-wavelength ratio, apart from a detrimental increase in the streamwise normal stresses for the shorter wavelength $λ+=612$⁠, starting within the buffer layer around $y+≈20$ and persisting up to $y+≈80$⁠. However, unlike wall-actuated Stokes layers,^1,39,40 which entail a larger forcing amplitude, the decline in the streamwise stress only persists in the present configuration up to $y+≈35$⁠, beyond which it exceeds the baseline value. Similarly, the shear stress is depressed up to about the same wall-normal location, beyond which it also exceeds the baseline level, unlike for the high-forcing-amplitude Stokes layers for which no increase is found.

The difference between the streamwise Reynolds-stress levels for the wavy walls and the baseline case is shown in Fig. 15 for G6W1A2. This brings to light two different regions showing distinct physical features. One region is close to the wall, where a material weakening of the streaks takes place, and another is further away above the trough, featuring enhanced streamwise turbulence intensity. The latter increase is stronger than the reduction above the crest, thus leading to a net increase in the mean streamwise turbulence intensity above $y+≈35$⁠, relative to the baseline [cf. Fig. 14(b)]. Again, the phase-variations along the wavy wall are much larger than in the SSL case where they are negligible.

Despite the similarity between the shear-strain phase variations of wavy walls and Stokes layers, highlighted in Sec. IV A, the effectiveness of the wavy wall is found to be lower. In order to understand the lower performance, two main mechanisms are considered as possible causes of the degradation.

First, beyond the observations made in Sec. IV B, an important difference is that the friction is increased on the windward side of the wave, as shown in Fig. 16. The overall skin-friction reduction arises as a balance between the depressed friction on the leeward side of the wave and the enhanced friction on the windward side, whereas in the SSL case, the friction is decreased at all phases. This variation is associated with the magnitude of the phase-varying streamwise velocity $ũ$⁠, which is greater than $w̃$⁠, whereas the former is almost negligible in the SSL (in optimum actuation conditions). This phase-variation of the mean longitudinal velocity was already identified by Chernyshenko⁹ as the main source of degradation of the performance of the wavy wall relative to that of the SSL.

Second, an additional mechanism, specific to the wavy wall, was revealed by the numerical calculations. As shown in Fig. 17(a), there exists a zone of intense production of turbulent kinetic energy $Πk=−ui′uj′¯∂u¯i/∂xj$ above the leeward side of the wave, reflecting a deeper penetration of the disturbance arising from the wavy wall into the boundary layer. As shown in Fig. 17(b), the increase in production relative to the baseline is quickly followed, in phase, by an increase in energy dissipation $𝜖k=−ν∂ui′/∂xj∂ui′/∂xj¯$⁠, at about the same wall-normal location. This phenomenon strengthens as the wave amplitude increases, and is stronger for G6W2* than for G6W1* at similar wave slopes.

An interesting question is to what extent the flow properties are sensitive to the Reynolds number, whilst the shape of the wall is kept constant in wall units. This has been investigated by reference to the flow configurations listed in Table V, the maximum friction Reynolds number being 1000.

Although net-drag-reduction levels of about 1%–2% were observed for the three Reynolds numbers tested, the value is not quantifiable with any greater degree of precision because it is extremely sensitive to various numerical issues, as demonstrated in Sec. III C through a grid-convergence study at $Reτ=360$⁠. However, it was also shown that the relative error of first- and second-order statistics is far less sensitive to the grid spacing and could therefore be compared with much greater confidence.

The phase-integrated streamwise-velocity and Reynolds-stress profiles are shown in Fig. 18. Figure 18(a) indicates, more clearly than for lower Reynolds-number values, that there is an upward shift in the log-law, associated with the friction reduction. In contrast, Fig. 14(a) shows, for $Reτ≈360$⁠, a change in slope. However, this needs to be considered against the background of the very short log-law portion combined with the fact that the large wave height, $Aw+≈20–30$⁠, provokes significant phase-varying streamwise-velocity perturbations that extend well into the log-law region, thus likely to cause deviations from the usual undisturbed channel-flow log-law.

Despite the wavy geometry in wall units not being exactly the same across the three flows in Table V, and the mesh being somewhat coarser for the $Reτ=1000$ case, the shear-strain profiles, scaled in wall units, are almost identical as demonstrated in Fig. 19. At $Reτ=180$⁠, a wave height of $Aw+=20$ represents a ratio $Aw/h$ of about 11%, which is significant, while it decreases to 2% for $Reτ=1000$⁠. Consequently, the wave height $Aw$ appears to be small enough relative to the channel height so that even at the lowest Reynolds number where the ratio $Aw/h$ is maximum, the distance between the walls remains sufficient to avoid an interference between the two solid boundaries, adding support to the findings reported in Sec. III B.

In the present study, skewed wavy-wall channels have been investigated by means of direct numerical simulation as a potential passive open-loop drag-reduction device. The spanwise shear-strain profiles generated by the wavy wall were shown to resemble closely those of the well-established method of drag reduction by in-plane wall motions and to depend only weakly on the Reynolds number when expressed in wall units. Various wavy-wall geometries for different combinations of flow angle $θ$⁠, wavelength $λ$⁠, and height $Aw$ were explored with the aim of seeking a configuration that minimises the total drag relative to a plane wall.

Unlike for the Stokes layer, there is no actuation power required, but the skin-friction reduction is militated against by the pressure drag arising from the wavy geometry. This drag increases quadratically with wave height, rapidly exceeding the friction-drag reduction beyond a modest wave height. Consequently, the cumulative effect on the drag is small. The corresponding accuracy requirements of the simulations, in terms of grid density and integration time, are, therefore, far beyond those generally adopted in DNS of channel flows, making the cost of optimising the wavy-wall parameters prohibitive. Even if relative changes in drag are quantified by reference to drag levels of a baseline plane channel, simulated with similar grid spacing, domain size, and flow-to-mesh angle, the net drag-reduction level is subject to a not-insignificant relative error.

Despite the above qualifications, a net drag-reduction value of about 0.7% (made up of a 2%-friction reduction and a pressure-drag penalty of 1.3%) was estimated for the configuration $Aw+≈20$⁠, $θ=70°$⁠, and $λ+≈920$⁠, at $Reτ≈360$⁠, at the finest grid resolution, with indications that this performance would drop slightly below 0.6% for an asymptotically fine mesh.

The generation of a significant phase variation of the mean longitudinal velocity, proposed in earlier studies as a mechanism accounting for the degradation of the performance of the wavy wall relative to the steady Stokes layer, was augmented by the identification of a new mechanism consisting of the intense localised production of turbulence kinetic energy above the leeward side of the wave.

This study was funded by Innovate UK (Technology Strategy Board), as part of the ALFET project, Project Reference No. 113022. The authors are grateful to the UK Turbulence Consortium (UKTC) for providing computational resources on the national supercomputing facility ARCHER under the EPSRC Grant No. EP/L000261/1. Access to Imperial College High Performance Computing Service, doi: 10.14469/hpc/2232, is also acknowledged.

The objective of finding the wall shape $hw$ (defined in Fig. 2) that maximises the drag reduction can be formulated as an optimisation problem, wherein the objective function to be minimised is the total drag, and the parameter space is three-dimensional, composed of the wave amplitude, wavelength, and flow-to-crest angle. In the work by Chernyshenko,⁹ a proxy for the objective function was used, based upon the assumption that some equivalence can be found between the flow distortions provoked by the wavy wall and those induced by the spatially varying steady wall actuation. The major ingredients and results of this approach are briefly summarised.

Following the work of Viotti et al.,⁴⁰ the boundary-layer equations are linearised around a linear profile for the mean streamwise velocity, resulting in ordinary differential equations for the perturbation fields $(u,v,w,p),$ expressed in wall units. The latter are solved both for the SSL and the wavy wall.

• In the SSL case, the streamwise and wall-normal momentum equations decouple from the spanwise direction so that the streamwise velocity is unchanged, and only distortions of the spanwise component are present. The spanwise-velocity perturbation satisfies the Airy equation, the solution of which is given by

  $wSSL=ASSLAi(0)ReikxxAi−iỹe−i4π/3,$

  (A1)

  where $R$ denotes the real part, $ỹ=kx1/3y$⁠, and Ai is the Airy function of the first kind.
• In the wavy-wall case, however, there are also perturbations in the streamwise velocity, arising from a combined effect of the wall-normal velocity induced by the spanwise inhomogeneity of the spanwise velocity and the longitudinal pressure gradient induced by the wall. Assuming a perturbation of the form $(u,v,w,p)=(û(y),v^(y),ŵ(y),p^)ei(kxx+kzz)$⁠, the linearised boun-dary-layer equations become

  $ikxyû+v^=−ikxp^+d2ûdy2,ikxyŵ=−ikzp^+d2ŵdy2,ikxû+dv^dy+ikzŵ=0.$

  (A2)

  Equations (A2) are then transformed into ordinary differential equations, which are solved numerically.

The L²-norm of the difference in the spanwise shear-strain profiles between a SSL of given forcing amplitude, A_SSL (denoted $Ŵ$ by Chernyshenko⁹), and the wavy wall, is minimised over $p^$⁠, resulting in the ratio of $p^$ to A_SSL yielding close shear-strain profiles. However, owing to the nature of the analysis, no explicit formula was obtained for the correspondent height of the wavy wall $Aw$⁠.

Next, it is assumed that the generation of a similar spanwise shear-strain layer causes the same level of weakening of the near-wall turbulence as in the wall-actuated case. In both cases, the net energy savings are decomposed as $Pnet=Psav+Preq$⁠, where $Psav$ is the power saved from the turbulence weakening, and $Preq$ is the power required to do so.

• Given the above-mentioned condition on the wave-induced pressure gradient, the value of $Psav$ is taken from a polynomial interpolation of the DNS results of Viotti et al.⁴⁰ for the SSL at various forcing amplitudes and wavelengths.

• In the SSL case, the power required $Preq$ is evaluated from the dissipation arising from the mean-velocity gradient induced by the control effects

  $PreqSSL=−Cf2kx+(1−Psav)1/3∫0∞12dwSSLdỹ2dỹ,$

  (A3)

  where the value of the friction coefficient is evaluated from the empirical relation $Cf=0.0336Reτ−0.273$ given by Pope.²⁸ Thus, in the case of a wavy wall, additional dissipation arises not only from the generation of the spanwise shear strain but also from the streamwise perturbation of the velocity. For a given wavelength $λx$ and spanwise pressure gradient (i.e., fixed $ASSL$⁠), this effect can be minimised by varying the angle of the wave through $λz$⁠, and the calculation of Chernyshenko⁹ yields an optimal angle $θopt≈52°$⁠. At $θ=θopt$⁠, the ratio

  $r=∫0∞dûdỹ2+dŵdỹ2dỹ∫0∞dŵdỹ2dỹ$

  (A4)

  between the total additional dissipation and that due to the generation of the spanwise motion is $rmin≈5.8$⁠, implying that the energy dissipation arising from the streamwise-velocity perturbation is almost 5 times as large as that associated with the generation of the spanwise shear strain.

Finally, at the optimal angle $θopt≈52°$⁠, the power required for the wavy wall to generate the spanwise motion is evaluated as $Preq=rminPreqSSL$—i.e., that of the equivalent SSL $PreqSSL$ multiplied by the optimal ratio $rmin$ accounting for the additional streamwise distortion. The resulting net energy saving $Pnet=Psav+rminPreqSSL$ is then maximised with respect to the forcing amplitude of the equivalent SSL and the associated wavelength, yielding $ASSL, opt+=2$ and $λx,opt+=1520$⁠.

The importance of the streamwise phase-varying velocity, relative to that of the induced spanwise velocity, was confirmed by the DNS results. Indeed, a simulation was carried out at $Reτ≈180$⁠, for $λx+≈1510$⁠, $θ=52°$⁠, and $Aw+≈20$⁠, resulting in a ratio

which is close to the estimate given in the analysis of Chernyshenko⁹ of $rmin≈5.8$⁠.

Data associated with this article can be accessed via the research website of S. Chernyshenko at https://www.imperial.ac.uk/aeronautics/fluiddynamics/ChernyshenkoResearch/misc.php.