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

Direct numerical simulations of fully-developed turbulent channel flows with wavy walls are undertaken. The wavy walls, skewed with respect to the mean flow direction, are introduced as a means of emulating a Spatial Stokes Layer (SSL) induced by in-plane wall motion. The transverse shear strain above the wavy wall is shown to be similar to that of a SSL, thereby affecting the turbulent flow, and leading to a reduction in the turbulent skin-friction drag. However, some important differences with respect to the SSL case are brought to light too. In particular, the phase variations of the turbulent properties are accentuated and, unlike in the SSL case, there is a region of increased turbulence production over a portion of the wall, above the leeward side of the wave, thus giving rise to a local increase in dissipation. The pressureand friction-drag levels are carefully quantified for various flow configurations, exhibiting a combined maximum overall-drag reduction of about 0.5%. The friction-drag reduction is shown to behave approximately quadratically for small wave slopes and then linearly for higher slopes, whilst the pressure-drag penalty increases quadratically. The transverse shear-strain layer is shown to be approximately Reynolds-number independent when the wave geometry is scaled in wall units. ∗ s.ghebali14@imperial.ac.uk † s.chernyshenko@imperial.ac.uk ‡ mike.leschziner@imperial.ac.uk


I. INTRODUCTION
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 [8], 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 in [4,21,27], but despite attempts to optimise the geometry, Bannier [4] 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. [20] 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 w w (x, t) = A sin(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. [3] in pipe flows, for which drag-reduction levels of up to 33% were achieved, while Bird et al. [7] designed a Kagome-lattice-based actuator for boundary layers, achieving a drag-reduction level of around 10%. 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 w w = A SSL sin(2π x/λ x ), as shown in figure 1. This method was studied by Viotti et al. [40] by means of DNS for various forcing amplitudes A SSL and wavelengths λ x . The 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 [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 figure 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 span-wise shear strain, somewhat away from the wall, that will weaken turbulence in a similar manner to that effected 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 wavychannel flows with the aim of examining the degree of Stokes-layer emulation achieved, and the degree to which the drag is reduced relative 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 figure 2, will be brought to light, including the impact on the near-wall turbulence.

A. Overall strategy
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 [9] 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 A + SSL = 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.

B. Computational simulations
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 in [13]. In the fractional-step procedure, the nonsolenoidal intermediate velocity field is projected onto the solenoidal space by solving a pressure-Poisson equation. The latter is solved by Successive Line Over-Relaxation [18, p. 510], the convergence of which is accelerated using a multigrid algorithm by [22]. 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 −9 in all the simulations presented.
Stable velocity-pressure coupling is ensured by use of the Rhie-and-Chow interpolation [31], preventing odd-even oscillations.
Simulations were initiated with a laminar Poiseuille profile to which random perturbations are 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][34][35][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 [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.
C. Spatial discretisation of the problem

Simulation of a wavy channel
The simulations were performed using surface-conformal meshing. The wavy geometry was created by adding an increment h w (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 (h w = 0) listed in table I, where all quantities are scaled by reference to the target friction Reynolds number. The labels G1 to 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 Re b = U b h/ν = 6200 (Re τ ≈ 360). Quantitative evidence for the level of refinement is given in figure 3     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. figure 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 section 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 shows 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 figure 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.

D. Flow decomposition
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 section 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 a doubling of the 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) = 1 0 q(x/λ, y + y w ) d(x/λ), y w = A w sin(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 Q 0 is the baseline plane-channel-flow value, and q is the difference to the planechannel-flow solution. These decompositions will be referred to as 'phase-integrated' (Q), 'phase-varying' ( q), and 'difference' ( q).

E. Calculation of the drag contributions
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 U b the bulk velocity.
As mentioned in section II C 2, in all the cases presented, the flow is driven at approximately constant flow rate, so that U b ≈ 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 = P x cos θ + P z sin θ.
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: III. SIMULATIONS

A. Overview of simulations
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 A w /h was kept relatively modest. The corresponding parameters are given in table II, along with the drag levels calculated, as explained in section 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 equation (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 in [41] was of about 3-4%.

B. Influence of the upper wall
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     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

Drag-reduction level
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 . 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 section 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 figure 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.5%. 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 figure 6 by way of the red cross, is in line with the simulations of the skewed configuration. The baseline drag level is taken at the same angle, except for G3 where the baseline drag is found by interpolation between θ = 0 • and θ = 90 • cases. Error bars: 90% confidence interval of the time-averaging error calculated using the method of batch means and batch correlations by Russo and Luchini [37], and assuming that the baseline and controlled drag levels are independent variables. 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 direction. The latter gradient generates a spanwise motion, shown in figure 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 [39], 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 figure 11 at a constant distance from the wall around y + ≈ 10. However, the effect is far less pronounced than observed in [39] because the forcing amplitude is much smaller in the present case.

E. Friction-drag reduction and pressure-drag increase
In [40], the reduction in skin friction was observed to increase linearly with forcing amplitudes up to A + SSL ≈ 5, and then to increase at a slower rate. Such an observation is not consistent with the expected symmetry of the problem around A + SSL = 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 figure 12 (a), is compatible with a symmetry condition, thus indicating that this is a local effect taking place below 2A w /λ ≈ 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 2A w /λ ≈ 0.04 corresponds to that of a SSL with a forcing amplitude of about A + SSL ≈ 1.1, thereby corroborating the observations of [40], since the smallest amplitude they considered was A + SSL = 1, which occurs at the intersection between the quadratic and linear behaviour for the present key simulations.
The decrease in λ x from G6W1* (λ + = 918, λ + x ≈ 2700) to G6W2* (λ + = 612, λ + x ≈ 1800) results in a reduced effectiveness of the wavy wall at the same wave slope, i.e. the same equivalent forcing amplitude A + SSL . 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 section IV C.

A. Shear strain
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 [39] 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 SSL or TSL 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. [42] 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 con-tributes 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 [40], 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 A + SSL = 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 [40], to the present value of Re τ ≈ 360 does not have a significant impact on the optimal wavelength. The simulation was run on a domain L + x ≈ 2500, L + z ≈ 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 section III C 3.
Results for some wavy channels are shown in figure 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 A + SSL -appears to be mostly dictated by the streamwise wave slope A w /λ x , as suggested by rescaling the amplitude of the transverse shear strain by this ratio, so as to compensate for the different forcing amplitude. This implies that, for A w and λ kept constant (same wave shape), increasing θ results in a decrease in the forcing amplitude, at least for angles in the range θ ∈ [50 • , 70 • ].
FIG. 13. Comparison of ∂w + /∂y + scaled by the streamwise-projected wave slope A w /λ x for various flow configurations. The SSL flow is for a wavelenth λ x ≈ 1250 and forcing amplitude A + SSL = 2 based on the unactuated friction velocity. The value of ∂w + /∂y + is based on the actual friction velocity, and the scaling factor A w /λ x based on the parameters of G2W1. In all cases, the profiles were shifted in the wall-normal direction in order to match the wave height of G2W1.
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, A + SSL = 1.25, and λ + x = 2700, approximately corresponding to the flow configuration that yields the maximum drag-reduction level, namely G6W1A2 (A + w = 18, θ = 70 • , λ + = 918). The resolution used was ∆x + = 4, ∆y + ∈ [0.6, 4.5], and ∆z + = 2, and the domain size was L x = 7.5h, 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 figure 13 may not be as close, and the previously noted observations of Yakeno et al. [42] on the importance of the shape of the Stokes-strain profile are pertinent.

B. Streamwise velocity and Reynolds stresses
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 uncon-trolled flow when actual scaling is used (i.e. with the actual friction velocity). This shift is shown in figure 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. section 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 wavelength, 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 figure 15 for G6W2A4. 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. figure 14b). Again, the phase-variations along the wavy wall are much larger than in the SSL case where they are negligible.

C. Detrimental effects of the wavy wall
Despite the similarity between the shear-strain phase variations of wavy walls and Stokes layers, highlighted in section 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 section IV B, an important difference is that the friction is increased on the windward side of the wave, as shown in figure 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 u, 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 in [9] 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 figure 17(a), there exists a zone of intense production of turbulent kinetic energy Π k = − u ′ i u ′ j ∂u i /∂x j 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 figure 17(b), the increase in production relative to the baseline is quickly followed, in phase, by an increase in energy dissipation ǫ k = −ν ∂u ′ i /∂x j ∂u ′ i /∂x j , 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. 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 section 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. it decreases to 2% for Re τ = 1000. Consequently, the wave height A w appears to be small enough relative to the channel height, so that, even at the lowest Reynolds number where the ratio A w /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 section III B.

V. CONCLUSIONS
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 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 A w 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 A + w ≈ 20, θ = 70 • , λ + ≈ 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. The objective of finding the wall shape h w (defined in figure 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 [9], 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. [40], 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 where ℜ denotes the real part,ỹ = k 1/3 x y, 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) e i(kxx+kzz) , the linearised boundary-layer equations become Equations (A2) are then transformed into ordinary differential equations, which are solved numerically.
The L 2 -norm of the difference in the spanwise shear-strain profiles between a SSL of given forcing amplitude, A SSL (denotedŴ in [9]), and the wavy wall, is minimised overp, resulting in the ratio ofp 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 A w .
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 P net = P sav + P req , where P sav is the power saved from the turbulence weakening, and P req is the power required to do so.
• Given the above-mentioned condition on the wave-induced pressure gradient, the value of P sav is taken from a polynomial interpolation of the DNS results of Viotti et al. [40] for the SSL at various forcing amplitudes and wavelengths.
• In the SSL case, the power required P req is evaluated from the dissipation arising from the mean-velocity gradient induced by the control effects where the value of the friction coefficient is evaluated from the empirical relation C f = 0.0336Re −0.273 τ given by Pope [28]. 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 A SSL ), this effect can be minimised by varying the angle of the wave through λ z , and the calculation of Chernyshenko [9] yields an optimal angle θ opt ≈ 52 • . At θ = θ opt , the ratio r = between the total additional dissipation and that due to the generation of the spanwise motion is r min ≈ 5.8, implying that the energy dissipation arising from the streamwisevelocity 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 P req = r min P SSL req , i.e. that of the equivalent SSL P SSL req multiplied by the optimal ratio r min accounting for the additional streamwise distortion. The resulting net energy saving P net = P sav + r min P SSL req is then maximised with respect to the forcing amplitude of the equivalent SSL and the associated wavelength, yielding A + SSL,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 A + w ≈ 20, resulting in a ratio which is close to the estimate given in the analysis of Chernyshenko [9] of r min ≈ 5.8.