Turbulent skin-friction reduction by wavy surfaces

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. The pressure- and 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. 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 transverse shear-strain layer is shown to be approximately Reynolds-number independent when the wave geometry is scaled in wall units.


I. INTRODUCTION
In the context of greener and more cost-effective aviation, industrial and academic researchers have proposed and studied a wide range of 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 [1], 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% [2][3][4][5]. 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 [6][7][8], but despite attempts to optimise the geometry, Bannier [8] showed that conventional (straight) riblets appear to be as effective.
On the active-control front, based upon the work of Jung et al. [9] 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 frictiondrag-reduction levels of up to about 45% [10,11] at low Reynolds numbers, the effectiveness being observed to reduce at higher Reynolds numbers [12][13][14]. An experimental confirmation of the numerical results for the travelling wave was undertaken by Auteri et al. [15] in pipe flows, for which drag-reduction levels of up to 33% were achieved, while Bird et al. [16] designed a Kagomelattice-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. [17] by means of DNS for various forcing amplitudes A SSL 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 arXiv:1705.01989v2 [physics.flu-dyn] 16 May 2017 of emulating spanwise in-plane wall motions, as proposed in [18]. 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 spanwise 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 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 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 'carpetbombing' 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 [18] to narrow down the exploration range within which the drag reduction might be maximised. A preliminary study of this type was undertaken by Chernyshenko who found an estimate for the streamwise-projected wavelength of the wave λ + x ≈ 1500, and the flow angle θ ≈ 52 • , but who did not provide an estimate for the height of the wave. Rather, a condition was given for the wave height, subject to the amplitude of the forcing of the emulated SSL A + SSL = 2. The present strategy was initiated with the configuration given in Chernyshenko [18]. The wave height was chosen in order to approximately satisfy the abovementioned condition on the emulated SSL. Other configurations, a selection of which will be presented below, were later simulated, and the exploration was mainly undertaken by trial an error.

B. Computational simulations
Direct Numerical Simulations are performed using an in-house code, that 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 [19]. 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 [20, p. 510], the convergence of which is accelerated using a multigrid algorithm by Lien and Leschziner [21]. 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 halfheight, is met. A typical value used for this criterion is 10 −10 . Stable velocity-pressure coupling is ensured by use of the Rhie-and-Chow interpolation [22], preventing odd-even oscillations.
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 [23][24][25][26][27], 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 Maaß and Schumann [28], provided in the ERCOF-TAC 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 surfaceconformal 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 simula- tions correspond to the finest grid G6, at a bulk Reynolds number of Re b = 6200 (Re τ ≈ 360). Quantitative evidence for the level of refinement is given in figure 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 ∆ = 3 √ ∆x∆y∆z, 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 figure 4. In the former, shown in figure 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 figure 4 b, the local height varies from 2h − A w to 2h + A w , where 2h is the mean channel height, and A w the amplitude of the sinusoidal wall undulations.

Computational implementation for skewed flow
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 allows 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 zdirection 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. 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 D. An approximately constant flow rate across the channel is 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 data shows that the target streamwise bulk velocity is satisfied 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-TABLE I. Properties of the grid configurations, on which almost all the simulations presented are based. Viscous units are based on the target value of Reτ for the plane channel. (Note that the angle of the flow may vary, e.g. for an angle of π/2, the roles of x and z are reversed.) 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).

Grid label
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 phasevarying 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 plane-channel-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 pre- hin less than 0.001%. Note that, for θ = 0, x is distinct from the streamwise direction. roughout the discussion to follow, only the incremental part of the pressure is reported, tive to a pressure-reference value located at one of the corners of the computational . Quantities expressed in the frame of reference of the flow will be denoted with a rlaying inverted circumflex accent (·), as in fig. 2, although this notation will be omitted r when not needed. Unless stated otherwise, all physical interpretations are given in frame of reference of the flow, i.e. with x-and z-directions referred to respectively as amwise and spanwise directions.
ere q = [u, v, w, p] T is the solution vector, an overbar denoting the time average, y) = q x,z,t , q is the phase-varying part of the mean field, and q is the turbulent tuation. The second approach targets the action of the wavy boundaries onto the seuille flow: ere Q 0 is a turbulent Poiseuille profile, and q p is the perturbation to the plane-channel tion. Similar decompositions can be considered for any other statistical quantity, and l be referred to as 'phase-integrated' (Q), 'phase-varying' ( q), and 'perturbation' (q p ). All statistical quantities can be averaged in the direction parallel to the wave crests sented, the flow is driven at approximately constant flow rate, so that U b ≈ 1, via the imposition of a spatiallyconstant (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:F 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 [29] was of about 3-4%.

B. Overall physical characteristics
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 figure 6. 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 direction. The latter gradient generates a spanwise motion, shown in figure 7, 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 [12], this orientation was observed to vary in time, with strong reduction in turbulence during the reorientation phase. For the wavy wall, TABLE II. Configurations simulated, all at Reτ ≈ 360. Each simulation is designated by a label consisting of a set of letters plus identifiers, meant to convey as much information as possible in a compact manner, and identifying each calculation uniquely. 'G' stands for 'Grid' and refers to the grid configurations detailed in table I, 'W' stands for 'Wavy', 'P' for 'Plane', and 'A' for 'Amplitude'. The figure following the letter 'A' identifies a particular value of the wave slope Aw/λ. The suffix 'f' indicates the presence of a flat upper wall instead of two in-phase wavy walls, and 'bis' reflects the fact that the wavelength of 'G2W1bis' is equal to that of 'G2W1', but the flow angle θ is different. TDR, PD and FDR respectively stand for Total-Drag Reduction, Pressure Drag and Friction-Drag Reduction.

Simulation
Flow Drag coefficients Relative label configuration (×10 6 ) drag variation the phase-modulation of the shear strain occurs in space, and also results in a reorientation following the shearstrain field, as well as in a weakening of the streaks, as shown in figure 8 at a constant distance from the wall around y + ≈ 10. However, the effect is far less pronounced than observed in [12] because the forcing amplitude is much smaller in the present case.

C. Influence of the upper wall
Previous studies on wavy channels often had only a single wavy wall, or featured varicose wall undulations [30][31][32][33][34][35]. 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.

D. Net drag reduction and grid convergence
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 an the mesh. However, this is not so computationally, owing to variations in so-lution 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 example of the influence of the domain size on the drag, for a flow aligned with the mesh, Ricco and Wu [36] 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 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 wallparallel 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 [14,37], 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 (A + w , θ, λ + ) = (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. The main outcome of the study, shown in figure 9, is that the total-drag reduction predicted decreases as the mesh is refined, with a quadratic dependence on the wall-parallel mesh spacing. The drag-reduction value appears to tend, asymp- totically, 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, although at a slightly coarser resolution, in a computational box aligned with the main flow direction, and the wavy wall at an angle to the grid. Two wavelengths were included in the streamwise and spanwise directions. The resulting point, shown in figure 9, is in line with the simulations of the skewed configuration.

E. Friction-drag reduction and pressure-drag increase
In [17], 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 10 (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 [17], 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 FIG. 9. Total-drag reduction for A + w = 18, θ = 70 • , and λ + = 918, for grids G2, G3 and G6. Inverted triangle : supplementary run with the wave at an angle and the flow aligned with the grid. 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 [38], and assuming that the baseline and controlled drag levels are independent variables. 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 [12] 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.
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 nearwall 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 simulation presented in this section, and those of [17], 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 ≈ 1250, subject to the assumption that the Reynolds-number change from Re τ = 200 in [17], 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.
Results for some wavy channels are shown in figure 11. 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 phasewise 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 • ].

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 uncontrolled flow when actual scaling is used (i.e. with the modified friction velocity). This shift is shown in figure 12 (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 configuration λ + = 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 [12,17,39], 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.
The difference between the streamwise Reynolds-stress levels for the wavy walls and the baseline case is shown in figure 13 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 FIG. 11. Comparison of ∂w + /∂y + scaled by the streamwise-projected wave slope Aw/λ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 Aw/λ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.

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 14. 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 phasevarying 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 [18] 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 15(a), there exists a zone of intense production of turbulent kinetic energy Π k = − u i u k ∂u i /∂x k 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 15(b), the increase in production relative to the baseline is quickly followed, in phase, by an increase in energy dissipa-  tion k = −ν ∂u i /∂x k ∂u i /∂x k , at about the same wallnormal location. This phenomenon strengthens as the wave amplitude increases, and is stronger for G6W2* than for G6W1* at similar wave slopes.

D. Reynolds-number effect
An interesting question is whether the flow properties remain similar when the shape of the wall is kept constant in viscous units as the Reynolds number is increased. This has been investigated by reference to the flow configurations listed in table V.
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 figure 16. At Re τ = 180, a wave height of A + w = 20 represents a ratio A w /h of about 11%, which is significant, while it decreases to 2% for Re τ = 1000. Consequently, the wave height A w is 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 FIG. 12. Comparison of (a) the mean streamwise velocity U + = u x,z /uτ and (b) Reynolds stresses u i u j x,z /u 2 τ for wavy walls with wave slopes of 2Aw/λ ∈ {0.05, 0.07} (corresponding to labels *A3 and *A4) and wavelengths of λ + ∈ {612, 918} (corresponding to labels *W2* and *W1*), and the baseline plane channel (G6P1) (scaled using the actual friction velocity).
avoid an interference between the two solid boundaries, adding support to the findings reported in section III C.
Although net-drag-reduction levels of about 1-2% were observed for the three Reynolds numbers tested, the value is not quantifiable with any degree of precision, because it is extremely sensitive to various numerical issues, as demonstrated in section III D through a gridconvergence study at Re τ = 360.

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 wellestablished 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 cumula-tive 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 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 to 0.5% 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.