Turbulent drag reduction by spanwise wall forcing. Part 1. Large-eddy simulations

Abstract Turbulent drag reduction (DR) through streamwise travelling waves of the spanwise wall oscillation is investigated over a wide range of Reynolds numbers. Here, in Part 1, wall-resolved large-eddy simulations in a channel flow are conducted to examine how the frequency and wavenumber of the travelling wave influence the DR at friction Reynolds numbers $Re_\tau = 951$ and $4000$. The actuation parameter space is restricted to the inner-scaled actuation (ISA) pathway, where DR is achieved through direct attenuation of the near-wall scales. The level of turbulence attenuation, hence DR, is found to change with the near-wall Stokes layer protrusion height $\ell _{0.01}$. A range of frequencies is identified where the Stokes layer attenuates turbulence, lifting up the cycle of turbulence generation and thickening the viscous sublayer; in this range, the DR increases as $\ell _{0.01}$ increases up to $30$ viscous units. Outside this range, the strong Stokes shear strain enhances near-wall turbulence generation leading to a drop in DR with increasing $\ell _{0.01}$. We further find that, within our parameter and Reynolds number space, the ISA pathway has a power cost that always exceeds any DR savings. This motivates the study of the outer-scaled actuation pathway in Part 2, where DR is achieved through actuating the outer-scaled motions.


Introduction
Flow control aims to reduce drag on vehicles, enhance their efficiency, manoeuvrability, and possibly modify the heat transfer.Techniques for flow control cover a variety of fields, and have been extensively reviewed (White & Mungal 2008;Dean & Bhushan 2010;Luchini & Quadrio 2022).Flow control devices are usually divided into two groups: passive devices that are fixed in place and do not change their shape or function in time, such as vortex generators (Lin 2002;Koike et al. 2004;Aider et al. 2010) and riblets (García-Mayoral & Jiménez 2011, 2012;Endrikat et al. 2021a,b;Modesti et al. 2021;Endrikat et al. 2022;Rouhi et al. 2022), and active devices that can be actuated in some way, such as targeted blowing (Abbassi et al. 2017) or intermittent blowing and suction (Segawa et al. 2007;Hasegawa & Kasagi 2011;Yamamoto et al. 2013;Schatzman et al. 2014;Kametani et al. 2015).
Here, we are interested in a particular form of active control for drag reduction in wallbounded flows based on spanwise oscillation of the surface, leading to the generation of a streamwise travelling wave (Jung et al. 1992;Quadrio et al. 2009;Viotti et al. 2009;Quadrio 2011;Quadrio & Ricco 2011;Gatti & Quadrio 2012, 2013, 2016;Ricco et al. 2021).The wall motion is described by where w s is the instantaneous spanwise velocity of the wall surface, A is the amplitude of the spanwise forcing, ω is the angular frequency of oscillation, and κ x = 2π/λ is the wavenumber of the travelling wave with wavelength λ.Negative frequencies result in an upstream travelling wave, and vice versa.With an appropriate choice of A, κ x and ω, turbulent drag reduction beyond 40% can be achieved (Quadrio & Sibilla 2000;Quadrio et al. 2009;Hurst et al. 2014;Gatti & Quadrio 2016).The actuation mechanism (1.1) has been mostly investigated in a turbulent channel flow.So far, the only studies that investigate this mechanism in a turbulent boundary layer are the numerical work by Skote (2022), and the experimental work by Bird et al. (2018) and Chandran et al. (2022) in Part 2. The amount of drag reduction, DR, is defined as where C f ≡ 2τ w /(ρU 2 b,∞ ) and C f0 ≡ 2τ w0 /(ρU 2 b,∞ ) are the skin-friction coefficients of the drag-reduced flow (with wall shear-stress τ w ) and the non-actuated flow (with wall shear-stress τ w0 ) and ρ is the fluid density.The overbar in τ w and τ w0 indicates averaging over the homogeneous directions and time.In a fully-developed channel flow (considered here in Part 1), the averaging dimensions are the streamwise and spanwise directions, as well as time, and in a boundary layer (considered in Part 2), the averaging dimensions are the spanwise direction and time.Further, in a channel flow the drag-reduced flow and the non-actuated flow are exposed to the same bulk velocity U b (present Part 1, Quadrio et al. 2009;Gatti & Quadrio 2013) or pressure gradient (Quadrio & Ricco 2011;Ricco et al. 2012), however, in a boundary layer the two flows are exposed to the same free-stream velocity U ∞ (Part 2, Bird et al. 2018).Accordingly, there are two friction velocities u τ ≡ τ w /ρ and u τ0 ≡ τ w0 /ρ, corresponding to the drag-reduced and nonactuated cases, respectively, leading to two choices of normalisation.In the current study, following Gatti & Quadrio (2016), the viscous-scaled quantities that are normalised by u τ0 are denoted by the '+' superscript, and those normalised by u τ are denoted by the ' * ' superscript.The friction Reynolds number Re τ in a channel flow (Part 1) is defined based on u τ0 and the channel half-height h (Re τ ≡ u τ0 h/ν).In a boundary layer (Part 2), is defined based on u τ0 and the boundary layer thickness δ (Re τ ≡ u τ0 δ/ν).By dimensional analysis (Gatti & Quadrio 2016;Marusic et al. 2021) we obtain DR = DR κ + x , ω + , A + , Re τ , (1.3) where κ + x = κ x ν/u τ0 , ω + = ων/u 2 τ0 and A + = A/u τ0 .Quadrio et al. (2009) studied this flow control problem using direct numerical simulations (DNS) of a turbulent channel flow.Their study acted as a proof of concept for (1.1) to demonstrate that the introduction of a streamwise travelling wave achieves higher DR than a purely oscillating wall mechanism (κ x = 0).They fixed Re τ = 200 and A + = 12, and populated a map of DR(ω + , κ + x ) for 0 ≤ κ + x ≤ +0.04 and −0.3 ≤ ω + ≤ +0.3.Gatti & Quadrio (2016) extended this work to Re τ = 1000 and a broader range of actuation parameters (0 ≤ κ + x ≤ +0.05, −0.6 ≤ ω + ≤ +0.6 and 3 ≤ A + ≤ 15) to construct isosurfaces of DR(ω + , κ + x , A + ) in the 3D actuation parameter space (figure 4 in Gatti & Quadrio 2016).They observed that this type of actuation appears to modify the mean velocity profile through a Reynolds number-invariant additive constant, ∆B, in the logarithmic region as: where U * ≡ U/u τ and y * ≡ yu τ /ν are the viscous-scaled velocity and wall distance, κ and B are the von Kármán and additive constants for the non-actuated channel.This behaviour in U * implies that the actuation is primarily acting on turbulent structures in the near wall region and that the outer flow effectively perceives the modified inner layer as one that has a lower stress.This behaviour is similar to the flows over riblets and rough surfaces (Chan et al. 2015;Squire et al. 2016;Endrikat et al. 2021b) and Gatti & Quadrio (2016) used this assumption to propose the modified friction law (hereafter called GQ's model) given by (1.5) In this framework, the Reynolds number dependence of the flow is captured by C f0 , provided that there is a well-defined logarithmic region in the mean velocity profile.The behaviour of the log-region is modified by the actuation solely through the offset parameter ∆B; this parameter is independent of Reynolds number and can be parameterised by the dimensionless actuation parameters so that ∆B = ∆B(κ * x , ω * , A * ).The model therefore predicts DR at arbitrarily high Reynolds numbers for a given set of actuation parameters.The model also predicts that DR decreases monotonically with increasing Re τ , regardless of the actuation parameters.To date, the predictions from this model have been found to be largely consistent with the existing low-Reynolds number simulations of travelling wave drag reduction (Baron & Quadrio 1995;Yudhistira & Skote 2011;Ricco et al. 2012;Touber & Leschziner 2012;Hurst et al. 2014).
The findings reported so far are based on DNS of turbulent channel flow.Experiments have also reported the efficacy of spanwise wall forcing for turbulent drag reduction.The configurations are mainly turbulent boundary layer (Choi et al. 1998;Choi & Clayton 2001;Ricco & Wu 2004;Bird et al. 2018) or pipe flow (Choi & Graham 1998;Auteri et al. 2010).The experiments mostly consider uniform spanwise wall oscillation (i.e.κ x = 0 in 1.1).The exceptions are Auteri et al. (2010) and Bird et al. (2018) that attempt to mimic the travelling wave motion.Auteri et al. (2010) subdivide the pipe wall into thin slabs that rotate independently, and Bird et al. (2018) pneumatically deform a compliant structure.The experimental findings are consistent with the DNS findings.They report DR between 21% (Bird et al. 2018) to 45% (Choi & Clayton 2001).They also observe the shift in the log region (1.4) that underlies GQ's model (Choi et al. 1998;Choi & Clayton 2001;Ricco & Wu 2004).The DNS and experimental studies reviewed so far consider Re τ 1500.Marusic et al. (2021) recently investigated the parameter space (1.3) at much higher Reynolds numbers by conducting experiments up to Re τ = 12800 and wall-resolved largeeddy simulations (LES) up to Re τ = 2000.By covering such a large Reynolds number range, they were able to explore the increasing contribution of turbulent scales in the logregion and beyond to the total drag (Marusic et al. 2010;Smits et al. 2011;Mathis et al. 2013;Chandran et al. 2020).In contrast to previous studies, the drag reduction was found to occur via two distinct physical pathways.The first pathway, which Marusic et al. (2021) referred to as the "small-eddy" actuation strategy, as was applied in previous studies.It will be more aptly termed inner-scaled actuation (ISA) in the present work because drag reduction is achieved by actuating at frequencies associated with the nearwall cycle and the near-wall peak in turbulent kinetic energy.For example, ω + ≈ −0.06 equates to a time period of oscillation of T + osc = 2π/ |ω + | = 100.The DR obtained under this pathway was found to follow GQ's model.The second pathway, which Marusic et al. (2021) referred to as the "large-eddy" actuation strategy, was new.It involved actuating at frequencies comparable to those of the inertia-carrying eddies in the logarithmic region and beyond (T + osc ≫ 100).It will be more aptly termed outer-scaled actuation (OSA) in the present work.Unlike the ISA pathway, the OSA pathway achieves drag reduction that increases with Reynolds number, and requires significantly less input power due to the lower actuation frequencies that are required to target the inertia-carrying eddies.Marusic et al. (2021) considered actuation frequencies with T + osc < 350 to be primarily along the ISA pathway, and those with T + osc > 350 to be primarily along the OSA pathway.
In conjunction with Part 2 (Chandran et al. 2022), we investigate the drag reduction (1.3) over a range of parameters that have not been investigated previously, covering both the ISA and OSA pathways, and explain the physics behind the variation of DR with Re τ , κ + x and ω + .In this Part 1, we focus on the ISA pathway and use wall-resolved LES to extend the parametric study of Gatti & Quadrio (2016) at Re τ ≈ 1000, generating a new map of DR at Re τ = 4000 over 0.002 ≤ κ + x ≤ 0.02 and −0.2 ≤ ω + ≤ +0.2 for A + = 12.Accurately populating the DR map required a careful study of the LES setup in terms of the subgrid-scale model, grid and computational domain size, to ensure the accuracy of the simulations and computational tractability.The resulting map at Re τ = 4000 is used to evaluate the predictive accuracy of GQ's model, and by using turbulence statistics, triple decompositions, spectrograms and flow visualisations, we identify and explain the regimes of the flow at different regions of the DR map.We find that the flow regimes change with the extent of the Stokes layer generated by the surface motion.As the Stokes layer grows in size, up to the optimal range of 20 − 30 viscous units, the nearwall turbulence is damped, and there is a corresponding increase in DR.In contrast, growth beyond 30 viscous units amplifies the near-wall turbulence, leading to a decrease in DR.Finally, we examine the power cost at Re τ = 4000 over the range of parameters considered here.

Governing equations and solution method
We solve the filtered equations for a channel flow (figure 1) of an incompressible fluid with constant density ρ and kinematic viscosity ν The hat (...) indicates the filtered quantity; x 1 , x 2 and x 3 (also referred to as x, y and z) are the streamwise, wall-normal, and spanwise directions, corresponding to the velocity components u 1 , u 2 and u 3 (or u, v and w), respectively.The pressure gradient that is matched between the actuated and non-actuated cases.The unresolved subgridscale (SGS) stresses τ ij = u i u j − u i u j are modelled using the dynamic Smagorinsky model (Germano et al. 1991) incorporating Lilly's improvement (Lilly 1992).For the model coefficient, we perform xz-plane averaging of the inner products of the identity stresses (equation 11 in Lilly 1992).
Equations (2.1a,b) are solved using an LES extension of the DNS code by Chung et al. (2014).We perform wall-resolved LES in a channel flow (figure 1) by applying periodic boundary conditions in the streamwise and spanwise directions.At the bottom wall we apply u = v = 0 and w(x, z, t) = A sin(κ x x − ωt), and at the top boundary we apply freeslip and impermeable conditions (∂ u/∂y = ∂ w/∂y = v = 0).The present channel flow with free-slip top boundary conditions and domain height h (also known as open channel flow) requires less computational cost to converge compared to the conventional channel flow with no-slip top boundary conditions and domain height 2h (also known as full channel flow).Except for a small outer region, the mean velocity profiles and turbulence statstics are very similar between the two channel configurations (figures 2 and 3 in Yao et al. 2022).Compared to the boundary layer, we speculate small differences with (open) channel flow when we focus on the ISA pathway.This is supported by extensive comparison of channel flow with the boundary layer (Monty et al. 2009;Mathis et al. 2009;Chin et al. 2014).The two configurations have identical mean velocity profiles up to the end of the logarithmic region (see figure 1a in Monty et al. 2009).Up to the fourth-order statistics are in agreement between the two configurations to a height of half the boundary layer thickness (half the channel height), e.g.see figure 3 in Mathis et al. (2009).However, differences appear in the outer region due to the differences in the large-scale motions.Nevertheless, in the ISA pathway these large-scale motions do not contribute to DR.

Simulation cases
Table 1 lists all the simulations completed for Re τ = 951 and Re τ = 4000, where Re τ ≡ u τ0 h/ν represents the friction Reynolds number of the non-actuated case.At each Re τ , a parametric sweep of 7 × 8 combinations of streamwise wavenumber (κ + x ) and oscillation frequency (ω + ) is conducted over 0.00238 ≤ κ + x ≤ 0.02 and −0.2 ≤ ω + ≤ +0.2.The spanwise velocity amplitude is fixed at A + = 12.The seven non-zero values of ω + give the oscillation time periods T + osc = 126, 63, 42, and 31 (all within the ISA pathway), and ω + = 0 corresponds to a time-invariant standing wave in the streamwise direction.In Table 1, N/A denotes the specifications of the non-actuated simulation which serves as the reference case for calculating DR.For each actuated case, the DR is computed by matching the bulk Reynolds number Re b = U b h/ν between the actuated and non-actuated cases and substituting the respective values of C f and C f0 into (1.2).We consider matched bulk Reynolds numbers Re b = 19700 and 94450, which correspond to Re τ = 951 and 4000 for the non-actuated channel flow.Quadrio & Ricco (2011) and Ricco et al. (2012) compute C f and C f0 at matched Re τ (instead of matched Re b ) by driving the actuated and non-actuated cases with a constant pressure-gradient.Several differences exist between matching Re b (constant flowrate) and matching Re τ (constant pressure-gradient), see Quadrio & Ricco (2011), Quadrio (2011) and Ricco et al. (2012).With matched Re b , C f and C f0 are obtained at different Re τ .However, for our considered parameter space, the maximum DR is about 30%, which leads to a maximum deviation of about 16% in Re τ between C f and C f0 .Another source of difference between matched Re b and matched Re τ is in the actuation amplitude A (1.1).With constant A + = 12, A * = 12 for the actuated cases with matched Re τ .However, with matched Re b , A * > 12 when DR > 0, and vice versa.Nevertheless, DR weakly depends on A * 12 (Quadrio et al. 2009;Gatti & Quadrio 2016;Chandran et al. 2022).Overall, we speculate marginal differences in DR between matched Re b and matched Re τ for our parameter space.
The grid resolutions were chosen based on extensive validation studies as presented in Appendices A and B. In these appendices, we compare our LES results with DNS data of Gatti & Quadrio (2016) at Re τ ≈ 1000, experimental data of Marusic et al. (2021) at Re τ = 6000, and our self-generated DNS data at Re τ = 590.For DR and the mean velocity profile, we used the same viscous-scaled grid resolution at Re τ = 951 and 4000, corresponding to the streamwise and spanwise grid sizes of ∆ + x × ∆ + z ≃ 21 × 31 (the first seven rows at each Re τ in table 1).At this grid resolution, the difference in DR between the LES and DNS was found to be within 2%, and similarly good agreement was found for the mean velocity profile.However, for the Reynolds stresses and spectra at Re τ = 4000 we used a finer grid resolution with ∆ + x × ∆ + z ≃ 14 × 21 (the last two rows in table 1).Our nominal LES filter width ∆ + = ∆ + x ∆ + y ∆ + z 1/3 is 7 ∆ + 34 for the coarser grid, and 5 ∆ + 22 for the finer grid.However, given our anisotropic grid, we estimate our effective filter width from the two-dimensional energy spectrograms (figures 16e,f ).Our maximum filter width is in the spanwise direction and is about 50 and 35 viscous units for the coarser and finer grids, respectively, equivalent to the cut-off wavenumbers k + ∆z ≃ 0.12 and 0.18.These wavenumbers are 6 and 9 times larger than our maximum actuation wavenumber κ + x = 0.02.We estimate our cut-off frequency from Taylor's frozen turbulence hypothesis (Taylor 1938).The most challenging zone in terms of resolution is the buffer region (y + ≃ 10) with the smallest energetic eddies.If we take the convective speed of 10u τ0 in this region, our cut-off frequencies are ω + ∆z ≃ 1.2 and 1.8 for the coarser and finer grids, respectively, which are 6 and 9 times larger than our maximum actuation frequency ω + = ±0.2.
In terms of the domain size, the cases at Re τ = 951 used a full-domain with L x × The last row at each Reτ indicated with N/A for κ + x and ω + , corresponds to the non-actuated reference case.For all the actuated cases, A + = 12.Each row of the actuated cases consists of a set of cases with equal domain size, Reτ , κ + x and grid size, but ω + is different for each case (as listed in the fifth column).Those values of ω + with ± sign indicate two separate simulations, one with a positive sign (downstream travelling wave) and one with a negative sign (upstream travelling wave).The first column indicates the domain size (see figure 1); at Reτ = 951 we use the full domain and at Reτ = 4000 we use the medium domain.The third column y + res is the maximum resolved height by the simulation domain (≃ 0.4L + z , Chung et al. 2015).The eighth row at Reτ = 4000 repeats some of the cases with κ + x = 0.007 (the third row at Reτ = 4000), but with a finer grid resolution.
L z ≃ 6.6h × 3.2h (figure 1c), which is sufficiently large to resolve the first and secondorder statistics across the entire channel (Lozano-Durán & Jiménez 2014).However, at Re τ = 4000 each full-domain calculation is about 500 times more expensive than that at Re τ = 951, and so the domain size was reduced to L x × L z ≃ 2.0h × 0.6h (figure 1a).As a consequence, the flow is only resolved up to a fraction of the channel height y + res ≃ 0.4L + z (Chung et al. 2015), shown by the grey shaded zones in figure 1.For a reduced domain calculation, the user decides the resolved height y + res , with the constraint that it must fall somewhere in the logarithmic region.Then the domain size is obtained from the prescriptions of Chung et al. (2015) and MacDonald et al. (2017).For the travelling wave actuation (1.1), the prescriptions are L + z ≃ 2.5y + res , L + x max(3L + z , 1000, λ + ), where λ is the travelling wavelength.MacDonald et al. (2017MacDonald et al. ( , 2018) ) used the reduced-domain approach with 60 y + res 250 for turbulent flows over roughness.Endrikat et al. (2021b) used the same approach with y + res ≃ 100 for turbulent flows over riblets.Jiménez & Moin (1991) who used this approach for the first time resolved the flow up to y + res ≃ 80.They named this approach "minimal flow unit".Here, with L x × L z ≃ 2.0h × 0.6h (figure 1a) at Re τ = 4000, we resolve a substantial fraction of the inner layer up to y + res ≃ 1000.
Therefore, we name our reduced domain the "medium domain" to highlight its relatively larger size compared to the minimal flow unit.Gatti & Quadrio (2016) also used the medium domain size of L x × L z ≃ 1.4h × 0.7h with y + res ≃ 250 to study the travelling wave (1.1).In Appendix C, we assess the suitability of the medium domain size (figure 1a) by comparing the results with those obtained using a larger domain size (figure 1b) for selected cases from table 1.

Calculation of the skin-friction coefficient
To compute DR (1.2), we need the skin-friction coefficient for both the actuated and the non-actuated cases.Here, U * b = h * 0 U * dy * /h * is the viscous-scaled bulk velocity.For the cases at Re τ = 951 with the full domain size the U * profile is resolved across the whole channel and U * b can be found directly.However, for the cases at Re τ = 4000 with the medium domain size, the U * profile is resolved only up to y * res ≃ 750 − 1000.Two of these high Reynolds number profiles are shown in figure 2(a): the actuated case with A + = 12, κ + x = 0.02 and ω + = −0.05(blue lines), and the non-actuated case (black lines).The resolved portion of the LES profile below y * res is shown with a solid line, and the unresolved portion above y * res with a dashed-dotted line.We also overlay the DNS of the non-actuated full-domain channel flow at Re τ = 4200 by Lozano-Durán & Jiménez (2014) (red squares).For the non-actuated LES, the resolved portion up to y * res ≃ 1000 (solid black line) accurately reproduces the non-actuated DNS.However, the unresolved portion beyond y * res (black dashed-dotted line) departs from the non-actuated DNS due to the reduced domain size.
This issue has been addressed previously by Chung et al. (2015) MacDonald et al. (2017), Endrikat et al. (2021a), andEndrikat et al. (2021b).For accurate prediction of U * b , hence C f , it was found that the resolved height y * res must fall inside the logarithmic region, and it needs to be larger than the extent of the disturbed flow due to the surface modification.If y * res satisfies these criteria, the U * profile is resolved up to a portion of the log region, similar to the LES cases shown in figure 2(a).Beyond y * res , the unresolved portion of the log region and the outer region is assumed to be universal and so it can be reconstructed based on previous work.Here, we reconstruct the unresolved portions using the composite profile for full-domain channel flow (Nagib & Chauhan 2008).Figure 2(a) demonstrates that for the non-actuated case at Re τ = 4000, we obtain good agreement between the reconstructed profile for LES (dashed black line) and DNS.Therefore, to obtain U * b we integrate the resolved U * profile up to y * res and the reconstructed profile beyond y * res .We find that C f0 using this corrected U * b is only 1% different than the value obtained from DNS.
We follow the same approach to reconstruct the actuated U * profile (dashed blue line in figure 2a).However, we need to add the log-law shift ∆B in the composite profile to make the resolved and reconstructed profiles continuous at y * res .We find ∆B by plotting the velocity difference between the actuated and non-actuated profiles ∆U * = U * act − U * non-actuated (figure 2b).As seen in figure 2(b), ∆U * reaches almost a plateau beyond y * ≃ 100.We set ∆B as the value of ∆U * at y * res ≃ 750.Note that since the actuated u τ is smaller than the non-actuated u τo , y * res for the actuated case is about 750, but for the non-actuated case is about 1000.We calculate U * b for the actuated case by integrating the resolved portion of the profile up to y * res (solid blue line) and the reconstructed portion beyond y * res (dashed blue line).Another way of calculating U * b (hence C f ) from the reduced domain is to integrate the composite profile from y = 0 to h (e.g.see 4.2 in MacDonald et al. 2019), which assumes that the viscous sublayer and buffer layer make a negligible contribution to U * b .We believe that our present approach is more accurate as it considers the complex variation of U * in the viscous sublayer and buffer layer.We only use the composite profile in the log-region and beyond.

Drag reduction map as a function of frequency and wavelength
Figures 3(a,b) display the maps of DR(ω + , κ + x ) at Re τ = 951 and 4000 from the computations listed in table 1.At each Re τ , we have 56 DR data points.To generate the maps, we perform bilinear interpolation of our DR data points onto a uniform 20 × 20 grid over the parameter space 0 ≤ κ + x ≤ 0.02 and −0.2 ≤ ω + ≤ +0.2.At Re τ = 951, the maximum DR of 35.4% at (ω + , κ + x ) = (0.05, 0.021) is in close agreement with the DNS of Gatti & Quadrio (2016) at Re τ ≃ 950, where the maximum DR was found to be 38.8% at (ω + , κ + x ) = (0.05, 0.0195).At Re τ = 4000, the maximum DR decreases to 27.5% at the same actuation parameters (ω + , κ + x ) = (0.05, 0.021).At each Reynolds number, DR changes more drastically by changing ω + than by changing κ + x .When κ + x = 0, there is no travelling wave (plane wall oscillation), and the variation of DR is symmetric between ω + < 0 and ω + > 0. In this case, two equal local maxima (at ω + ≃ ±0.05) and a local minimum (at ω = 0) emerge.When κ + x > 0, a travelling wave is generated, and the variation of DR is asymmetric between ω + < 0 and ω + > 0. In this case, at each κ + x only one local maximum (blue dashed-dotted curve in figure 3) and one local minimum (black dashed curve in figure 3) appear in DR.These observations are in agreement with Quadrio et al. (2009) and Gatti & Quadrio (2016).Overall, within our parameter space, the map of DR consists of three distinct regions.Region I to the left of local maximum DR (blue dashed-dotted curve) where ω + 0 (upstream travelling wave); in this region DR > 0, hence drag reduction.Region II represents the crossover from the local maximum to the local minimum DR (between the blue dashed-dotted curve and the black dashed curve).For κ + x 0.007, the local minimum DR is positive, however, for κ + x 0.007 the local minimum DR becomes negative (hence a drag increase).

Increase in κ +
x beyond 0.007 leads to a larger drag increase area, and the local minimum DR becomes more negative; Quadrio et al. (2009) and Gatti & Quadrio (2016) observe similar trends.Quadrio et al. (2009) find that the local minimum DR follows the line ω + /κ + x ≃ 10.In other words, maximum drag increase occurs when the travelling wave speed is about 10u τ0 , which is nearly the same as the convective speed of the near-wall flow structures.Similarly, in figures 3(a,b) the black dashed curve that marks the local minimum DR follows ω + /κ + x ≃ 10.Region III covers the right of local minimum DR (black dashed curve) where ω + > 0 (downstream travelling wave); in this region, increase in ω + increases DR.
In figure 3(c), we display the difference in DR as the Reynolds number changes from 4000 to 951.For most of the (ω + , κ + x ) space, DR is lower at the higher Reynolds number.Only within the range 0.005 κ + x 0.020, +0.1 ω + +0.2 we observe the opposite trend.This region coincides with the drag increasing range (DR < 0) with ω + > 0. This observation is consistent with GQ's model (1.5),where DR < 0 (hence ∆B < 0) predicts an increase in DR as Reynolds number increases.To make these comparisons more quantitative, in figure 3(d ) we show the difference in DR between our results and GQ's model at Re τ = 4000.To predict DR, the model (1.5) requires C f0 and the value of the log-law shift ∆B for each set of actuation parameters (A + , κ + x , ω + ).For C f0 , we use Dean's power-law correlation (Dean 1978) which agrees well with the DNS data given by MacDonald et al. (2019).We can obtain ∆B from a low Reynolds number simulation for the same set of (A + , κ + x , ω + ) because ∆B is assumed to be Reynolds number independent.Therefore, we use our results at Re τ = 951, where for each (ω + , κ + x ) we find ∆B from the velocity difference ∆U * at y * = 200 (similar to figure 2b).We choose y * = 200 as it is far enough from the wall to fall into the log region, but not too far to fall into the wake region (y/h 0.3 according to Pope 2000).By having ∆B at each (ω + , κ + x ) and having C f0 at Re τ = 4000, we can reconstruct the DR map based on GQ's model.
Figure 3(d ) shows the overall good performance of GQ's model for this range of Reynolds numbers.In region I, the difference in DR between LES and GQ's model is less than 2%, i.e. |DR% LES,Reτ =4000 − DR% GQ,Reτ =4000 | 2%.This is a very good agreement considering that DR varies between 15% and 28% in region I.In regions II and III, we observe some slight differences in DR between LES and GQ's model, especially in region II in the drag increasing range.In this range, the difference in DR between LES and GQ's model reaches 4%, which is the same order as DR (see figure 3b).In region III, for ω + +0.1 again we observe good agreement between LES and GQ's model (less than 2% difference).In the following sections, we investigate the reasons behind the different performance of GQ's model in regions I, II and III related to the changes in the Stokes layer dynamics and the near-wall turbulence in each of these regions.

Mean velocity profiles
To obtain an overall picture of the mean velocity behaviour in regions I, II and III (see figure 4), we consider the seven runs conducted at Re τ = 4000, A + = 12 and κ + x = 0.007 for ω + ranging from −0.2 to +0.2.In figure 4(a), we identify the selected values of ω + (filled squares) on the DR map along with the local maximum DR ( ) and the local minimum DR ( ).The corresponding velocity profiles are shown in figures 4(b,c) for ω + ≤ 0 (upstream travelling waves) up to the local maximum DR (region I) and ω + ≥ 0 (downstream travelling waves) beyond the local maximum DR (regions II, III), respectively.
When ω + ≤ 0, the log region of the actuated profiles is shortened and shifted above the non-actuated counterpart (figure 4b), corresponding to a positive DR.The shortening of the log region is due to the thickening of the viscous sublayer.We show the viscous sublayer thickening in the inset of figure 4(b), in that the actuated profiles of U * /y * are closer to unity for a greater wall distance compared to their non-actuated counterpart.We show the shortening of the log region in figure 4(d ) by plotting the diagnostic function y * dU * /dy * .The log region appears as a plateau with the value of κ −1 ≃ 2.5.For the non-actuated case, the plateau appears for 100 y * 600.This range is consistent with the DNS of channel flow by Lozano-Durán & Jiménez (2014) and Lee & Moser (2015) at Re τ = 4200 and 5200, respectively (see figure 3a in Lee & Moser 2015).For the actuated cases, the plateau is narrowed further (i.e.log region is shortened) as DR increases.We quantify the shift in the log region by plotting ∆U * (inset of figure 4d ).The magnitude of the shift increases as DR increases.These observations are also reported in the previous turbulent drag reduction studies, including turbulent flow with the spanwise wall oscillation (Di Cicca et al. 2002;Touber & Leschziner 2012;Hurst et al. 2014), turbulent flow with the streamwise travelling wave (Hurst et al. 2014;Gatti & Quadrio 2016), turbulent flow of a polymer solution (Ptasinski et al. 2003;White & Mungal 2008), and turbulent flow over piezoelectrically excited travelling waves (Musgrave & Tarazaga 2019).Gatti & Quadrio (2016) derived their predictive model (1.5) based on similar observations of the velocity profiles in region I, and as a result GQ's prediction works well in this region (figure 3d ).The behaviour of the profiles in region I is consistent with the ISA pathway, where only the inner-scale eddies up to the buffer region are actuated.
In region II, we observe a sudden drop in DR as ω + changes from 0 to +0.1 (figure 4a), with a corresponding decrease in the logarithmic shift (figure 4c).A distinct feature of region II is the high level of distortion in the U * profile, which is particularly severe at ω + = +0.05.For this case, the diagnostic function tends towards the plateau κ −1 , but does not reach it.Similarly, ∆U * for this case approaches a plateau of 1.7 by the resolved height y * res ≃ 750, but does not reach it (inset in figure 4e).This is our most challenging case for computing DR using our approach in § 2.3 (figure 2).For accurate calculation of DR, ∆U * needs to reach a plateau by the resolved height y * res ≃ 750, i.e. the resolved height must fall into the logarithmic region.In Appendix C, we deliberately consider this challenging case for domain size study.We double the domain length and width compared to the medium domain (figure 1b), extending the resolved height to y * res ≃ 1500.The difference in DR is 1.4% between the medium domain and the large domain (table 4).Further, the large domain reinforces the approach of ∆U * to a plateau of 1.7 (the inset in figure 18b).To our knowledge, such significant levels of distortion in the U * profile have not been seen before in previous studies of flows over drag-reducing or drag-increasing surfaces.For example, in rough wall turbulent flows ∆U * is almost constant for y * 30 (e.g.figure 6  100 (e.g.figure 2 in Endrikat et al. 2021b).In § 3.4 and 3.5, we discuss the physics behind the highly distorted mean velocity profiles (figures 4c,e).
In region III, when ω + increases to +0.2 (figure 4c), DR increases to 13% and the U * profile behaves similarly to that seen in region I.A well-defined logarithmic shift appears beyond y * ≃ 100 with viscous sublayer thickening.

Turbulence statistics
We now assess the behaviour of the Reynolds stress distributions at Re τ = 4000 (figures 5a-d ) and the turbulent kinetic energy production P = − u ′ v ′ xzt dU/dy (figures 5e,f ), where ... xzt denotes averaging over xz-plane and time.We highlight four cases from figure 4 (A + = 12, κ + x = 0.007), where we vary ω + from −0.05 to +0.20.As indicated earlier, we employ a finer grid resolution for these cases to properly resolve the Reynolds stresses (see Appendix B).We plot the profiles scaled by the non-actuated u τ0 (dashed-dotted lines, figures 5a,c) and by the actuated u τ (solid lines, figures 5b,d ).Scaling by u τ0 is comparable to scaling by the bulk velocity U b (Gatti & Quadrio 2016) because the bulk velocity U b is the same between the actuated and non-actuated cases.Any difference between the outer-scaled actuated and non-actuated profiles reflects the overall response of turbulence to the wall oscillation (1.1).
Scaling by the non-actuated u τ0 ('+' superscript), as in figures 5(a,c,e), indicates that the wall oscillation attenuates the u ′2 + xzt levels up to the resolved height y + res ≃ 1000.The cases with the highest DR (ω + = −0.05,0 in figure 5a) show the highest level of attenuation in their u ′2 + xzt .Additionally, for these cases the inner peak of u ′2 + xzt is farther from the wall.Consistently, the viscous sublayer is thickened and the buffer layer is shifted away from the wall (figure 4b).In contrast to the behaviour of u ′2 + xzt , the w ′2 + xzt profiles are amplified near the wall.According to Quadrio & Ricco (2011) and Touber & Leschziner (2012), this amplification is due to the Stokes layer that forms as a result of the spanwise wall motion.The pre-multiplied turbulent kinetic energy production y + P + (figure 5e) also displays the attenuation of turbulence that accompanies increasing DR.All these trends are similar to previous studies on spanwise wall oscillation at lower Reynolds numbers (Quadrio & Ricco 2011;Touber & Leschziner 2012).
Scaling by the actuated u τ (' * ' superscript) is equivalent to inner scaling, which highlights the extent up to which the actuated profiles depart from the non-actuated profile.For u ′2 * xzt and w ′2 * xzt (figures 5b,d ), the actuated cases agree with the nonactuated case at distances far from the wall, but near the wall the actuated u ′2 * xzt levels are attenuated, while the w ′2 * xzt levels are amplified.For ω + = +0.05(in region II), the point where the actuated profiles begin to depart from the non-actuated counterpart occurs at y * ≃ 100, considerably farther than for the other cases (y * 30).The same case yields the strongest level of near-wall amplification for w ′2 * xzt (the red profile in figures 5c,d ) and the highest level of distortion in mean velocity (red profile in figures 4c,e).
Regardless of the scaling used, as w ′2 xzt is amplified near the wall, u ′2 xzt is attenuated, the viscous sublayer is thickened and DR is increased.This trend occurs in regions I (ω + = −0.05,0) and III (ω + = +0.2).In region II (ω + = +0.05),however, there is an excessive amplification of w ′2 xzt near the wall, a thinning of the viscous sublayer and a drop in DR.

Stokes layer: an important source of inner-scaled actuation
As indicated earlier, the near-wall amplification of w ′2 xzt is related to the growth of the Stokes layer.We now apply triple decomposition to more precisely uncover how the strength of the Stokes layer modifies the near-wall turbulence, which in turn affects the wall drag.We primarily consider u τ scaling, as we are interested in the level of departure from the non-actuated behaviour.In Part 2, we mostly use u τ0 scaling, as we are interested to study the overall response of turbulence to the wall actuation.Nevertheless, the conclusions from Parts 1 and 2 are valid regardless of the scaling.
Because the flow is subjected to a harmonic forcing (1.1), the instantaneous flow can be triply decomposed similar to Touber & Leschziner (2012), as in where f indicates the quantity of interest, i.e. u, v or w.In (3.1a), the total fluctuations f ′ is decomposed into the harmonic contribution f and the stochastic (turbulent) contribution f ′′ .The harmonic contribution f is obtained by phase averaging the spanwise averaged field f z in time over the number of periods N , and then subtracting the mean vertical profile f xzt .Accordingly, the total Reynolds stress f ′2 xzt is decomposed into its harmonic component f 2 xt associated with the Stokes layer dynamics and its turbulent (stochastic) component f ′′2 xzt (3.1c).In figure 6, we plot these two components for the cases given in figures 4 and 5 (A + = 12, κ + x = 0.007, Re τ = 4000).For reference, figures 6(a,b) show the considered U * profiles (as in figures 4b,c).Figures 6(c,d 6f ) the decay rate in w2 * xt is noticeably slower compared to regions I and III, implying the presence of a more protrusive Stokes layer.Accordingly, the U * profile in region II is distorted to the highest level.The turbulent stress profiles are also shown in figure 8, where they are accompanied by the turbulent kinetic energy profiles K xzt , which follow the same trends.
To quantify the protrusion of the Stokes layer (figures 6e,f ), we calculate two length scales from the spanwise Reynolds stress profiles.The first is the laminar Stokes layer thickness δ * S that is featured in Stokes' second problem (Batchelor 2000).Following Quadrio & Ricco (2011), we define δ * S as the height y * where the amplitude of w decays to Ae −1 (i.e.where w2 * xt = 1 2 A * 2 e −2 ).In figure 6, we mark δ * S on each profile with a cross symbol.The second length scale ℓ * 0.01 is new, and it is defined as the height where w2 * xt = 0.01.Our choice for the threshold of w2 xt is a small fraction of its maximum value at the wall A * 2 /2.Thus, we ignore the background turbulence in this definition.However, we mark ℓ * 0.01 where the Stokes layer stress w2 * xt is a small fraction of the turbulent stress w ′′2 * xzt , hence considering the background turbulence in this definition.In figures 6(c-f ), ℓ * 0.01 coincides well with the distance where the actuated u ′′2 * xzt and w ′′2 * xzt profiles depart from the nonactuated counterpart.However, δ * S underestimates the actual protrusion by the Stokes layer due to its ignorance of the background turbulence.For instance, for the case with ω + = 0 at y * = δ * S (black cross symbol in figure 6e), w2 * xt ≃ 8 w ′′2 * xzt , i.e. the Stokes layer is 8 times stronger than the background turbulence.However, at y * = ℓ * 0.01 (black bullet symbol) w2 * xt ≃ 0.01 w ′′2 * xzt , i.e. the Stokes layer is 100 times weaker than the background turbulence.We propose, therefore, that ℓ * 0.01 is a more suitable measure for reflecting the entire penetration of the Stokes layer into the turbulent field.
In regions I and III, the level of protrusion by the Stokes layer ℓ * 0.01 , as well as the departure height in the u ′′2 * xzt and w ′′2 * xzt profiles, stay below 20−30 viscous units.As a result, the mean velocity profiles in regions I and III (figures 6a,b) yield a well-defined logarithmic shift beyond y * ≃ 100 with viscous sublayer thickening.However, in region II there is a large increase in ℓ * 0.01 and the departure in the u ′′2 * xzt and w ′′2 * xzt profiles also starts at larger distance from the wall.For example, for ω + = +0.05 in region II ( figures 6d,f ), ℓ * 0.01 ≃ 80, which also closely marks the point where the actuated u ′′2 * xzt and w ′′2 * xzt profiles depart from their non-actuated counterpart.As a result, the mean velocity profile for ω + = +0.05 in region II (figure 6b) is highly distorted up to y * ≃ 200 − 300.
Furthermore, we can draw a connection between the protrusion height ℓ * 0.01 and the level of drag reduction.In figures 7(a,b), we overlay the map of ℓ * 0.01 onto the maps of DR and δ * S .In region I (left side of the blue dashed-dotted line), ℓ * 0.01 30 and δ * S 7. In this region, an increase in ℓ * 0.01 and δ * S leads to an increase in DR.For upstream travelling waves (ω + < 0), therefore, the growing protrusion of the Stokes layer has a favourable effect on DR.In contrast, in region II (between the blue dashed-dotted line and the black dashed line), DR drops by increasing ℓ * profile by phase averaging the simulation data.The variation in DR versus ℓ * 0.01 is similar to DR versus δ * S up to region I, in terms of the linear trends, an optimal thickness for the maximum DR and a minimum thickness for the occurrence of DR (δ * S,min , ℓ * 0.01,min ).However, in region II when the linear trends are broken, we observe noticeable differences between DR versus δ * S and DR versus ℓ * 0.01 .In region II, there does not appear to be a consistent relation between DR and δ * S (the red symbols in figure 7d ).In other words, we cannot find a threshold for δ * S beyond which DR drops.For instance, for the case with κ + x = 0.02, ω + = +0.2,DR drops to −10% but δ * S ≃ 5 which is within the optimal range (5 δ * S 7).In contrast, there is a much stronger connection between DR and ℓ * 0.01 , even in region II (figure 7c).For all cases, increasing ℓ * 0.01 beyond 30 decreases DR.As a result, the value of ℓ * 0.01 can be used to determine whether we are in region I (ℓ * 0.01 30) or region II (ℓ * 0.01 30).

Interaction between the Stokes layer and the near-wall turbulence
In a turbulent flow with spanwise wall oscillation, Touber & Leschziner (2012) similarly report that an overly protrusive Stokes layer leads to the degradation of DR.They proposed that the attenuation of u ′′2 xzt and amplification of w ′′2 xzt are based on the periodic realignment of the near-wall streaks.To examine this proposal further, we consider energy spectrograms and near-wall flow visualisations (figures 8 and 9).We focus on the same cases as in figure 6, where Re τ = 4000, A + = 12, and κ + x = 0.007.For the cases with ℓ *  layer).For these cases, increasing ℓ * 0.01 (hence strengthening the Stokes layer) attenuates u ′′ over a wider range of wavelength λ * z and height y * (e.g.compare figure 8k with 8l ).At the same time, the energetic peak in k * z φ * u ′′ u ′′ is shifted to a higher y * and a higher λ * z .This attenuation is apparent in the visualisations of the instantaneous velocity fields of u ′′ at y * = 10 for the cases with ℓ * 0.01 30 (figures 9e,f ).On the w ′′ fields, we overlay the spanwise and phase-averaged w * (solid black curves) as a measure of the Stokes motion.As ℓ * 0.01 increases from 20 at ω + = −0.05(figure 9e) to 30 at ω + = 0 (figure 9f ), the Stokes motion becomes stronger and the energy level in the u ′′ field is decreased compared to the non-actuated counterpart (figure 9i).At the same time, the spanwise spacing between the high-speed streaks increases by increasing ℓ * 0.01 .Overall, for ℓ * 0.01 30, the Stokes layer dampens the level of turbulence within y * ℓ * 0.01 , hence acting favourably towards increasing DR.
For the cases with ℓ * 0.01 > 30, the Stokes layer is excessively strong and protrusive.As a result, the near-wall flow structures meander, following the Stokes motion.This meandering is observed in the u ′′ and w ′′ fields for the cases with ℓ * 0.01 > 30 at y * = 10 (figures 9g,h).Even at y * = 50 for the same cases (figures 9l,m), we can see the protrusion of the Stokes motion (solid black curves) and evidence of meandering in the u ′′ fields.This meandering is also evident in the spectrograms.For instance, for ω + = +0.05 with ℓ * 0.01 ≃ 90 (figures 8m,r ), the meandering flow structures at y * ≃ 10 (visualised in figure 9g) manifest as an energetic peak in the k * z φ * w ′′ w ′′ spectrogram (figure 8r ) at (λ * z , y * ) ≃ (100, 10); this peak coincides with the peak in the k * z φ * u ′′ u ′′ spectrogram (figure 8m).Touber & Leschziner (2012) relate the attenuation of u ′′2 xzt and amplification w ′′2 xzt (e.g.figures 9c,d ) to this meandering behaviour and argue that the strong Stokes shear strain periodically re-orients the streaks.As a result, energy is transferred from u ′′ to w ′′ , and the anisotropy between u ′′ and w ′′ is reduced.Considering the flow visualisation at y * = 10 for ω + = +0.05 with ℓ * S ≃ 90 (figure 9g), we see a strong resemblance between the u ′′ and w ′′ fields in terms of the energy level and structure, which support the reduction in anisotropy.
A noticeable difference between the cases with ℓ * 0.01 30 and those with ℓ * 0.01 > 30 is the wall distance of the maximum turbulence activity given by the location of the energetic peaks in k  In each of (e-n), the upper field shows the turbulent streamwise velocity u ′′ * and the lower field shows the turbulent spanwise velocity w ′′ * .We overlay the spanwise and phase-averaged spanwise velocity w * (as solid curves) onto the w ′′ * field.

Power performance analysis
While drag reduction is an important performance parameter for many applications, the efficiency of the flow control effort is often even more important.Here we use the concept of net power saving N P S: where P + = (1 − DR) U + b is the pumping power required to drive the flow through the actuated channel, P + 0 is the non-actuated analogue of P + , and P + in is the input power required to oscillate the wall actuation mechanism (1.1) while neglecting any mechanical losses.A positive N P S indicates that the total power cost of the actuated case is less than the total cost of its non-actuated counterpart.We are also interested in assessing the accuracy of generalised Stokes layer (GSL) theory (Quadrio & Ricco 2011) for estimating P + in .In Part 2 (Chandran et al. 2022), we use this theory to estimate N P S for our experimental data.Here, in Part I, our actuation frequencies fall into the ISA regime.In Part 2, the data fall into both the ISA and OSA regimes.
The input power is given as follows, as first proposed by Baron & Quadrio (1995) for an oscillating plane, and then used by Quadrio & Ricco (2011), Gatti & Quadrio (2013) and Marusic et al. (2021) for a travelling wave.
where all the quantities are normalised by ν and the non-actuated u τo (hence superscripted with a cross symbol).In (3.2), T avg is the averaging time, w s is the instantaneous wall velocity (1.1) and ∂w + /∂y + | y + =0 is the instantaneous wall-normal gradient of the spanwise velocity at the wall.
In figure 10(a), we present the map of P + in /P + 0 as computed over our parameter space of (ω + , κ + x ) at Re τ = 4000.The map is much more symmetric about ω + = 0 compared to DR (figure 3).We also see that substantially more power is required at higher actuation frequencies.For example, P + in /P + 0 % can reach up to 100% when ω + ≃ ±0.2.In region II, between the local maximum and the local minimum DR (between the blue dashed-dotted line and the black dashed line), P + in /P + 0 decreases to about 30% − 35%.We can use (3.2) only if we have an estimate for ∂w + /∂y + | y + =0 .In most experimental studies, including Part 2 of the present study, this quantity is unavailable and some estimate needs to be made instead.In Part 2 we use GSL theory, which gives the instantaneous spanwise velocity for a laminar flow with wall actuation.That is, , Ai is the Airy function of the first kind, and R{...} is the real part of the argument.To use (3.3) for a turbulent flow, one needs to assume that 1) the Stokes layer preserves its laminar structure near the wall, i.e., w is the same in the laminar and turbulent flow, and 2) the turbulent spanwise velocity is negligible near the wall, i.e., w ′′ ≃ 0. We now compare the results for P + in using GSL to the results obtained using LES, so as to verify the validity of using GSL estimates in experiments.Figure 10(b) shows the difference between the pumping power obtained from the LES at Re τ = 4000 (P + in /P + 0 % LES,Reτ =4000 ) and that estimated using GSL theory at the same Re τ (P + in /P + 0 % GSL,Reτ =4000 ).Overall, the differences are small, especially in the range ω + < 0 where they differ less than 3%.Only in region II with ω + > 0 (between the blue dashed-dotted line and the black dashed line) do the differences approach 10%, especially along the minimum drag reduction line (black dashed line).The overlay of the Stokes layer protrusion height ℓ * 0.01 (contour lines) indicates that the region where the power differences are significant coincides with the region where ℓ * 0.01 is large.In other words, the error in using GSL theory (laminar Stokes layer assumption) is largest when the Stokes layer is most protrusive.
Our observations regarding the differences between P in from the simulation and that from the GSL theory (figure 10) are similar to those reported by Quadrio & Ricco (2011) (their figure 7); they report close agreement between the GSL theory and the turbulence simulation in the drag-decreasing range, but report noticeable differences in the dragincreasing range.They explain this behaviour through the timescale T + ≡ 2π/(ω + − κ + x U + w ), which represents the period of oscillation as observed by the near-wall eddies with the convective speed U + w ≃ 10.As discussed in § 3.1, in the drag-increasing range ω + /κ + x → U + w leading to T + → ∞.In other words, the spanwise oscillation becomes too slow that close to the wall, the u-and w-momentum equations are coupled together.However, GSL theory assumes that these equations are decoupled.Here, we add a new explanation based on the protrusion of the Stokes layer.As discussed in § 3.5 (figure 9), in the drag-increasing range, the Stokes layer is too protrusive and a near-wall cycle of turbulence is embedded within the Stokes layer.As a result, near the wall, all the terms of the momentum equation (2.1b) are active.However, the GSL theory neglects the advection (non-linear) terms from the w-momentum equation.The departure of the w2 * xt profiles from the GSL solution (figure 11b) supports the activation of these terms.We can now explore the net power saving (N P S). Figure 12(a) demonstrates that for our considered parameter space N P S is mostly negative.The highest (best) N P S is 0.5% at (ω + , κ + x ) ≃ (0.05, 0.012).In figure 12(b), we plot the map of the difference between N P S from LES at Re τ = 4000 and its counterpart at Re τ = 951.If this difference is positive, N P S increases with Re τ .Over a large portion of our parameter space, the difference is negative, i.e., N P S becomes more negative with Re τ .However, for a small portion of region I with ω + < 0 and κ + x 0.0025, the difference is positive.One experimental case reported by Marusic et al. (2021) falls into this region, with ω + = −0.044and κ + x = 0.0014 (see their figures 3a,b).The N P S of this case was negative, but it increased with Re τ in accordance with our analysis.Quadrio et al. (2009), similar to figure 12(a), generate a map of N P S for their travelling wave study at Re τ = 200 (their figure 5).They report N P S > 0 within the range 0 ω + +0.05, 0.002 κ + x 0.025.This range coincides with the range where DR 40%.Gatti & Quadrio (2013) generate a similar map at Re τ = 1000 (their figure 9).They also observe N P S > 0 within the same range of (ω + , κ + x ).However, the level of N P S > 0 is lower at Re τ = 1000 compared to Re τ = 200.Considering (3.4), we speculate that the decrease in N P S > 0 from Re τ = 200 to 1000 is due to the decrease in DR.We observe a similar trend in figure 12(b).Within the range of (ω + , κ + x ) where DR is maximum (blue dashed-dotted line), N P S decreases by increasing Re τ from 1000 to 4000.
Overall, for our considered parameter space, N P S is negative and predominantly decreases with Reynolds number.As discussed in § 1, our parameter space falls into the inner-scaled actuation (ISA) regime.In Part 2, we conduct experiments with some actuation parameters in the outer-scaled actuation (OSA) regime, which yield positive values for the N P S that actually increase with Reynolds number.

Conclusions
Turbulent drag reduction was considered using spanwise wall oscillation based on streamwise travelling waves at friction Reynolds numbers Re τ = 951 and 4000 using wall-resolved large-eddy simulation in a channel flow.We conducted parametric studies at both Reynolds numbers with a fixed actuation amplitude A + = 12, for wavenumbers and frequencies within the range 0.002 ≤ κ + x ≤ 0.02 and −0.2 ≤ ω + ≤ +0.2, covering upstream (ω + < 0) and downstream (ω + > 0) travelling waves.Our actuation parameters fall into the inner-scaled actuation (ISA) regime, where only the near-wall scales are actuated.
We find that GQ's model for the variation of drag reduction with Reynolds number performs well if the logarithmic shift in the velocity profile is accurately calculated.The present travelling wave actuation can highly distort the mean velocity profile and extend the beginning of the logarithmic region beyond 200 viscous units above the surface.We find that such a high level of distortion is related to the protrusive Stokes layer.Accordingly, we propose a length scale ℓ 0.01 for the protrusion height, where the Reynolds stress due to the Stokes layer drops to 1% of the Reynolds stress due to the background turbulence.We find that depending on ℓ 0.01 , hence the Stokes layer protrusion, the DR map over the parameter space of (ω + , κ + x ) can be categorised into two regions.When ℓ 0.01 is less than 30 viscous units, increasing ℓ 0.01 leads to an increase in DR.In this regime, the viscous sublayer is thickened and the logarithmic region appears at a point about 100 viscous units above the wall.The Stokes layer acts to attenuate the turbulence below ℓ 0.01 and lifts the cycle of turbulence generation away from the wall.Increasing ℓ 0.01 in this regime, further attenuates the turbulence and leads to higher DR.When ℓ 0.01 exceeds 30 viscous units, however, increasing the Stokes layer thickness leads to a drop in DR.In this regime, the logarithmic region appears beyond 200 viscous units above the wall.The decrease of DR in this regime is due to the Stokes layer becoming strong enough to cause a meandering of the near-wall turbulence, rather than attenuating it.That is, a cycle of near-wall streaks appear within 10 viscous units that follow the Stokes oscillatory motion.
Our power cost analysis showed that generalised Stokes layer theory agrees reasonably well with the LES data, so that it can be used with some confidence in cases where the gradient of the velocity at the wall is not accessible, as in most experiments.In addition, for our considered range of ω + and κ + x at Re τ = 4000 the net-power savings (N P S) was always negative.In other words, the power cost necessary to oscillate the near-wall fluid exceeds the power savings by the drag reduction.We speculate that negative N P S is inevitable in the ISA pathway at least at high Reynolds numbers.We confirm this speculation in Part 2, where we investigate the ISA and OSA pathways experimentally at Re τ up to O(10 4 ).
We conclude that with the viscous-scaled grid resolution of ∆ + x × ∆ + z ≃ 23 × 31 (blue diamond) we can study DR with high confidence.Therefore, we adopt this grid resolution to study DR (table 1).Table 3. Simulation cases for assessing the LES grid for studying the Reynolds stresses and their spectra (Appendix B).All the cases have the same Reτ , actuation parameters (A + , κ + x , ω + ) and domain size Lx, Lz, where h is the open-channel height.The top two cases are LES with coarse and fine grid resolutions, respectively.The third case is DNS.We perform two LES calculations that match the DNS case in terms of the domain size, Re τ and actuation parameters, but have different grid resolutions (table 3).We name the LES case with a coarser grid (∆ + x × ∆ + z ≃ 22 × 31) "coarse LES", and the case with a finer grid (∆ + x × ∆ + z ≃ 15 × 21) "fine LES".Note that the coarse LES case still has a fine grid for wall-resolved LES.Previous LES studies have employed a similar grid size to study a turbulent wall jet (Banyassady & Piomelli 2014) or separating turbulent boundary layer (Wu & Piomelli 2018).Further, the coarse LES grid predicts DR quite well (Appendix A).
In figures 15 to 17, we compare coarse and fine LES cases with DNS in terms of various parameters of interest.In figure 15, our comparison is based on the mean velocity profiles U * and DR (figure 15a), as well as the Reynolds stress profiles due to the phase-averaged spanwise velocity w2 * xt (figure 15b).We use w2 * xt to calculate the protrusion height by the Stokes layer ( § 3.4).Figure 15 shows that U * , DR and w2 * xt are predicted reasonably well with the coarse LES grid (∆ + x × ∆ + z ≃ 22 × 31).We also concluded in Appendix A that the coarse LES grid predicts DR quite well.Therefore, we employ the coarse LES grid (∆ + x × ∆ + z ≃ 22 × 31) to produce the maps of DR (figure 3), and study the mean velocity profiles (figure 4), and the protrusion height by the Stokes layer (figure 7).
However, studying the turbulent stresses and their spectra, requires the fine LES grid

(∆ +
x × ∆ + z ≃ 14 × 21) as evidenced by figures 16 and 17.In figures 16(a,b), we compare coarse LES with DNS (figure 16a), and fine LES with DNS (figure 16b).Our comparison is based on the one-dimensional pre-multiplied spectrogram for the fluctuating streamwise velocity k * z φ * u ′′ u ′′ (λ * z , y * ).The coarse LES spectrogram (red contour lines in figure 16a) is highly distorted for λ * z 100.This is due to the aliasing error that energises the scales near the cut-off wavelength (Kravchenko & Moin 1997;Park et al. 2004).The aliasing error is clearer from the two-dimensional premultiplied spectrogram k * x k * z φ * u ′′ u ′′ (λ * x , λ * z ) at y * ≃ 20 (figure 16e); the coarse LES spectrogram (red contour lines) agrees well with the DNS spectrogram (filled contour) above the breaking grey line.However, below the grey line, the energy in the LES spectrogram starts to rise, while it must fall following the DNS spectrogram.
raw LES data (figure 17a) with the dealiased u ′′2 * xzt profiles (figure 17b), highlights the improvement due to dealiasing, especially for the fine LES case (blue line in figure 17b).
Overall, we conclude that the coarse LES grid (∆ + x × ∆ + z ≃ 22 × 31) is suitable for studying DR, mean velocity profiles U * and w2 * xt (for the Stokes layer dynamics).The fine LES resolution (∆ + x × ∆ + z ≃ 14 × 21) with dealiasing is more suitable for studying the Reynolds stress profiles and their spectrograms.

Appendix C. Domain size study
In figure 13 we obtained very good agreement in DR between the medium-domain simulation and the full-domain simulation for the case at Re τ = 4000 with A + = 12, κ + x = 0.0014, ω + = −0.044.Here, we further study the domain size effect for some of our production cases at Re τ = 4000 (table 4).We aim to show that the medium domain size is suitable for our parameter space of interest.We select three cases with (κ + x , ω + ) = (0.021, −0.1), (0.021, +0.1), (0.007, +0.05).The cases with κ + x = 0.021 fall at the upper bound of our range of interest for κ + x , and the case with κ + x = 0.007 falls within this range.Also, we consider cases with upstream travelling waves (ω + < 0) and downstream travelling waves (ω + > 0).We deliberately choose the case with (κ + x , ω + ) = (0.007, +0.05), because the wall actuation disturbs the flow to the highest extent ( § 3, figures 4, 6, 7).In fact, this is the most challenging case for the application of the medium domain size among our production cases (table 1).For each case, we perform LES with the medium domain size (figure 1a, 2.0h × 0.6h, y + res ≃ 1000) and the large domain size (figure 1b, 4.0h × 1.2h, y + res ≃ 2000).We report the obtained DR for each case in table 4. The agreement in DR between the medium domain and the large domain is quite good for all cases (within 1.6% difference).We compute DR (hence C f ), following § 2.3.First, we reconstruct the U * profile beyond y * res using Nagib & Chauhan (2008)'s composite profile, indicated with dashed line in figure 18.Then, we obtain C f ≡ 2/U * b 2 by integrating the resolved portion of the U * profile up to y * res (solid line in figure 18) and its reconstructed portion beyond y * res (dashed line in figure 18).Therefore, for the medium domain, C f is obtained by integrating the resolved U * profile up to y * res ≃ 1000 and the reconstructed part beyond that.However, for the large domain size, the integrated U * profile consists of the resolved portion up to y * res ≃ 2000 and the reconstructed portion beyond that.The close agreement in DR between the medium domain and the large domain (table 4), indicates the suitability of the medium domain size (hence sufficiency of resolving up to y * res ≃ 1000).Beyond y * res ≃ 1000 can be accurately reconstructed with the composite profile.Further support for the suitability of the medium domain size is provided in figure 18.We compare the profiles of the mean velocity U * and the velocity difference ∆U * between the medium domain size (red solid line) and the large domain size (blue solid line) for two cases from table 4; κ + x = 0.021, ω + = +0.1 (figure 18a) and κ + x = 0.007, ω + = +0.05(figure 18b).For both actuated cases, the resolved portion of the profiles agree well between the medium domain and the large domain.We observe this agreement in the U * and ∆U * profiles (the insets).For both cases, the logarithmic U * profile appears by y * ≃ 200.This allows to use the composite profile beyond y * ≃ 200.Overall, we conclude that the medium domain size (L x × L z ≃ 2.0h × 0.6h, figure 1a) is suitable for our production simulations at Re τ = 4000 (table 1).

Figure 1 .
Figure 1.Various domain sizes for LES in a channel configuration -(a) medium 2.0h × 0.6h, (b) large 4.0h × 1.2h and (c) full 6.6h × 3.2h.For each domain size, the instantaneous streamwise velocity (u) field is visualised at about 15 viscous units above the bottom wall.The grey shaded zones indicate the wall-heights up to which the flow is resolved for each domain size (yres ≃ 0.4Lz, Chung et al. 2015).

Figure 2 .
Figure 2. (a) Profiles of the mean velocity U * for the LES of the actuated case at Reτ = 4000, A + = 12, κ + x = 0.02 and ω + = −0.05( , ), and LES of the non-actuated case at Reτ = 4000 ( , ).The viscous-scaled quantities U * and y * are scaled by the actual values of uτ for each case.The resolved portion of each LES profile (y * 750) is shown with a solid line, and the unresolved portion (y * 750) is shown with a dashed-dotted line.The unresolved portion of each profile appears as a fictitious wake and is due to the medium domain size (figure 1a).We reconstruct the unresolved portion using the composite profile for channel flow by Nagib & Chauhan (2008) (the dashed lines for y * 750).We compare the resolved ( ) and reconstructed ( ) portions of the non-actuated LES with the DNS of Lozano-Durán & Jiménez (2014) at Reτ = 4200 ( ).(b) Difference between the actuated and non-actuated profiles ∆U * = U * act − U * non-act (blue and black profiles in a) up to the maximum resolved height y * res ≃ 750.To reconstruct the actuated profile beyond y * res ≃ 750 using the composite profile suggested by Nagib & Chauhan (2008), we set the log-law shift ∆B as the value of ∆U * at y * res .

Figure 3 .
Figure 3. (a,b) Maps of DR for A + = 12 at (a) Reτ = 951 and (b) Reτ = 4000.The local maximum DR ( ) and the local minimum DR ( ) for κ + x > 0 are indicated for clarity.We label the region on the left side of ( ) with I, between ( ) and ( ) with II and the right side of ( ) with III.(c) Map of the difference in DR between Reτ = 4000 and Reτ = 951.(d ) Map of the difference in DR between Reτ = 4000 and GQ's prediction (Gatti & Quadrio 2016) at the same Reynolds number.In all (a,b,c,d ), the contour fields and the contour lines show the same quantity.For (a,b) the contour lines grow from −20% to 40%, and for (c,d ) the contour lines grow from −7% to +7%.

Figure 4 .
Figure 4. Variation of DR and the mean velocity profiles U * at Reτ = 4000, A + = 12, κ + x = 0.007, and −0.2 ≤ ω + ≤ +0.2.(a) Variation of DR with ω + ; the inset shows the location of the data points on the DR map.The lines ( ) and ( ) are the local maximum and minimum DR.(b,c) Variation of the U * profiles with ω + for (b) upstream travelling wave (ω + ≤ 0) and (c) downstream travelling wave (ω + ≥ 0); the profile ( ) corresponds to the non-actuated case and the profiles with no symbol correspond to the actuated cases.For each profile, the solid line is the resolved portion and the dashed line is the reconstructed portion following Nagib & Chauhan (2008).For each case, the colour of its U * profile in (b,c) is consistent with the colour of its DR datapoint in (a).In (b,c), the inset plots the same profiles in terms of U * /y * versus y * .(d,e) Diagnostic function y * dU * /dy * for the profiles in (b,c); the inset shows the velocity difference ∆U * = U * act − U * non-act between each actuated profile U * act and the non-actuated profile U * non-act .
in Chan et al. 2015 or figure 3 in MacDonald et al. 2017), while in turbulent flows over riblets ∆U * is almost constant for y *

Figure 5 .
Figure 5. Profiles of Reynolds stresses and turbulence production for four cases from figure 4 at Reτ = 4000, A + = 12, κ + x = 0.007 and ω + = −0.05,0, +0.05, +0.20.The insets in (a,b) indicate the considered values of ω + and their drag reduction values.Line and symbol colours are consistent with figure 4. In each panel, only the resolved portion of the profiles are shown (y + 1000, y * 700).The black lines with symbols correspond to the non-actuated case.(a,c,e) plot the actuated profiles (dashed-dotted lines) scaled by the non-actuated uτ 0 (superscripted with +); (b,d,f ) plot the actuated profiles (solid lines) scaled by the actuated uτ (superscripted with * ).(a-d ) Reynolds stress profiles for (a,b) the streamwise velocity u ′2 xzt and (c,d ) the spanwise velocity w ′2 xzt .(e,f ) Premultiplied production of turbulent kinetic energy.

Figure 6 .
Figure 6.Profiles of velocity statistics at Reτ = 4000 for given cases as in figures 4 and 5 (A + = 12, κ + x = 0.007).Line legends are consistent with figures 4 and 5.In each panel, only the resolved portion of the profiles are shown corresponding to y * ≤ 1000.(a,c,e) correspond to ω + ≤ 0, and (b,d,f ) correspond to ω + > 0. (a,b) U * profiles; the insets indicate the value of ω + and its DR for each profile.In (c-f ) the Reynolds stress profiles are presented in terms of the turbulent component (solid lines) and the harmonic component (dashed lines) following (3.1a,3.1b).(c,d ) turbulent component of the streamwise velocity u ′′2 * xzt .(e,f ) turbulent component w ′′2 * xzt and harmonic component w2 * xt for the spanwise velocity.On each actuated profile, the cross symbol (+) marks the Stokes layer thickness δ * S , and the bullet symbol (•) marks the protrusion height ℓ * 0.01 due to the Stokes layer.

Figure 7 .
Figure 7. (a) Comparison between the map of drag reduction DR (contour field) and the protrusion height by the Stokes layer ℓ * 0.01 (contour lines) for our considered parameter space at Reτ = 4000.(b) Comparison between the map of Stokes layer thickness δ * S (contour field) and ℓ * 0.01 (contour lines) for the same cases as in (a).The lines ( ) and ( ) are the local maximum and minimum DR (same as in figure 3b).(c,d ) plot DR versus ℓ * 0.01 and DR versus δ * S , respectively, for the same data as in (a,b); κ + x = 0.00238 ( ), 0.004 ( ), 0.007 ( ), 0.010 (◮), 0.012 (◭), 0.017 ( ), 0.021 (•).At each κ + x , we plot the cases only in regions I and II (see the map in d ), with the maximum DR case highlighted with a green outline.The grey regions in (c,d ) (20 ≤ ℓ * 0.01 ≤ 30, 5 ≤ δ * S ≤ 7) shade the range of maximum DR at each κ + x .The linear dotted lines in (c,d ) fit the data for ℓ * 0.01 20 (c) and δ * S 5 (d ).The fitting lines also locate the minimum values for ℓ * 0.01,min ≃ 5 (c) and δ * S,min ≃ 1 (d ) to achieve drag reduction.

Figure 10 .
Figure 10.(a) Map of P + in at Reτ = 4000.The filled contour and the line contours show the same quantity.(b) Filled contour is the difference in calculation of P +in at Reτ = 4000 between LES and its theoretical estimation from the generalised Stokes layer (GSL) theory(Quadrio & Ricco 2011); line contours give the Stokes layer protrusion height ℓ * 0.01 (same as in figure7a).

Figure 12 .
Figure 12.(a) Net power saving N P S for LES at Reτ = 4000.(b) Difference in N P S between LES at Reτ = 4000 and LES at Reτ = 951.In each of (a,b), the filled contour and the line contours show the same quantity.All the quantities with + superscript are scaled by ν and the non-actuated uτ o .( ) and ( ) locate the local maximum and the local minimum DR, respectively.

Figure 16 .
Figure 16.Comparison between the coarse LES (∆ + x × ∆ + z ≃ 22 × 31, ), fine LES (∆ + x × ∆ + z ≃ 14 × 21, ) and DNS (∆ + x × ∆ + z ≃ 7 × 4, filled contour).All cases have matched Reτ = 590 and actuation parameters (A + , κ + x , ω + ) = (12, 0.0014, −0.044), see table 3. The comparison is made in terms of (a,b,c,d ) one-dimensional pre-multiplied spectrograms of the turbulent part of the streamwise velocity k * z φ * u ′′ u ′′ (λ * z , y * ), and (e,f ) two-dimensional pre-multiplied spectrograms of the turbulent part of the streamwise velocity k * x k * z φ * u ′′ u ′′ (λ * x , λ * z ) at y * ≃ 20.(a,c,e) are the comparison between the coarse LES and DNS, and (b,d,f ) are the comparison between the fine LES and DNS.(a,b) compare the original spectrograms from the raw LES data (contour lines) with the DNS spectrogram (contour field).(c,d ) compare the dealiased spectrograms from LES (contour lines) with the DNS spectrogram (contour field).Dealiasing is performed through the two-dimensional spectrograms, e.g. by removing the scales below ( ) in (e,f ).See the text for the details.The colourbar for (a-d ) is next to (b), and the colourbar for (e,f ) is next to (f ).

Table 1 .
Summary of the parameters of computational runs.Cases above the separating line are conducted at Reτ = 951 (Re b = 19700), and below the line at Reτ = 4000 (Re b = 94450).
) display u ′′2 * xzt , the stochastic component of the streamwise Reynolds stress.By comparing figure 5b with figures 6c,d, we see that u ′2 * xzt ≃ u ′′2 * xzt , indicating that the harmonic (Stokes layer) component makes a negligible contribution.For the spanwise velocity, however, the harmonic component w2 * xt contributes significantly to the total spanwise Reynolds stress w ′2 * xzt close to the wall (see figures 6e,f ).At y * ∼ O(1), the harmonic component is about three orders of magnitude larger than the turbulent component, while at y * ∼ O(10) they have comparable magnitudes.Figures 6(e,f ) indicate that the rate of decay in w2 * xt , hence the protrusion of the Stokes layer, strongly depends on ω + .Further, the level of distortion in the U * profiles (figures 6a,b) strongly depends on the rate of decay in w2 * xt .Interestingly, in region II (figure Quadrio & Ricco (2011)d lines in figures 7c,d ).FollowingQuadrio & Ricco (2011), if we extrapolate the linear fits to DR = 0, we obtain ℓ * 0.01,min ≃ 5 and δ * S,min ≃ 1; these values indicate the minimum limits for drag reduction to occur.The linear relation between DR, ℓ * 0.01 and δ * S is limited to region I.In region II when DR drops, this linear relation is broken.The observed trends for DR versus δ * ′′ u ′′ is lifted away from the wall to a y * distance that coincides with ℓ * 0.01 .However, for the cases with ℓ * 0.01 > 30, the energetic peaks in k * z φ * u ′′ u ′′ (figures 8m,n) and k * z φ * w ′′ w ′′ (figures 8r,s) instead reside near the wall at y * ≃ 10, well below ℓ * 0.01 ≃ 90.It appears that when ℓ * 0.01 > 30, a near-wall cycle of streaks with high turbulence activity is generated within the Stokes layer.Contrast this behaviour to the case when ℓ * 30, the Stokes layer becomes excessively strong.In this situation, a near-wall cycle of turbulence is generated at y * ≃ 10 that meanders following the Stokes motion.As a result, DR drops by increasing ℓ * 0.01 .

Table 2 .
Summary of the LES cases for validation.The cases above the separating line have fixed actuation parameters A + , κ +x , ω + and grid resolution ∆ + Table 2 lists the simulation details for LES.Both DNS and LES cases are compared at matched Reτ = 951, A + = 12, and κ + x , ω + .(a,b)show the comparison at κ + x = 0.0347 and 0.0208, respectively.At each value of κ +x , comparison is made at six values of ω + (listed in table 2).