[

colorlinks=true, linkcolor=blue, urlcolor=blue, citecolor=blue]hyperref

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

A. Rouhi1<EMAIL_ADDRESS>M. K. Fu2 D. Chandran3 A. Zampiron4 A. J.
Smits5 I. Marusic3 1Department of Engineering, School of Science and
Technology
Nottingham Trent University, Nottingham NG11 8NS, United Kingdom
2Graduate Aerospace Laboratories (GALCIT), Caltech, Pasadena, CA 91125, USA
3Department of Mechanical Engineering, University of Melbourne, Victoria 3010,
Australia
4School of Engineering, University of Aberdeen,
King’s College, Aberdeen AB24 3FX, United Kingdom
5Department of Mechanical and Aerospace Engineering, Princeton University,
Princeton, NJ 08544, USA

###### Abstract

Turbulent drag reduction through streamwise travelling waves of 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 drag reduction at friction Reynolds numbers $Re_{\tau}=951$ and $4000$.
The actuation parameter space is restricted to the inner-scaled actuation
(ISA) pathway, where drag reduction is achieved through direct attenuation of
the near-wall scales. The level of turbulence attenuation, hence drag
reduction, 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 drag reduction 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 drag reduction 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 drag reduction savings. This motivates the study
of the outer-scaled actuation (OSA) pathway in Part 2, where drag reduction is
achieved through actuating the outer-scaled motions.

###### keywords:

turbulence simulation, turbulence control, drag reduction

##  1 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 wall-bounded 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

$w_{s}(x,t)=A\sin{(\kappa_{x}x-\omega t)},$ (1)

where $w_{s}$ is the instantaneous spanwise velocity of the wall surface, $A$
is the amplitude of the spanwise forcing, $\omega$ is the angular frequency of
oscillation, and $\kappa_{x}=2\pi/\lambda$ is the wavenumber of the travelling
wave with wavelength $\lambda$. Negative frequencies result in an upstream
travelling wave, and vice versa. With an appropriate choice of $A,\kappa_{x}$
and $\omega$, 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) 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

$DR=\frac{C_{f_{0}}-C_{f}}{C_{f_{0}}},$ (2)

where $C_{f}\equiv 2\overline{\tau_{w}}/(\rho U^{2}_{b,\infty})$ and
$C_{f_{0}}\equiv 2\overline{\tau_{w_{0}}}/(\rho U^{2}_{b,\infty})$ are the
skin-friction coefficients of the drag-reduced flow (with wall shear-stress
$\overline{\tau_{w}}$) and the non-actuated flow (with wall shear-stress
$\overline{\tau_{w_{0}}}$) and $\rho$ is the fluid density. The overbar in
$\overline{\tau_{w}}$ and $\overline{\tau_{w_{0}}}$ 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_{\infty}$ (Part 2, Bird et al. 2018).
Accordingly, there are two friction velocities
$u_{\tau}\equiv\sqrt{\overline{\tau_{w}}/\rho}$ and
$u_{\tau_{0}}\equiv\sqrt{\overline{\tau_{w_{0}}}/\rho}$, corresponding to the
drag-reduced and non-actuated 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_{\tau_{0}}$ are denoted by
the ‘$+$’ superscript, and those normalised by $u_{\tau}$ are denoted by the
‘$*$’ superscript. The friction Reynolds number $Re_{\tau}$ in a channel flow
(Part 1) is defined based on $u_{\tau_{0}}$ and the channel half-height $h$
$(Re_{\tau}\equiv u_{\tau_{0}}h/\nu)$. In a boundary layer (Part 2), is
defined based on $u_{\tau_{0}}$ and the boundary layer thickness $\delta$
$(Re_{\tau}\equiv u_{\tau_{0}}\delta/\nu)$. By dimensional analysis (Gatti &
Quadrio, 2016; Marusic et al., 2021) we obtain

$DR=DR\left(\kappa^{+}_{x},\omega^{+},A^{+},Re_{\tau}\right),$ (3)

where
$\kappa^{+}_{x}=\kappa_{x}\nu/u_{\tau_{0}},\,\omega^{+}=\omega\nu/u^{2}_{\tau_{0}}$
and $A^{+}=A/u_{\tau_{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) to demonstrate that the introduction of a streamwise
travelling wave achieves higher $DR$ than a purely oscillating wall mechanism
($\kappa_{x}=0$). They fixed $Re_{\tau}=200$ and $A^{+}=12$, and populated a
map of $DR(\omega^{+},\kappa^{+}_{x})$ for $0\leq\kappa^{+}_{x}\leq+0.04$ and
$-0.3\leq\omega^{+}\leq+0.3$. Gatti & Quadrio (2016) extended this work to
$Re_{\tau}=1000$ and a broader range of actuation parameters
($0\leq\kappa^{+}_{x}\leq+0.05,-0.6\leq\omega^{+}\leq+0.6$ and $3\leq
A^{+}\leq 15$) to construct isosurfaces of
$DR(\omega^{+},\kappa^{+}_{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, $\Delta B$, in the logarithmic region as:

$U^{*}=\frac{1}{\kappa}\ln(y^{*})+B+\Delta B,$ (4)

where $U^{*}\equiv U/u_{\tau}$ and $y^{*}\equiv yu_{\tau}/\nu$ are the
viscous-scaled velocity and wall distance, $\kappa$ 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

$\Delta
B=\sqrt{\frac{2}{C_{f_{0}}}}\left[\left(1-DR\right)^{-1/2}-1\right]-\frac{1}{2\kappa}\ln{\left(1-DR\right)}.$
(5)

In this framework, the Reynolds number dependence of the flow is captured by
$C_{f_{0}}$, 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 $\Delta B$; this parameter is
independent of Reynolds number and can be parameterised by the dimensionless
actuation parameters so that $\Delta B=\Delta
B(\kappa^{*}_{x},\,\omega^{*},\,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_{\tau}$, 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. $\kappa_{x}=0$ in 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 (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_{\tau}\lesssim 1500$.

Marusic et al. (2021) recently investigated the parameter space (3) at much
higher Reynolds numbers by conducting experiments up to $Re_{\tau}=12800$ and
wall-resolved large-eddy simulations (LES) up to $Re_{\tau}=2000$. By covering
such a large Reynolds number range, they were able to explore the increasing
contribution of turbulent scales in the log-region 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 near-wall cycle and the near-wall
peak in turbulent kinetic energy. For example, $\omega^{+}\approx-0.06$
equates to a time period of oscillation of
$T^{+}_{osc}=2\pi/\left|\omega^{+}\right|=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}\gg 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 (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_{\tau}$, $\kappa^{+}_{x}$ and
$\omega^{+}$. 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_{\tau}\approx 1000$, generating a new map of $DR$ at $Re_{\tau}=4000$ over
$0.002\leq\kappa^{+}_{x}\leq 0.02$ and $-0.2\leq\omega^{+}\leq+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_{\tau}=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 near-wall 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_{\tau}=4000$ over the range of parameters
considered here.

##  2 Numerical flow setup

###  2.1 Governing equations and solution method

We solve the filtered equations for a channel flow (figure 1) of an
incompressible fluid with constant density $\rho$ and kinematic viscosity
$\nu$

$\frac{\partial\widehat{u}_{i}}{\partial
x_{i}}=0,\quad\frac{\partial\widehat{u}_{i}}{\partial
t}+\frac{\partial\widehat{u}_{i}\widehat{u}_{j}}{\partial
x_{j}}=-\frac{1}{\rho}\frac{\partial\widehat{p}}{\partial
x_{i}}+\nu\frac{\partial^{2}\widehat{u}_{i}}{\partial
x^{2}_{j}}-\frac{\partial\tau_{ij}}{\partial x_{j}}+G\delta_{i1}$ (2.1a,b)

The hat $\widehat{(...)}$ 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
$\widehat{u}_{1},\widehat{u}_{2}$ and $\widehat{u}_{3}$ (or
$\widehat{u},\widehat{v}$ and $\widehat{w}$), respectively. The pressure
gradient in (2.1b) is decomposed into the domain and time-averaged driving
part $-\rho G$, and the periodic (fluctuating) part
$\partial\widehat{p}/\partial x_{i}$. By averaging (2.1b) in time and over the
entire fluid domain, we obtain $G=\overline{\tau_{w}}/(\rho h)=u^{2}_{\tau}/h$
where $h$ is the (open) channel height. $G$ is adjusted based on a target
flowrate (i.e. target bulk Reynolds number $Re_{b}\equiv U_{b}h/\nu$), that is
matched between the actuated and non-actuated cases. The unresolved subgrid-
scale (SGS) stresses
$\tau_{ij}=\widehat{u_{i}u_{j}}-\widehat{u}_{i}\widehat{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).

Figure 1: Various domain sizes for LES in a channel configuration - (a) medium
$2.0h\times 0.6h$, (b) large $4.0h\times 1.2h$ and (c) full $6.6h\times 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 ($y_{\mathrm{res}}\simeq 0.4L_{z}$, Chung et al. 2015).

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 $\widehat{u}=\widehat{v}=0$ and
$\widehat{w}(x,z,t)=A\sin(\kappa_{x}x-\omega t)$, and at the top boundary we
apply free-slip and impermeable conditions ($\partial\widehat{u}/\partial
y=\partial\widehat{w}/\partial y=\widehat{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$.

###  2.2 Simulation cases

Table 1 lists all the simulations completed for $Re_{\tau}=951$ and
$Re_{\tau}=4000$, where $Re_{\tau}\equiv u_{\tau_{0}}h/\nu$ represents the
friction Reynolds number of the non-actuated case. At each $Re_{\tau}$, a
parametric sweep of $7\times 8$ combinations of streamwise wavenumber
($\kappa^{+}_{x}$) and oscillation frequency ($\omega^{+}$) is conducted over
$0.00238\leq\kappa^{+}_{x}\leq 0.02$ and $-0.2\leq\omega^{+}\leq+0.2$. The
spanwise velocity amplitude is fixed at $A^{+}=12$. The seven non-zero values
of $\omega^{+}$ give the oscillation time periods $T_{osc}^{+}=126$, 63, 42,
and 31 (all within the ISA pathway), and $\omega^{+}=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/\nu$ between the
actuated and non-actuated cases and substituting the respective values of
$C_{f}$ and $C_{f_{0}}$ into (2). We consider matched bulk Reynolds numbers
$Re_{b}=19700$ and $94450$, which correspond to $Re_{\tau}=951$ and $4000$ for
the non-actuated channel flow. Quadrio & Ricco (2011) and Ricco et al. (2012)
compute $C_{f}$ and $C_{f_{0}}$ at matched $Re_{\tau}$ (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_{\tau}$ (constant pressure-gradient), see
Quadrio & Ricco (2011), Quadrio (2011) and Ricco et al. (2012). With matched
$Re_{b}$, $C_{f}$ and $C_{f_{0}}$ are obtained at different $Re_{\tau}$.
However, for our considered parameter space, the maximum $DR$ is about $30\%$,
which leads to a maximum deviation of about $16\%$ in $Re_{\tau}$ between
$C_{f}$ and $C_{f_{0}}$. Another source of difference between matched $Re_{b}$
and matched $Re_{\tau}$ is in the actuation amplitude $A$ (1). With constant
$A^{+}=12$, $A^{*}=12$ for the actuated cases with matched $Re_{\tau}$.
However, with matched $Re_{b}$, $A^{*}>12$ when $DR>0$, and vice versa.
Nevertheless, $DR$ weakly depends on $A^{*}\gtrsim 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_{\tau}$ for our
parameter space.

domain | $Re_{\tau},Re_{b}$ | $y^{+}_{\mathrm{res}}$ | $\kappa^{+}_{x}$ | $\omega^{+}$ | $L_{x}/h,,L_{z}/h$ | $N_{x},N_{y},N_{z}$ | $\Delta^{+}_{x},\Delta^{+}_{z}$
---|---|---|---|---|---|---|---
full | $951,19700$ | $951$ | $0.00238$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $8.33,3.14$ | $360,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.004$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.61,3.14$ | $288,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.007$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.61,3.14$ | $288,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.010$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.61,3.14$ | $288,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.012$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.61,3.14$ | $288,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.017$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.61,3.14$ | $288,48,96$ | $22,31$
full | $951,19700$ | $951$ | $0.021$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $6.29,3.14$ | $256,48,96$ | $23,31$
full | $951,19700$ | $951$ | N/A | N/A | $6.28,3.14$ | $192,64,128$ | $47,23$
medium | $4000,94450$ | $1000$ | $0.00238$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $1.99,0.63$ | $384,192,80$ | $21,31$
medium | $4000,94450$ | $1000$ | $0.004$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $1.96,0.63$ | $384,192,80$ | $20,31$
medium | $4000,94450$ | $1000$ | $0.007$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $2.04,0.63$ | $384,192,80$ | $21,31$
medium | $4000,94450$ | $1000$ | $0.010$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $2.04,0.63$ | $384,192,80$ | $21,31$
medium | $4000,94450$ | $1000$ | $0.012$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $1.96,0.63$ | $384,192,80$ | $20,31$
medium | $4000,94450$ | $1000$ | $0.017$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $2.03,0.63$ | $384,192,80$ | $21,31$
medium | $4000,94450$ | $1000$ | $0.021$ | $0,\pm 0.05,\pm 0.1,+0.15,\pm 0.2$ | $2.04,0.63$ | $384,192,80$ | $21,31$
medium | $4000,94450$ | $1000$ | $0.007$ | $0,\pm 0.05,+0.1,+0.2$ | $2.04,0.63$ | $576,288,120$ | $14,21$
medium | $4000,94450$ | $1000$ | N/A | N/A | $2.04,0.63$ | $576,288,120$ | $14,21$
Table 1: Summary of the parameters of computational runs. Cases above the
separating line are conducted at $Re_{\tau}=951$ ($Re_{b}=19700$), and below
the line at $Re_{\tau}=4000$ ($Re_{b}=94450$). The last row at each
$Re_{\tau}$ indicated with N/A for $\kappa^{+}_{x}$ and $\omega^{+}$,
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_{\tau}$, $\kappa^{+}_{x}$ and grid size, but
$\omega^{+}$ is different for each case (as listed in the fifth column). Those
values of $\omega^{+}$ with $\pm$ 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_{\tau}=951$ we use the full domain and at $Re_{\tau}=4000$
we use the medium domain. The third column $y^{+}_{\mathrm{res}}$ is the
maximum resolved height by the simulation domain ($\simeq 0.4L^{+}_{z}$, Chung
et al. 2015). The eighth row at $Re_{\tau}=4000$ repeats some of the cases
with $\kappa^{+}_{x}=0.007$ (the third row at $Re_{\tau}=4000$), but with a
finer grid resolution.

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_{\tau}\approx 1000$,
experimental data of Marusic et al. (2021) at $Re_{\tau}=6000$, and our self-
generated DNS data at $Re_{\tau}=590$. For $DR$ and the mean velocity profile,
we used the same viscous-scaled grid resolution at $Re_{\tau}=951$ and $4000$,
corresponding to the streamwise and spanwise grid sizes of
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 21\times 31$ (the first seven rows
at each $Re_{\tau}$ 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_{\tau}=4000$ we used a finer grid resolution with
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 21$ (the last two rows in
table 1). Our nominal LES filter width
$\Delta^{+}=\left(\Delta^{+}_{x}\Delta^{+}_{y}\Delta^{+}_{z}\right)^{1/3}$ is
$7\lesssim\Delta^{+}\lesssim 34$ for the coarser grid, and
$5\lesssim\Delta^{+}\lesssim 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^{+}_{\Delta_{z}}\simeq 0.12$ and $0.18$. These wavenumbers are $6$ and $9$
times larger than our maximum actuation wavenumber $\kappa^{+}_{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^{+}\simeq 10$) with the smallest energetic eddies. If we take the
convective speed of $10u_{\tau_{0}}$ in this region, our cut-off frequencies
are $\omega^{+}_{\Delta_{z}}\simeq 1.2$ and $1.8$ for the coarser and finer
grids, respectively, which are $6$ and $9$ times larger than our maximum
actuation frequency $\omega^{+}=\pm 0.2$.

In terms of the domain size, the cases at $Re_{\tau}=951$ used a full-domain
with $L_{x}\times L_{z}\simeq 6.6h\times 3.2h$ (figure 1c), which is
sufficiently large to resolve the first and second-order statistics across the
entire channel (Lozano-Durán & Jiménez, 2014). However, at $Re_{\tau}=4000$
each full-domain calculation is about $500$ times more expensive than that at
$Re_{\tau}=951$, and so the domain size was reduced to $L_{x}\times
L_{z}\simeq 2.0h\times 0.6h$ (figure 1a). As a consequence, the flow is only
resolved up to a fraction of the channel height $y^{+}_{\mathrm{res}}\simeq
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^{+}_{\mathrm{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), the prescriptions are $L^{+}_{z}\simeq
2.5y^{+}_{\mathrm{res}},L^{+}_{x}\gtrsim\max(3L^{+}_{z},1000,\lambda^{+})$,
where $\lambda$ is the travelling wavelength. MacDonald et al. (2017, 2018)
used the reduced-domain approach with $60\lesssim y^{+}_{\mathrm{res}}\lesssim
250$ for turbulent flows over roughness. Endrikat et al. (2021b) used the same
approach with $y^{+}_{\mathrm{res}}\simeq 100$ for turbulent flows over
riblets. Jiménez & Moin (1991) who used this approach for the first time
resolved the flow up to $y^{+}_{\mathrm{res}}\simeq 80$. They named this
approach “minimal flow unit”. Here, with $L_{x}\times L_{z}\simeq 2.0h\times
0.6h$ (figure 1a) at $Re_{\tau}=4000$, we resolve a substantial fraction of
the inner layer up to $y^{+}_{\mathrm{res}}\simeq 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}\times L_{z}\simeq 1.4h\times 0.7h$ with
$y^{+}_{\mathrm{res}}\simeq 250$ to study the travelling wave (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.

###  2.3 Calculation of the skin-friction coefficient

To compute $DR$ (2), we need the skin-friction coefficient $C_{f}\equiv
2\overline{\tau_{w}}/(\rho U^{2}_{b})\equiv 2/{U^{*}_{b}}^{2}$ for both the
actuated and the non-actuated cases. Here,
$U^{*}_{b}=\int_{0}^{h^{*}}U^{*}dy^{*}/h^{*}$ is the viscous-scaled bulk
velocity. For the cases at $Re_{\tau}=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_{\tau}=4000$ with the medium
domain size, the $U^{*}$ profile is resolved only up to
$y^{*}_{\mathrm{res}}\simeq 750-1000$. Two of these high Reynolds number
profiles are shown in figure 2(a): the actuated case with $A^{+}=12$,
$\kappa^{+}_{x}=0.02$ and $\omega^{+}=-0.05$ (blue lines), and the non-
actuated case (black lines). The resolved portion of the LES profile below
$y^{*}_{\mathrm{res}}$ is shown with a solid line, and the unresolved portion
above $y^{*}_{\mathrm{res}}$ with a dashed-dotted line. We also overlay the
DNS of the non-actuated full-domain channel flow at $Re_{\tau}=4200$ by
Lozano-Durán & Jiménez (2014) (red squares). For the non-actuated LES, the
resolved portion up to $y^{*}_{\mathrm{res}}\simeq 1000$ (solid black line)
accurately reproduces the non-actuated DNS. However, the unresolved portion
beyond $y^{*}_{\mathrm{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), and Endrikat et al. (2021b). For accurate
prediction of $U^{*}_{b}$, hence $C_{f}$, it was found that the resolved
height $y^{*}_{\mathrm{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^{*}_{\mathrm{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^{*}_{\mathrm{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) Profiles of the mean velocity $U^{*}$ for the LES of the
actuated case at $Re_{\tau}=4000,A^{+}=12,\kappa^{+}_{x}=0.02$ and
$\omega^{+}=-0.05$ ( ,  ), and LES of the non-actuated case at
$Re_{\tau}=4000$ ( ,  ). The viscous-scaled quantities $U^{*}$ and $y^{*}$ are
scaled by the actual values of $u_{\tau}$ for each case. The resolved portion
of each LES profile $(y^{*}\lesssim 750)$ is shown with a solid line, and the
unresolved portion $(y^{*}\gtrsim 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^{*}\gtrsim 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_{\tau}=4200$ (  ). (b) Difference between the
actuated and non-actuated profiles $\Delta
U^{*}=U^{*}_{\mathrm{act}}-U^{*}_{\mathrm{non{\text{-}}act}}$ (blue and black
profiles in a) up to the maximum resolved height $y^{*}_{\mathrm{res}}\simeq
750$. To reconstruct the actuated profile beyond $y^{*}_{\mathrm{res}}\simeq
750$ using the composite profile suggested by Nagib & Chauhan (2008), we set
the log-law shift $\Delta B$ as the value of $\Delta U^{*}$ at
$y^{*}_{\mathrm{res}}$.

Figure 2(a) demonstrates that for the non-actuated case at $Re_{\tau}=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^{*}_{\mathrm{res}}$ and the reconstructed
profile beyond $y^{*}_{\mathrm{res}}$. We find that $C_{f_{0}}$ 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
$\Delta B$ in the composite profile to make the resolved and reconstructed
profiles continuous at $y^{*}_{\mathrm{res}}$. We find $\Delta B$ by plotting
the velocity difference between the actuated and non-actuated profiles $\Delta
U^{*}=U^{*}_{\mathrm{act}}-U^{*}_{\mathrm{non{\text{-}}actuated}}$ (figure
2b). As seen in figure 2(b), $\Delta U^{*}$ reaches almost a plateau beyond
$y^{*}\simeq 100$. We set $\Delta B$ as the value of $\Delta U^{*}$ at
$y^{*}_{\mathrm{res}}\simeq 750$. Note that since the actuated $u_{\tau}$ is
smaller than the non-actuated $u_{\tau_{o}}$, $y^{*}_{\mathrm{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^{*}_{\mathrm{res}}$ (solid blue line) and the
reconstructed portion beyond $y^{*}_{\mathrm{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.

##  3 Results

###  3.1 Drag reduction map as a function of frequency and wavelength

Figures 3(a,b) display the maps of $DR(\omega^{+},\kappa^{+}_{x})$ at
$Re_{\tau}=951$ and $4000$ from the computations listed in table 1. At each
$Re_{\tau}$, we have $56$ $DR$ data points. To generate the maps, we perform
bilinear interpolation of our $DR$ data points onto a uniform $20\times 20$
grid over the parameter space $0\leq\kappa^{+}_{x}\leq 0.02$ and
$-0.2\leq\omega^{+}\leq+0.2$. At $Re_{\tau}=951$, the maximum $DR$ of $35.4\%$
at $(\omega^{+},\kappa^{+}_{x})=(0.05,0.021)$ is in close agreement with the
DNS of Gatti & Quadrio (2016) at $Re_{\tau}\simeq 950$, where the maximum $DR$
was found to be $38.8\%$ at $(\omega^{+},\kappa^{+}_{x})=(0.05,0.0195)$. At
$Re_{\tau}=4000$, the maximum $DR$ decreases to $27.5\%$ at the same actuation
parameters $(\omega^{+},\kappa^{+}_{x})=(0.05,0.021)$. At each Reynolds
number, $DR$ changes more drastically by changing $\omega^{+}$ than by
changing $\kappa^{+}_{x}$.

Figure 3: (a,b) Maps of $DR$ for $A^{+}=12$ at (a) $Re_{\tau}=951$ and (b)
$Re_{\tau}=4000$. The local maximum $DR$ ( ) and the local minimum $DR$ ( )
for $\kappa^{+}_{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_{\tau}=4000$ and
$Re_{\tau}=951$. (d) Map of the difference in $DR$ between $Re_{\tau}=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\%$.

When $\kappa^{+}_{x}=0$, there is no travelling wave (plane wall oscillation),
and the variation of $DR$ is symmetric between $\omega^{+}<0$ and
$\omega^{+}>0$. In this case, two equal local maxima (at $\omega^{+}\simeq\pm
0.05$) and a local minimum (at $\omega=0$) emerge. When $\kappa^{+}_{x}>0$, a
travelling wave is generated, and the variation of $DR$ is asymmetric between
$\omega^{+}<0$ and $\omega^{+}>0$. In this case, at each $\kappa^{+}_{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 $\omega^{+}\lesssim 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 $\kappa^{+}_{x}\lesssim 0.007$, the local
minimum $DR$ is positive, however, for $\kappa^{+}_{x}\gtrsim 0.007$ the local
minimum $DR$ becomes negative (hence a drag increase). Increase in
$\kappa^{+}_{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 $\omega^{+}/\kappa^{+}_{x}\simeq 10$. In
other words, maximum drag increase occurs when the travelling wave speed is
about $10u_{\tau_{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
$\omega^{+}/\kappa^{+}_{x}\simeq 10$. Region III covers the right of local
minimum $DR$ (black dashed curve) where $\omega^{+}>0$ (downstream travelling
wave); in this region, increase in $\omega^{+}$ 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 $(\omega^{+},\kappa^{+}_{x})$ space,
$DR$ is lower at the higher Reynolds number. Only within the range
$0.005\lesssim\kappa^{+}_{x}\lesssim 0.020,+0.1\lesssim\omega^{+}\lesssim+0.2$
we observe the opposite trend. This region coincides with the drag increasing
range ($DR<0$) with $\omega^{+}>0$. This observation is consistent with GQ’s
model (5), where $DR<0$ (hence $\Delta 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_{\tau}=4000$. To predict $DR$, the model (5) requires $C_{f_{0}}$ and
the value of the log-law shift $\Delta B$ for each set of actuation parameters
$(A^{+},\kappa^{+}_{x},\omega^{+})$. For $C_{f_{0}}$, 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 $\Delta B$ from a low Reynolds number
simulation for the same set of $(A^{+},\kappa^{+}_{x},\omega^{+})$ because
$\Delta B$ is assumed to be Reynolds number independent. Therefore, we use our
results at $Re_{\tau}=951$, where for each $(\omega^{+},\kappa^{+}_{x})$ we
find $\Delta B$ from the velocity difference $\Delta 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\gtrsim 0.3$ according to Pope 2000). By having $\Delta B$ at each
$(\omega^{+},\kappa^{+}_{x})$ and having $C_{f_{0}}$ at $Re_{\tau}=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\%_{\mathrm{LES,Re_{\tau}=4000}}-DR\%_{\mathrm{GQ,Re_{\tau}=4000}}|\lesssim
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 $\omega^{+}\gtrsim+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.

###  3.2 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_{\tau}=4000,A^{+}=12$ and $\kappa^{+}_{x}=0.007$ for $\omega^{+}$ ranging
from $-0.2$ to $+0.2$. In figure 4(a), we identify the selected values of
$\omega^{+}$ (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 $\omega^{+}\leq 0$ (upstream travelling waves)
up to the local maximum $DR$ (region I) and $\omega^{+}\geq 0$ (downstream
travelling waves) beyond the local maximum $DR$ (regions II, III),
respectively.

Figure 4: Variation of $DR$ and the mean velocity profiles $U^{*}$ at
$Re_{\tau}=4000$, $A^{+}=12$, $\kappa^{+}_{x}=0.007$, and
$-0.2\leq\omega^{+}\leq+0.2$. (a) Variation of $DR$ with $\omega^{+}$; 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 $\omega^{+}$ for (b) upstream travelling wave ($\omega^{+}\leq
0$) and (c) downstream travelling wave ($\omega^{+}\geq 0$); the profile (
$\square$ ) 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 $\Delta
U^{*}=U^{*}_{\mathrm{act}}-U^{*}_{\mathrm{non{\text{-}}act}}$ between each
actuated profile $U^{*}_{\mathrm{act}}$ and the non-actuated profile
$U^{*}_{\mathrm{non{\text{-}}act}}$.

When $\omega^{+}\leq 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 $\kappa^{-1}\simeq 2.5$. For the non-actuated case, the
plateau appears for $100\lesssim y^{*}\lesssim 600$. This range is consistent
with the DNS of channel flow by Lozano-Durán & Jiménez (2014) and Lee & Moser
(2015) at $Re_{\tau}=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 $\Delta 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 (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 $\omega^{+}$ 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 $\omega^{+}=+0.05$.
For this case, the diagnostic function tends towards the plateau
$\kappa^{-1}$, but does not reach it. Similarly, $\Delta U^{*}$ for this case
approaches a plateau of $1.7$ by the resolved height
$y^{*}_{\mathrm{res}}\simeq 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$, $\Delta U^{*}$ needs to
reach a plateau by the resolved height $y^{*}_{\mathrm{res}}\simeq 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^{*}_{\mathrm{res}}\simeq 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 $\Delta 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 $\Delta U^{*}$ is almost
constant for $y^{*}\gtrsim 30$ (e.g. figure 6 in Chan et al. 2015 or figure 3
in MacDonald et al. 2017), while in turbulent flows over riblets $\Delta
U^{*}$ is almost constant for $y^{*}\gtrsim 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 $\omega^{+}$ 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^{*}\simeq 100$
with viscous sublayer thickening.

###  3.3 Turbulence statistics

Figure 5: Profiles of Reynolds stresses and turbulence production for four
cases from figure 4 at $Re_{\tau}=4000,A^{+}=12,\kappa^{+}_{x}=0.007$ and
$\omega^{+}=-0.05,0,+0.05,+0.20$. The insets in (a,b) indicate the considered
values of $\omega^{+}$ 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^{+}\lesssim 1000,y^{*}\lesssim 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_{\tau_{0}}$ (superscripted with $+$); (b,d,f) plot the actuated profiles
(solid lines) scaled by the actuated $u_{\tau}$ (superscripted with $*$).
(a–d) Reynolds stress profiles for (a,b) the streamwise velocity
$\left<u^{\prime 2}\right>_{xzt}$ and (c,d) the spanwise velocity
$\left<w^{\prime 2}\right>_{xzt}$. (e,f) Premultiplied production of turbulent
kinetic energy.

We now assess the behaviour of the Reynolds stress distributions at
$Re_{\tau}=4000$ (figures 5a–d) and the turbulent kinetic energy production
$P=-\left<u^{\prime}v^{\prime}\right>_{xzt}dU/dy$ (figures 5e,f), where
$\left<...\right>_{xzt}$ denotes averaging over $xz$-plane and time. We
highlight four cases from figure 4 ($A^{+}=12,\ \kappa^{+}_{x}=0.007$), where
we vary $\omega^{+}$ 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_{\tau_{0}}$ (dashed-dotted lines, figures 5a,c) and by the actuated
$u_{\tau}$ (solid lines, figures 5b,d). Scaling by $u_{\tau_{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).

Scaling by the non-actuated $u_{\tau_{0}}$ (‘$+$’ superscript), as in figures
5(a,c,e), indicates that the wall oscillation attenuates the $\langle
u^{\prime 2}\rangle^{+}_{xzt}$ levels up to the resolved height
$y^{+}_{\mathrm{res}}\simeq 1000$. The cases with the highest $DR$
($\omega^{+}=-0.05,0$ in figure 5a) show the highest level of attenuation in
their $\langle u^{\prime 2}\rangle^{+}_{xzt}$. Additionally, for these cases
the inner peak of $\langle u^{\prime 2}\rangle^{+}_{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
$\langle u^{\prime 2}\rangle^{+}_{xzt}$, the $\langle w^{\prime
2}\rangle^{+}_{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_{\tau}$ (‘$*$’ superscript) is equivalent to inner
scaling, which highlights the extent up to which the actuated profiles depart
from the non-actuated profile. For $\langle u^{\prime 2}\rangle^{*}_{xzt}$ and
$\langle w^{\prime 2}\rangle^{*}_{xzt}$ (figures 5b,d), the actuated cases
agree with the non-actuated case at distances far from the wall, but near the
wall the actuated $\langle u^{\prime 2}\rangle^{*}_{xzt}$ levels are
attenuated, while the $\langle w^{\prime 2}\rangle^{*}_{xzt}$ levels are
amplified. For $\omega^{+}=+0.05$ (in region II), the point where the actuated
profiles begin to depart from the non-actuated counterpart occurs at
$y^{*}\simeq 100$, considerably farther than for the other cases
$(y^{*}\lesssim 30)$. The same case yields the strongest level of near-wall
amplification for $\langle w^{\prime 2}\rangle^{*}_{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 $\langle w^{\prime 2}\rangle_{xzt}$ is
amplified near the wall, $\langle u^{\prime 2}\rangle_{xzt}$ is attenuated,
the viscous sublayer is thickened and $DR$ is increased. This trend occurs in
regions I ($\omega^{+}=-0.05,0$) and III ($\omega^{+}=+0.2$). In region II
($\omega^{+}=+0.05$), however, there is an excessive amplification of
$\left<w^{\prime 2}\right>_{xzt}$ near the wall, a thinning of the viscous
sublayer and a drop in $DR$.

###  3.4 Stokes layer: an important source of inner-scaled actuation

As indicated earlier, the near-wall amplification of $\left<w^{\prime
2}\right>_{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_{\tau}$ scaling, as we are interested in the level of
departure from the non-actuated behaviour. In Part 2, we mostly use
$u_{\tau_{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), the instantaneous
flow can be triply decomposed similar to Touber & Leschziner (2012), as in

$\displaystyle f(x,y,z,t)$
$\displaystyle=\left<f\right>_{xzt}(y)+\underbrace{\tilde{f}(x,y,t)+f^{\prime\prime}(x,y,z,t)}_{f^{\prime}(x,y,z,t)}$
(3.1a) $\displaystyle\tilde{f}(x,y,t)$
$\displaystyle=\frac{1}{N}\sum_{n=0}^{N-1}\left<f\right>_{z}(x,y,t+nT_{osc})-\left<f\right>_{xzt}(y)$
(3.1b) $\displaystyle\left<f^{\prime 2}\right>_{xzt}$
$\displaystyle=\left<\tilde{f}^{2}\right>_{xt}+\left<f^{\prime\prime
2}\right>_{xzt}$ (3.1c)

where $f$ indicates the quantity of interest, i.e. $u,v$ or $w$. In (3.1a),
the total fluctuations $f^{\prime}$ is decomposed into the harmonic
contribution $\tilde{f}$ and the stochastic (turbulent) contribution
$f^{\prime\prime}$. The harmonic contribution $\tilde{f}$ is obtained by phase
averaging the spanwise averaged field $\left<f\right>_{z}$ in time over the
number of periods $N$, and then subtracting the mean vertical profile
$\left<f\right>_{xzt}$. Accordingly, the total Reynolds stress
$\left<f^{\prime 2}\right>_{xzt}$ is decomposed into its harmonic component
$\left<\tilde{f}^{2}\right>_{xt}$ associated with the Stokes layer dynamics
and its turbulent (stochastic) component $\left<f^{\prime\prime
2}\right>_{xzt}$ (3.1c).

Figure 6: Profiles of velocity statistics at $Re_{\tau}=4000$ for given cases
as in figures 4 and 5 ($A^{+}=12,\kappa^{+}_{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^{*}\leq 1000$. (a,c,e) correspond
to $\omega^{+}\leq 0$, and (b,d,f) correspond to $\omega^{+}>0$. (a,b) $U^{*}$
profiles; the insets indicate the value of $\omega^{+}$ 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
$\left<u^{\prime\prime 2}\right>^{*}_{xzt}$. (e,f) turbulent component
$\left<w^{\prime\prime 2}\right>^{*}_{xzt}$ and harmonic component
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ for the spanwise velocity. On each
actuated profile, the cross symbol (+) marks the Stokes layer thickness
$\delta^{*}_{S}$, and the bullet symbol ($\bullet$) marks the protrusion
height $\ell^{*}_{0.01}$ due to the Stokes layer.

In figure 6, we plot these two components for the cases given in figures 4 and
5 ($A^{+}=12,\kappa^{+}_{x}=0.007$, $Re_{\tau}=4000$). For reference, figures
6(a,b) show the considered $U^{*}$ profiles (as in figures 4b,c). Figures
6(c,d) display $\left<u^{\prime\prime 2}\right>^{*}_{xzt}$, the stochastic
component of the streamwise Reynolds stress. By comparing figure 5b with
figures 6c,d, we see that $\left<u^{\prime
2}\right>^{*}_{xzt}\simeq\left<u^{\prime\prime 2}\right>^{*}_{xzt}$,
indicating that the harmonic (Stokes layer) component makes a negligible
contribution. For the spanwise velocity, however, the harmonic component
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ contributes significantly to the total
spanwise Reynolds stress $\left<w^{\prime 2}\right>^{*}_{xzt}$ close to the
wall (see figures 6e,f). At $y^{*}\sim\mathcal{O}(1)$, the harmonic component
is about three orders of magnitude larger than the turbulent component, while
at $y^{*}\sim\mathcal{O}(10)$ they have comparable magnitudes. Figures 6(e,f)
indicate that the rate of decay in $\left<\tilde{w}^{2}\right>^{*}_{xt}$,
hence the protrusion of the Stokes layer, strongly depends on $\omega^{+}$.
Further, the level of distortion in the $U^{*}$ profiles (figures 6a,b)
strongly depends on the rate of decay in
$\left<\tilde{w}^{2}\right>^{*}_{xt}$. Interestingly, in region II (figure 6f)
the decay rate in $\left<\tilde{w}^{2}\right>^{*}_{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
$\left<\mathcal{K}\right>_{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 $\delta^{*}_{S}$ that is featured in Stokes’
second problem (Batchelor, 2000). Following Quadrio & Ricco (2011), we define
$\delta^{*}_{S}$ as the height $y^{*}$ where the amplitude of $\tilde{w}$
decays to $Ae^{-1}$ (i.e. where
$\left<\tilde{w}^{2}\right>^{*}_{xt}=\textstyle{\frac{1}{2}}{A^{*}}^{2}e^{-2}$).
In figure 6, we mark $\delta^{*}_{S}$ on each profile with a cross symbol. The
second length scale $\ell^{*}_{0.01}$ is new, and it is defined as the height
where $\left<\tilde{w}^{2}\right>^{*}_{xt}=0.01$. Our choice for the threshold
of $\left<\tilde{w}^{2}\right>^{*}_{xt}=0.01$ is based on the observation that
$\left<w^{\prime\prime 2}\right>^{*}_{xzt}\sim\mathcal{O}(1)$ in the buffer
and log regions (also reported by Lee & Moser 2015 and Baidya et al. 2021). In
other words, we define $\ell^{*}_{0.01}$ as the height where the Stokes layer
stress $\left<\tilde{w}^{2}\right>^{*}_{xt}$ drops to about $1\%$ of the
spanwise turbulent stress $\left<w^{\prime\prime 2}\right>^{*}_{xzt}$. In
figure 6, we mark $\ell^{*}_{0.01}$ on each profile with a bullet symbol.

The key difference between $\delta^{*}_{S}$ and $\ell^{*}_{0.01}$ is that we
mark $\delta^{*}_{S}$ where the Stokes layer stress
$\left<\tilde{w}^{2}\right>^{*}_{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 $\ell^{*}_{0.01}$ where the Stokes layer stress
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ is a small fraction of the turbulent
stress $\left<w^{\prime\prime 2}\right>^{*}_{xzt}$, hence considering the
background turbulence in this definition. In figures 6(c–f), $\ell^{*}_{0.01}$
coincides well with the distance where the actuated $\left<u^{\prime\prime
2}\right>^{*}_{xzt}$ and $\left<w^{\prime\prime 2}\right>^{*}_{xzt}$ profiles
depart from the non-actuated counterpart. However, $\delta^{*}_{S}$
underestimates the actual protrusion by the Stokes layer due to its ignorance
of the background turbulence. For instance, for the case with $\omega^{+}=0$
at $y^{*}=\delta^{*}_{S}$ (black cross symbol in figure 6e),
$\left<\tilde{w}^{2}\right>^{*}_{xt}\simeq 8\left<w^{\prime\prime
2}\right>^{*}_{xzt}$, i.e. the Stokes layer is $8$ times stronger than the
background turbulence. However, at $y^{*}=\ell^{*}_{0.01}$ (black bullet
symbol) $\left<\tilde{w}^{2}\right>^{*}_{xt}\simeq 0.01\left<w^{\prime\prime
2}\right>^{*}_{xzt}$, i.e. the Stokes layer is $100$ times weaker than the
background turbulence. We propose, therefore, that $\ell^{*}_{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
$\ell^{*}_{0.01}$, as well as the departure height in the
$\left<u^{\prime\prime 2}\right>^{*}_{xzt}$ and $\left<w^{\prime\prime
2}\right>^{*}_{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^{*}\simeq 100$ with viscous sublayer
thickening. However, in region II there is a large increase in
$\ell^{*}_{0.01}$ and the departure in the $\left<u^{\prime\prime
2}\right>^{*}_{xzt}$ and $\left<w^{\prime\prime 2}\right>^{*}_{xzt}$ profiles
also starts at larger distance from the wall. For example, for
$\omega^{+}=+0.05$ in region II ( figures 6d,f), $\ell^{*}_{0.01}\simeq 80$,
which also closely marks the point where the actuated $\left<u^{\prime\prime
2}\right>^{*}_{xzt}$ and $\left<w^{\prime\prime 2}\right>^{*}_{xzt}$ profiles
depart from their non-actuated counterpart. As a result, the mean velocity
profile for $\omega^{+}=+0.05$ in region II (figure 6b) is highly distorted up
to $y^{*}\simeq 200-300$.

Figure 7: (a) Comparison between the map of drag reduction $DR$ (contour
field) and the protrusion height by the Stokes layer $\ell^{*}_{0.01}$
(contour lines) for our considered parameter space at $Re_{\tau}=4000$. (b)
Comparison between the map of Stokes layer thickness $\delta^{*}_{S}$ (contour
field) and $\ell^{*}_{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 $\ell^{*}_{0.01}$ and $DR$ versus
$\delta^{*}_{S}$, respectively, for the same data as in (a,b);
$\kappa^{+}_{x}=0.00238$ (

$\blacksquare$

), $0.004$ ($\blacktriangle$), $0.007$ (  ), $0.010$ ($\blacktriangleright$),
$0.012$ ($\blacktriangleleft$), $0.017$ ($\blacktriangledown$), $0.021$
($\bullet$). At each $\kappa^{+}_{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\leq\ell^{*}_{0.01}\leq
30,5\leq\delta^{*}_{S}\leq 7)$ shade the range of maximum $DR$ at each
$\kappa^{+}_{x}$. The linear dotted lines in (c,d) fit the data for
$\ell^{*}_{0.01}\lesssim 20$ (c) and $\delta^{*}_{S}\lesssim 5$ (d). The
fitting lines also locate the minimum values for
$\ell^{*}_{0.01,\mathrm{min}}\simeq 5$ (c) and
$\delta^{*}_{S,\mathrm{min}}\simeq 1$ (d) to achieve drag reduction.

Furthermore, we can draw a connection between the protrusion height
$\ell^{*}_{0.01}$ and the level of drag reduction. In figures 7(a,b), we
overlay the map of $\ell^{*}_{0.01}$ onto the maps of $DR$ and
$\delta^{*}_{S}$. In region I (left side of the blue dashed-dotted line),
$\ell^{*}_{0.01}\lesssim 30$ and $\delta^{*}_{S}\lesssim 7$. In this region,
an increase in $\ell^{*}_{0.01}$ and $\delta^{*}_{S}$ leads to an increase in
$DR$. For upstream travelling waves ($\omega^{+}<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 $\ell^{*}_{0.01}$ and $\delta^{*}_{S}$. Another
difference between regions I and II, is in the relation between
$\ell^{*}_{0.01}$ and $\delta^{*}_{S}$. In region I, $\ell^{*}_{0.01}$ and
$\delta^{*}_{S}$ are proportional to each other with $\ell^{*}_{0.01}\approx
4\delta^{*}_{S}$. However, in region II this proportional relation is broken
and $\ell^{*}_{0.01}$ can reach as high as $8\delta^{*}_{S}$. At each
$\kappa^{+}_{x}$, the maximum $DR$ (the blue dashed-dotted line) coincides
with the optimal range $20\lesssim\ell^{*}_{0.01}\lesssim 30$
$(5\lesssim\delta^{*}_{S}\lesssim 7)$. In figure 7(c), we plot $DR$ versus
$\ell^{*}_{0.01}$ for our simulation cases in regions I and II. Also,
following Quadrio & Ricco (2011) (their figure 9), we plot $DR$ versus
$\delta^{*}_{S}$ for the same cases (figure 7d). These plots confirm that the
maximum $DR$ coincides with $20\lesssim\ell^{*}_{0.01}\lesssim 30$ and
$5\lesssim\delta^{*}_{S}\lesssim 7$ (shaded in grey). Further, for
$\ell^{*}_{0.01}\lesssim 20$ ($\delta^{*}_{S}\lesssim 5$),
$\ell^{*}_{0.01}\simeq 4\delta^{*}_{S}$ and $DR$ increases linearly with
$\ell^{*}_{0.01}$ and $\delta^{*}_{S}$ (see the fitting dotted lines in
figures 7c,d). Following Quadrio & Ricco (2011), if we extrapolate the linear
fits to $DR=0$, we obtain $\ell^{*}_{0.01,\mathrm{min}}\simeq 5$ and
$\delta^{*}_{S,\mathrm{min}}\simeq 1$; these values indicate the minimum
limits for drag reduction to occur.

The linear relation between $DR,\ell^{*}_{0.01}$ and $\delta^{*}_{S}$ is
limited to region I. In region II when $DR$ drops, this linear relation is
broken. The observed trends for $DR$ versus $\delta^{*}_{S}$ and the minimum
limit of $\delta^{*}_{S,\mathrm{min}}\simeq 1$ (figure 7d) are also reported
by Quadrio & Ricco (2011) (their figure 9). Quadrio & Ricco (2011) calculated
$\delta^{*}_{S}$ from the laminar $\left<\tilde{w}^{2}\right>^{*}_{xt}$
profile based on Stokes layer solution. Here, however, we calculate
$\delta^{*}_{S}$ from the actual $\left<\tilde{w}^{2}\right>^{*}_{xt}$ profile
by phase averaging the simulation data. The variation in $DR$ versus
$\ell^{*}_{0.01}$ is similar to $DR$ versus $\delta^{*}_{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$
($\delta^{*}_{S,\mathrm{min}},\ell^{*}_{0.01,\mathrm{min}}$). However, in
region II when the linear trends are broken, we observe noticeable differences
between $DR$ versus $\delta^{*}_{S}$ and $DR$ versus $\ell^{*}_{0.01}$. In
region II, there does not appear to be a consistent relation between $DR$ and
$\delta^{*}_{S}$ (the red symbols in figure 7d). In other words, we cannot
find a threshold for $\delta^{*}_{S}$ beyond which $DR$ drops. For instance,
for the case with $\kappa^{+}_{x}=0.02,\omega^{+}=+0.2$, $DR$ drops to $-10\%$
but $\delta^{*}_{S}\simeq 5$ which is within the optimal range
$(5\lesssim\delta^{*}_{S}\lesssim 7)$. In contrast, there is a much stronger
connection between $DR$ and $\ell^{*}_{0.01}$, even in region II (figure 7c).
For all cases, increasing $\ell^{*}_{0.01}$ beyond $30$ decreases $DR$. As a
result, the value of $\ell^{*}_{0.01}$ can be used to determine whether we are
in region I ($\ell^{*}_{0.01}\lesssim 30$) or region II
($\ell^{*}_{0.01}\gtrsim 30$).

###  3.5 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
$\left<u^{\prime\prime 2}\right>_{xzt}$ and amplification of
$\left<w^{\prime\prime 2}\right>_{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_{\tau}=4000$, $A^{+}=12$, and
$\kappa^{+}_{x}=0.007$.

For the cases with $\ell^{*}_{0.01}\lesssim 30$, the streamwise pre-multiplied
spectrograms $k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ (figures
8k,l,o) show the attenuation of $u^{\prime\prime 2}$ below
$y^{*}\simeq\ell^{*}_{0.01}$ (i.e. within the Stokes layer). For these cases,
increasing $\ell^{*}_{0.01}$ (hence strengthening the Stokes layer) attenuates
$u^{\prime\prime}$ over a wider range of wavelength $\lambda^{*}_{z}$ and
height $y^{*}$ (e.g. compare figure 8k with 8l). At the same time, the
energetic peak in $k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ is
shifted to a higher $y^{*}$ and a higher $\lambda^{*}_{z}$. This attenuation
is apparent in the visualisations of the instantaneous velocity fields of
$u^{\prime\prime}$ at $y^{*}=10$ for the cases with $\ell^{*}_{0.01}\lesssim
30$ (figures 9e,f). On the $w^{\prime\prime}$ fields, we overlay the spanwise
and phase-averaged $\tilde{w}^{*}$ (solid black curves) as a measure of the
Stokes motion. As $\ell^{*}_{0.01}$ increases from $20$ at $\omega^{+}=-0.05$
(figure 9e) to $30$ at $\omega^{+}=0$ (figure 9f), the Stokes motion becomes
stronger and the energy level in the $u^{\prime\prime}$ 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
$\ell^{*}_{0.01}$. Overall, for
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>30$, the Stokes layer dampens the level of turbulence within
$y^{*}\lesssim{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color<EMAIL_ADDRESS>hence acting favourably towards increasing $DR$.

For the cases with $\ell^{*}_{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^{\prime\prime}$ and $w^{\prime\prime}$ fields for the cases with
$\ell^{*}_{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^{\prime\prime}$
fields. This meandering is also evident in the spectrograms. For instance, for
$\omega^{+}=+0.05$ with $\ell^{*}_{0.01}\simeq 90$ (figures 8m,r), the
meandering flow structures at $y^{*}\simeq 10$ (visualised in figure 9g)
manifest as an energetic peak in the
$k^{*}_{z}\phi^{*}_{w^{\prime\prime}w^{\prime\prime}}$ spectrogram (figure 8r)
at $(\lambda^{*}_{z},y^{*})\simeq(100,10)$; this peak coincides with the peak
in the $k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ spectrogram
(figure 8m). Touber & Leschziner (2012) relate the attenuation of
$\left<u^{\prime\prime 2}\right>_{xzt}$ and amplification
$\left<w^{\prime\prime 2}\right>_{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^{\prime\prime}$ to $w^{\prime\prime}$, and the anisotropy between
$u^{\prime\prime}$ and $w^{\prime\prime}$ is reduced. Considering the flow
visualisation at $y^{*}=10$ for $\omega^{+}=+0.05$ with $\ell^{*}_{S}\simeq
90$ (figure 9g), we see a strong resemblance between the $u^{\prime\prime}$
and $w^{\prime\prime}$ fields in terms of the energy level and structure,
which support the reduction in anisotropy.

A noticeable difference between the cases with $\ell^{*}_{0.01}\lesssim 30$
and those with $\ell^{*}_{0.01}>30$ is the wall distance of the maximum
turbulence activity given by the location of the energetic peaks in
$k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ and
$k^{*}_{z}\phi^{*}_{w^{\prime\prime}w^{\prime\prime}}$. For the cases with
$\ell^{*}_{0.01}\lesssim 30$ (figures 8k,l,o), the energetic peak in
$k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ is lifted away from the
wall to a $y^{*}$ distance that coincides with $\ell^{*}_{0.01}$. However, for
the cases with $\ell^{*}_{0.01}>30$, the energetic peaks in
$k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ (figures 8m,n) and
$k^{*}_{z}\phi^{*}_{w^{\prime\prime}w^{\prime\prime}}$ (figures 8r,s) instead
reside near the wall at $y^{*}\simeq 10$, well below $\ell^{*}_{0.01}\simeq
90$. It appears that when $\ell^{*}_{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 $\ell^{*}_{0.01}\lesssim 30$ where the
turbulence is damped within the Stokes layer and the cycle of turbulence
generation is lifted away from the wall.

Overall, through flow visualisations and spectrograms we could explain the
physics behind the trends in $DR$ versus $\ell^{*}_{0.01}$ (figure 7). When
$\ell^{*}_{0.01}\lesssim 30$, turbulence is damped within
$y^{*}\lesssim\ell^{*}_{0.01}$. The level of damping increases by increasing
$\ell^{*}_{0.01}$. As a result, $DR$ increases by increasing
$\ell^{*}_{0.01}$, with the maximum $DR$ attained when $\ell^{*}_{0.01}\simeq
30$. However, when $\ell^{*}_{0.01}>30$, the Stokes layer becomes excessively
strong. In this situation, a near-wall cycle of turbulence is generated at
$y^{*}\simeq 10$ that meanders following the Stokes motion. As a result, $DR$
drops by increasing $\ell^{*}_{0.01}$.

Figure 8: Profiles of the turbulence statistics (a–j) and pre-multiplied
spectrograms (k–t) at $Re_{\tau}=4000$ for the non-actuated case and the
actuated cases with $A^{+}=12,\kappa^{+}_{x}=0.007$ and different values of
$\omega^{+}$ (same cases as in figure 6); $\omega^{+}=-0.05$ (a,f,k,p),
$\omega^{+}=0$ (b,g,l,q), $\omega^{+}=+0.05$ (c,h,m,r), $\omega^{+}=+0.10$
(d,i,n,s) and $\omega^{+}=+0.20$ (e,j,o,t). (a–e) profiles of turbulent
stresses for the streamwise and spanwise velocity $\left<u^{\prime\prime
2}\right>_{xzt},\left<w^{\prime\prime 2}\right>_{xzt}$. (f–j) profiles of
turbulent kinetic energy
$\left<\mathcal{K}\right>_{xzt}=\left(\left<u^{\prime\prime
2}\right>_{xzt}+\left<v^{\prime\prime 2}\right>_{xzt}+\left<w^{\prime\prime
2}\right>_{xzt}\right)/2$. Throughout (a–j), black lines with symbols
correspond to the non-actuated case; lines with no symbol correspond to the
actuated case scaled by the actuated $u_{\tau}$ (solid line) and the non-
actuated $u_{\tau_{o}}$ (dashed line). Pre-multiplied spectrograms for the
turbulent part of the streamwise velocity
$k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ (k–o) and spanwise
velocity $k^{*}_{z}\phi^{*}_{w^{\prime\prime}w^{\prime\prime}}$ (p–t). The
filled contours correspond to the actuated cases and the line contours
correspond to the non-actuated case. The contour lines for
$k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}$ (k–o) change from $0.6$
to $4.8$ with an increment of $0.6$, and for
$k^{*}_{z}\phi^{*}_{w^{\prime\prime}w^{\prime\prime}}$ (p–t) change from $0.2$
to $1.8$ with an increment of $0.2$. We locate $\ell^{*}_{0.01}$ with a bullet
point.

$y^{*}=50$$y^{*}=10$

Figure 9: Profiles of the turbulence statistics (a–d) and visualisation of the
near-wall turbulence (e–n) at $Re_{\tau}=4000$ for the non-actuated case (i,n)
and the actuated cases (e–h, j–m) with $A^{+}=12,\kappa^{+}_{x}=0.007$ and
different vales of $\omega^{+}$ (same cases as in figures 6 and 8);
$\omega^{+}=-0.05$ (a,e,j), $\omega^{+}=0$ (b,f,k), $\omega^{+}=+0.05$ (c,g,l)
and $\omega^{+}=+0.10$ (d,h,m). (a–d) profiles of turbulent stresses for the
streamwise and spanwise velocity components $\left<u^{\prime\prime
2}\right>^{*}_{xzt},\left<w^{\prime\prime 2}\right>^{*}_{xzt}$ (same profiles
as in figures 8a–d); lines correspond to the actuated cases and lines with
symbols correspond to the non-actuated case. Visualisation of the near-wall
turbulence at (e–i) $y^{*}=10$ and (j–n) $y^{*}=50$ (located with vertical
dashed-dotted lines in a–d). In each of (e–n), the upper field shows the
turbulent streamwise velocity $u^{\prime\prime*}$ and the lower field shows
the turbulent spanwise velocity $w^{\prime\prime*}$. We overlay the spanwise
and phase-averaged spanwise velocity $\tilde{w}^{*}$ (as solid curves) onto
the $w^{\prime\prime*}$ field.

###  3.6 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 $NPS$:

$NPS=\frac{P^{+}_{0}-\left(P^{+}+P^{+}_{\mathrm{in}}\right)}{P^{+}_{0}}=DR-\frac{P^{+}_{\mathrm{in}}}{P^{+}_{0}}$
(3.4)

where $P^{+}=\left(1-DR\right)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^{+}_{\mathrm{in}}$ is the input power required to
oscillate the wall actuation mechanism (1) while neglecting any mechanical
losses. A positive $NPS$ 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^{+}_{\mathrm{in}}$. In Part 2
(Chandran et al., 2022), we use this theory to estimate $NPS$ 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.

$P^{+}_{\mathrm{in}}=\frac{1}{T^{+}_{\mathrm{avg}}L^{+}_{x}L^{+}_{z}}\int_{t^{+}}^{t^{+}+T^{+}_{\mathrm{avg}}}\int_{0}^{L^{+}_{x}}\int_{0}^{L^{+}_{z}}w^{+}_{s}\left(\frac{\partial
w^{+}}{\partial y^{+}}\Big{|}_{y^{+}=0}\right)dx^{+}dz^{+}dt^{+}$ (3.2)

where all the quantities are normalised by $\nu$ and the non-actuated
$u_{\tau_{o}}$ (hence superscripted with a cross symbol). In (3.2),
$T_{\mathrm{avg}}$ is the averaging time, $w_{s}$ is the instantaneous wall
velocity (1) and $\partial w^{+}/\partial 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^{+}_{\mathrm{in}}/P^{+}_{0}$ as
computed over our parameter space of $(\omega^{+},\kappa^{+}_{x})$ at
$Re_{\tau}=4000$. The map is much more symmetric about $\omega^{+}=0$ compared
to $DR$ (figure 3). We also see that substantially more power is required at
higher actuation frequencies. For example, $P^{+}_{\mathrm{in}}/P^{+}_{0}\%$
can reach up to $100\%$ when $\omega^{+}\simeq\pm 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^{+}_{\mathrm{in}}/P^{+}_{0}$ decreases to
about $30\%-35\%$.

Figure 10: (a) Map of $P^{+}_{\mathrm{in}}$ at $Re_{\tau}=4000$. The filled
contour and the line contours show the same quantity. (b) Filled contour is
the difference in calculation of $P^{+}_{\mathrm{in}}$ at $Re_{\tau}=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 $\ell^{*}_{0.01}$ (same as in figure 7a).

We can use (3.2) only if we have an estimate for $\partial w^{+}/\partial
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,

$\displaystyle{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}=}w^{+}(x^{+},y^{+},t^{+})=$
(3.3)
$\displaystyle{\color[rgb]{1,1,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,1,1}\pgfsys@color@gray@stroke{1}\pgfsys@color@gray@fill{1}=}A^{+}\mathcal{R}\left\\{Ce^{i(\kappa^{+}_{x}x^{+}-\omega^{+}t^{+})}\mbox{Ai}\left[e^{\pi
i/6}\left(\kappa^{+}_{x}[1-DR]\right)^{1/3}\left(y^{+}-\frac{\omega^{+}}{\kappa^{+}_{x}[1-DR]}-\frac{i\kappa^{+}_{x}}{1-DR}\right)\right]\right\\}$

where
$C=\left\\{\mbox{Ai}\left[ie^{i\pi/3}\left(\kappa^{+}_{x}\left[1-DR\right]\right)^{1/3}(\omega^{+}/\kappa^{+}_{x}+i\kappa^{+}_{x})/[1-DR]\right]\right\\}^{-1}$,
Ai is the Airy function of the first kind, and $\mathcal{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., $\tilde{w}$ is the same in the laminar and turbulent flow, and 2) the
turbulent spanwise velocity is negligible near the wall, i.e.,
$w^{\prime\prime}\simeq 0$.

We now compare the results for $P^{+}_{\mathrm{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_{\tau}=4000$
($P^{+}_{\mathrm{in}}/P^{+}_{0}\%_{\mathrm{LES},Re_{\tau}=4000}$) and that
estimated using GSL theory at the same $Re_{\tau}$
($P^{+}_{\mathrm{in}}/P^{+}_{0}\%_{\mathrm{GSL},Re_{\tau}=4000}$). Overall,
the differences are small, especially in the range $\omega^{+}<0$ where they
differ less than $3\%$. Only in region II with $\omega^{+}>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 $\ell^{*}_{0.01}$ (contour
lines) indicates that the region where the power differences are significant
coincides with the region where $\ell^{*}_{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.

This behaviour is further substantiated by figure 11, which compares the
phase-averaged (harmonic) Reynolds stress profiles
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ between LES (solid lines) and its
laminar solution (3.3) from GSL theory (dashed lines with symbols). We see
that in region I with $\omega^{+}\leq 0$ (figure 11a), the agreement is
reasonably good; we obtain better agreement with $\omega^{+}=-0.05$
($\ell^{*}_{0.01}\simeq 20$) than with $\omega^{+}=0$ ($\ell^{*}_{0.01}\simeq
30$). However, in region II (figure 11b), when $\omega^{+}=+0.05$
($\ell^{*}_{0.01}\simeq 90$) and $\omega^{+}=+0.10$ ($\ell^{*}_{0.01}\simeq
80$), we observe significant departures between LES and the GSL theory. For
instance, for $\omega^{+}=+0.05$ at $y^{*}\simeq 20$,
$\left<\tilde{w}^{2}\right>^{*}_{xz}$ from LES is $2.2$ but from GSL is $0.2$.
This is a significant difference considering that the background turbulent
stress $\left<w^{\prime\prime 2}\right>^{*}_{xzt}\sim\mathcal{O}(1)$. In
region III with $\omega^{+}=+0.20$ where $\ell^{*}_{0.01}\simeq 17$, we see a
return of the good agreement between LES and GSL theory.

Our observations regarding the differences between $P_{\mathrm{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 drag-increasing
range. They explain this behaviour through the timescale
$\mathcal{T}^{+}\equiv 2\pi/(\omega^{+}-\kappa^{+}_{x}\mathcal{U}^{+}_{w})$,
which represents the period of oscillation as observed by the near-wall eddies
with the convective speed $\mathcal{U}^{+}_{w}\simeq 10$. As discussed in §
3.1, in the drag-increasing range
$\omega^{+}/\kappa^{+}_{x}\rightarrow\mathcal{U}^{+}_{w}$ leading to
$\mathcal{T}^{+}\rightarrow\infty$. 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 $\left<\tilde{w}^{2}\right>^{*}_{xt}$ profiles from the GSL solution
(figure 11b) supports the activation of these terms.

Figure 11: Comparison of the phase-averaged (harmonic) Reynolds stress
profiles $\left<\tilde{w}^{2}\right>^{*}_{xt}$ between LES (solid lines) and
the laminar solution from the generalised Stokes layer (GSL) theory (dashed
lines with symbols) ($A^{+}=12,\ \kappa^{+}_{x}=0.007$ at $Re_{\tau}=4000$).
(a) $\omega^{+}=-0.05,0$. (b) $\omega^{+}=+0.05,+0.10,+0.20$. The insets plot
$DR$ for the selected cases.

We can now explore the net power saving ($NPS$). Figure 12(a) demonstrates
that for our considered parameter space $NPS$ is mostly negative. The highest
(best) $NPS$ is $0.5\%$ at $(\omega^{+},\kappa^{+}_{x})\simeq(0.05,0.012)$. In
figure 12(b), we plot the map of the difference between $NPS$ from LES at
$Re_{\tau}=4000$ and its counterpart at $Re_{\tau}=951$. If this difference is
positive, $NPS$ increases with $Re_{\tau}$. Over a large portion of our
parameter space, the difference is negative, i.e., $NPS$ becomes more negative
with $Re_{\tau}$. However, for a small portion of region I with $\omega^{+}<0$
and $\kappa^{+}_{x}\lesssim 0.0025$, the difference is positive. One
experimental case reported by Marusic et al. (2021) falls into this region,
with $\omega^{+}=-0.044$ and $\kappa^{+}_{x}=0.0014$ (see their figures 3a,b).
The $NPS$ of this case was negative, but it increased with $Re_{\tau}$ in
accordance with our analysis.

Quadrio et al. (2009), similar to figure 12(a), generate a map of $NPS$ for
their travelling wave study at $Re_{\tau}=200$ (their figure 5). They report
$NPS>0$ within the range
$0\lesssim\omega^{+}\lesssim+0.05,0.002\lesssim\kappa^{+}_{x}\lesssim 0.025$.
This range coincides with the range where $DR\gtrsim 40\%$. Gatti & Quadrio
(2013) generate a similar map at $Re_{\tau}=1000$ (their figure 9). They also
observe $NPS>0$ within the same range of $(\omega^{+},\kappa^{+}_{x})$.
However, the level of $NPS>0$ is lower at $Re_{\tau}=1000$ compared to
$Re_{\tau}=200$. Considering (3.4), we speculate that the decrease in $NPS>0$
from $Re_{\tau}=200$ to $1000$ is due to the decrease in $DR$. We observe a
similar trend in figure 12(b). Within the range of
$(\omega^{+},\kappa^{+}_{x})$ where $DR$ is maximum (blue dashed-dotted line),
$NPS$ decreases by increasing $Re_{\tau}$ from $1000$ to $4000$.

Figure 12: (a) Net power saving $NPS$ for LES at $Re_{\tau}=4000$. (b)
Difference in $NPS$ between LES at $Re_{\tau}=4000$ and LES at
$Re_{\tau}=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
$\nu$ and the non-actuated $u_{\tau_{o}}$. ( ) and ( ) locate the local
maximum and the local minimum $DR$, respectively.

Overall, for our considered parameter space, $NPS$ 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 $NPS$ that
actually increase with Reynolds number.

##  4 Conclusions

Turbulent drag reduction was considered using spanwise wall oscillation based
on streamwise travelling waves at friction Reynolds numbers $Re_{\tau}=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\leq\kappa^{+}_{x}\leq 0.02$ and $-0.2\leq\omega^{+}\leq+0.2$, covering
upstream ($\omega^{+}<0$) and downstream ($\omega^{+}>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 $\ell_{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 $\ell_{0.01}$, hence
the Stokes layer protrusion, the $DR$ map over the parameter space of
$(\omega^{+},\kappa^{+}_{x})$ can be categorised into two regions. When
$\ell_{0.01}$ is less than $30$ viscous units, increasing $\ell_{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 $\ell_{0.01}$
and lifts the cycle of turbulence generation away from the wall. Increasing
$\ell_{0.01}$ in this regime, further attenuates the turbulence and leads to
higher $DR$. When $\ell_{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 $\omega^{+}$ and
$\kappa^{+}_{x}$ at $Re_{\tau}=4000$ the net-power savings ($NPS$) 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 $NPS$ 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_{\tau}$ up to $\mathcal{O}(10^{4})$.

To afford the parametric study conducted here, we employed a reduced
simulation domain size. This setup was found to be suitable for the ISA
pathway, especially for studying $DR$, Stokes layer dynamics and the near-wall
turbulence. However, the present configuration cannot resolve the outer scale
eddies, which become important in the OSA pathway. This aspect will also be
investigated in Part 2, where the inner and outer scale eddies are captured
through experimental techniques and for higher Reynolds numbers.

Acknowledgements
The research was funded through the Deep Science Fund of Intellectual
Ventures. We acknowledge Dr. Daniel Chung for providing insightful comments,
and sharing his DNS solver and computing resources during the early stages of
this work. Computing resources were provided through the Spartan High-
Performance Computing service at The University of Melbourne, ARCHER2 UK
National Supercomputing Service (https://www.archer2.ac.uk), the Pawsey
Supercomputing Centre with funding from the Australian Government and the
Government of Western Australia, and the National Computing Infrastructure
(NCI), which is supported by the Australian Government.

Declaration of interests. The authors report no conflict of interest.

##  Appendix A Validation of LES and grid resolution study for drag reduction

We perform several validation studies for LES. In this appendix, we focus on
the accuracy of the dynamic Smagorinsky subgrid-scale model (Germano et al.,
1991) in predicting $DR$. We also assess the proper grid resolution for
predicting $DR$. In Appendix B, we perform a grid resolution study for the
Reynolds stresses and their spectra.

Our first validation study is summarised in figure 13. We compare $DR$ between
LES, the experimental data of Marusic et al. (2021) and the DNS data of Gatti
& Quadrio (2016). All sets of data have matched actuation parameters
$A^{+}\simeq 12,\kappa^{+}_{x}\simeq 0.0014,\omega^{+}\simeq-0.044$. For the
LES cases (table 2 above the separating line), we change $Re_{\tau}$ from
$951$ to $6000$. The LES cases at $Re_{\tau}=951$ and $6000$ are comparable
with the DNS of Gatti & Quadrio (2016) and experiments of Marusic et al.
(2021), respectively. All the LES cases have the viscous-scaled grid size
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 60\times 30$. We use the full domain
size (figure 1c) at $Re_{\tau}=951,2000$ and $4000$ (red bullet), and the
medium domain size (figure 1a) at $Re_{\tau}=4000$ and $6000$ (black circle).

$Re_{\tau}$ | $\kappa^{+}_{x}$ | $\omega^{+}$ | $L_{x},L_{z}$ | $N_{x},N_{y},N_{z}$ | $\Delta^{+}_{x},\Delta^{+}_{z}$ | Sym.
---|---|---|---|---|---|---
$951$ | $0.0014$ | $-0.044$ | $9.44h,3.14h$ | $144,48,96$ | $62,33$ | $\bullet$
$2000$ | $0.0014$ | $-0.044$ | $6.73h,3.14h$ | $224,96,192$ | $60,33$ | $\bullet$
$4000$ | $0.0014$ | $-0.044$ | $6.73h,3.14h$ | $448,192,384$ | $60,33$ | $\bullet$
$4000$ | $0.0014$ | $-0.044$ | $2.24h,0.63h$ | $144,192,72$ | $62,35$ | $\circ$
$6000$ | $0.0014$ | $-0.044$ | $2.24h,0.63h$ | $216,288,108$ | $62,34$ | $\circ$
$951$ | $0.0347$ | $-0.20,-0.05,-0.01,+0.06,+0.12,+0.28$ | $6.86h,3.14h$ | $108,48,96$ | $60,31$ | $\blacklozenge$
$951$ | $0.0347$ | $-0.20,-0.05,-0.01,+0.06,+0.12,+0.28$ | $6.86h,3.14h$ | $288,48,96$ | $23,31$ | $\blacklozenge$
$951$ | $0.0347$ | $-0.20,-0.05,-0.01,+0.06,+0.12,+0.28$ | $6.86h,3.14h$ | $384,64,128$ | $17,23$ | $\blacklozenge$
$951$ | $0.0208$ | $-0.10,0,+0.05,+0.07,+0.15,+0.20$ | $6.35h,3.14h$ | $256,48,96$ | $24,31$ | $\blacklozenge$
$951$ | $0.0208$ | $-0.10,0,+0.05,+0.07,+0.15,+0.20$ | $6.35h,3.14h$ | $384,64,128$ | $16,23$ | $\blacklozenge$
Table 2: Summary of the LES cases for validation. The cases above the
separating line have fixed actuation parameters
$A^{+},\kappa^{+}_{x},\omega^{+}$ and grid resolution
$\Delta^{+}_{x}\times\Delta^{+}_{z}$, but $Re_{\tau}$ changes from $951$
(first row) to $6000$ (fifth row). These cases are compared with the DNS of
Gatti & Quadrio (2016) and experiments of Marusic et al. (2021) at matched
actuation parameters and Reynolds number (figure 13). The cases below the
separating line have fixed $Re_{\tau}=951$ and $A^{+}=12$, but
$\kappa^{+}_{x},\omega^{+}$ and grid resolution change. Each row consists of
six cases with fixed $Re_{\tau},A^{+},\kappa^{+}_{x}$ and grid resolution, but
$\omega^{+}$ is different for each case. These cases are for validation
against the DNS of Gatti & Quadrio (2016) at selected actuation parameters
(figure 14). Figure 13: Comparison of $DR$ from the LES runs in table 2
($\bullet$, $\circ$) with the DNS of Gatti & Quadrio (2016) at $Re_{\tau}=951$
($\diamond$), and experiment of Marusic et al. (2021) at $Re_{\tau}=6000$
using hot-wire anemometry
(${\color[rgb]{0,1,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,1,0}\square}$)
and drag balance ($\triangle$ ). All the data points from different techniques
have matched actuation parameters
$A^{+}=12,\kappa^{+}_{x}=0.0014,\omega^{+}=-0.044$. For the LES, we use full
domain (figure 1c) at $Re_{\tau}=951,2000,4000$ ($\bullet$), and medium domain
(figure 1a) at $Re_{\tau}=4000,6000$ ($\circ$). We overlay GQ’s predictive
model for $DR$ ( ).

Considering figure 13, at $Re_{\tau}=951$ we obtain good agreement between LES
(red bullet) and DNS of Gatti & Quadrio (2016) (blue diamond), and at
$Re_{\tau}=6000$ we obtain good agreement between LES (black circle) and the
experimental data of Marusic et al. (2021) from the hot-wire anemometry (green
square) and drag balance (green triangle). These agreements support the
accuracy of the dynamic Smagorinsky model (Germano et al., 1991) for LES. At
$Re_{\tau}=4000$, we obtain less than $1\%$ difference between the LES case
with the medium domain (black circle) and the case with the full domain (red
bullet). This agreement supports the suitability of the medium domain size for
the actuation parameters considered here. We further demonstrate the accuracy
of the medium domain size in Appendix C. All the data points from DNS, LES and
experiments agree well with GQ’s predictive model for $DR$ (dashed-dotted
line). This agreement is because the actuation frequency $\omega^{+}=-0.044$
($T^{+}_{osc}\simeq 142$) falls into the inner-scaled actuation pathway
($T^{+}_{osc}<350$). As discussed in Marusic et al. (2021) and § 1, GQ’s model
performs accurately in this pathway.

Our second validation study is shown in figure 14. We compare the present LES
with DNS dataset of Gatti & Quadrio (2016) at matched $Re_{\tau}=951$ over a
range of actuation parameters within our parameter space of interest. We
compare at $A^{+}=12,\kappa^{+}_{x}=0.0347$ (figure 14a) and
$A^{+}=12,\kappa^{+}_{x}=0.0208$ (figure 14b) over the range
$-0.20\lesssim\omega^{+}\lesssim+0.28$. Table 2 lists the LES cases for this
validation study. For
$A^{+}=12,\kappa^{+}_{x}=0.0347,-0.20\leq\omega^{+}\leq+0.28$, we perform LES
with three grids $(\Delta^{+}_{x}\times\Delta^{+}_{z})=(60\times 31)$,
$(23\times 31)$, $(17\times 23)$. For
$A^{+}=12,\kappa^{+}_{x}=0.0208,-0.10\leq\omega^{+}\leq+0.20$, we perform LES
with two grids $(\Delta^{+}_{x}\times\Delta^{+}_{z})=(24\times 31)$,
$(16\times 23)$. Figure 14 shows that the LES grid
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 23\times 31$ (blue diamond) yields
good agreement with DNS for all the compared cases. Also, this LES grid yields
grid convergence. Further grid refinement to
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 16\times 23$ (green diamond) does
not significantly change $DR$. In our first validation study with the
experiments (figure 13), we employed the LES grid
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 60\times 31$. We also employed this
grid for our second validation study with
$A^{+}=12,\kappa^{+}_{x}=0.0347,-0.20\leq\omega^{+}\leq+0.28$ (red diamond in
table 2 and figure 15a). We observe that this grid performs accurately for the
upstream travelling wave ($\omega^{+}<0$ in figure 15a). This observation is
consistent with our first validation study with $\omega^{+}=-0.044$ (figure
13). However, for the downstream travelling wave ($\omega^{+}>0$), the LES
grid $\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 60\times 31$ (red diamond)
under-predicts $DR$. Further refinement to
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 23\times 31$ (blue diamond) improves
the prediction of $DR$ for all the values of $\omega^{+}$.

We conclude that with the viscous-scaled grid resolution of
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 23\times 31$ (blue diamond) we can
study $DR$ with high confidence. Therefore, we adopt this grid resolution to
study $DR$ (table 1).

Figure 14: Comparison between the DNS of Gatti & Quadrio (2016) ($\square$)
and the LES of the present study at different grid resolutions:
$\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 60\times 31$ ($\blacklozenge$),
$23\times 31$ ($\blacklozenge$), $16\times 23$ ($\blacklozenge$). Table 2
lists the simulation details for LES. Both DNS and LES cases are compared at
matched $Re_{\tau}=951$, $A^{+}=12$, and $\kappa^{+}_{x},\omega^{+}$. (a,b)
show the comparison at $\kappa^{+}_{x}=0.0347$ and $0.0208$, respectively. At
each value of $\kappa^{+}_{x}$, comparison is made at six values of
$\omega^{+}$ (listed in table 2).

##  Appendix B Grid resolution study for Reynolds stresses and spectra

Where the previous section determined the adequate grid resolution for
calculating the drag reduction, we conduct a similar analysis to assess the
proper grid spacing for resolving the Reynolds stresses and velocity spectra.
These are the quantities that we investigate to explain the flow physics (§
3).

case | $Re_{\tau}$ | $A^{+}$ | $\kappa^{+}_{x}$ | $\omega^{+}$ | $L_{x},L_{z}$ | $N_{x},N_{y},N_{z}$ | $\Delta^{+}_{x},\Delta^{+}_{z}$ | line
---|---|---|---|---|---|---|---|---
coarse LES | $590$ | $12$ | $0.0014$ | $-0.044$ | $7.6h,3.14h$ | $204,30,60$ | $22,31$ |
fine LES | $590$ | $12$ | $0.0014$ | $-0.044$ | $7.6h,3.14h$ | $306,45,90$ | $15,21$ |
DNS | $590$ | $12$ | $0.0014$ | $-0.044$ | $7.6h,3.14h$ | $608,240,480$ | $7,4$ |
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_{\tau}$, actuation parameters $(A^{+},\kappa^{+}_{x},\omega^{+})$ and
domain size $L_{x},L_{z}$, 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.

To evaluate the accuracy of LES for the Reynolds stresses and spectra, we
generate a DNS dataset $(\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 7\times 4)$
in a full-domain open-channel flow with wall actuation (table 3). To afford
the DNS, we consider $Re_{\tau}=590$ with the actuation parameters
$(A^{+},\kappa^{+}_{x},\omega^{+})=(12,0.0014,-0.044)$. We perform two LES
calculations that match the DNS case in terms of the domain size, $Re_{\tau}$
and actuation parameters, but have different grid resolutions (table 3). We
name the LES case with a coarser grid
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$) “coarse LES”, and the
case with a finer grid ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 15\times
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).

Figure 15: Comparison between the coarse LES
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$,  ), fine LES
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 21$,  ) and DNS
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 7\times 4$,  ). All cases have the
same $Re_{\tau}=590$ and actuation parameters
$(A^{+},\kappa^{+}_{x},\omega^{+})=(12,0.0014,-0.044)$, see table 3. The
comparison is based on (a) mean velocity profiles $U^{*}$ and $DR$, and (b)
profiles of the Reynolds stress by the phase-averaged spanwise velocity
$\left<\tilde{w}^{2}\right>^{*}_{xt}$. Figure 16: Comparison between the
coarse LES ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$,  ), fine
LES ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 21$,  ) and DNS
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 7\times 4$, filled contour). All
cases have matched $Re_{\tau}=590$ and actuation parameters
$(A^{+},\kappa^{+}_{x},\omega^{+})=(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}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}(\lambda^{*}_{z},y^{*})$,
and (e,f) two-dimensional pre-multiplied spectrograms of the turbulent part of
the streamwise velocity
$k^{*}_{x}k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}(\lambda^{*}_{x},\lambda^{*}_{z})$
at $y^{*}\simeq 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).

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
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ (figure 15b). We use
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ to calculate the protrusion height by
the Stokes layer (§ 3.4). Figure 15 shows that $U^{*},DR$ and
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ are predicted reasonably well with the
coarse LES grid ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$). We
also concluded in Appendix A that the coarse LES grid predicts $DR$ quite
well. Therefore, we employ the coarse LES grid
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 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 ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 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}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}(\lambda^{*}_{z},y^{*})$.
The coarse LES spectrogram (red contour lines in figure 16a) is highly
distorted for $\lambda^{*}_{z}\lesssim 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}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}\left(\lambda^{*}_{x},\lambda^{*}_{z}\right)$
at $y^{*}\simeq 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.

Figure 17: Comparison between the coarse LES
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$,  ), fine LES
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 21$,  ) and DNS
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 7\times 4$,  ) in terms of the
Reynolds stress profiles due to the turbulent part of the streamwise velocity
$\left<u^{\prime\prime 2}\right>^{*}_{xzt}$. All cases have the same
$Re_{\tau}=590$ and actuation parameters
$(A^{+},\kappa^{+}_{x},\omega^{+})=(12,0.0014,-0.044)$, see table 3. The LES
profiles ( ,  ) in (a) are obtained from the raw LES data, and in (b) are
obtained by integrating the dealiased spectrograms (figures 16c,d).

Refining the LES grid improves the spectrograms (figures 16b,d,f). In figure
16(b), we compare the one-dimensional spectrogram of the fine LES (blue
contour lines) with DNS (filled contour). The range of scales affected by the
aliasing error is narrowed to $\lambda^{*}_{z}\lesssim 50$. Attenuation of the
aliasing error by the grid refinement is also evident in the two-dimensional
spectrograms (compare figure 16e with 16f). Further improvement is achieved by
removing the aliased scales (dealiasing). We perform dealiasing through the
$k^{*}_{x}k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}\left(\lambda^{*}_{x},\lambda^{*}_{z}\right)$
spectrogram at each $y^{*}$. We can explain the dealiasing process through
figures 16(e,f). At each $\lambda^{*}_{x}$, if aliasing error occurs, a local
minimum appears in
$k^{*}_{x}k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}\left(\lambda^{*}_{x},\lambda^{*}_{z}\right)$.
In figures 16(e,f), we mark the local minima at all values of
$\lambda^{*}_{x}$ and connect them together with a grey line. Thus, the grey
line separates the healthy scales from the aliased scales. For dealiasing, we
remove the aliased scales below the grey line. After dealiasing
$k^{*}_{x}k^{*}_{z}\phi^{*}_{u^{\prime\prime}u^{\prime\prime}}\left(\lambda^{*}_{x},\lambda^{*}_{z}\right)$
at each $y^{*}$, we integrate it to reconstruct the dealiased one-dimensional
spectrograms (figures 16c,d). Accordingly, we integrate the dealiased one-
dimensional spectrograms to reconstruct the dealiased Reynolds stress profiles
$\left<u^{\prime\prime 2}\right>^{*}_{xzt}$ (figure 17b). Comparing the
original spectrograms from the raw LES data (figures 16a,b) with the dealiased
spectrograms (figures 16c,d), highlights the improvement due to dealiasing.
Similarly, comparing the original $\left<u^{\prime\prime 2}\right>^{*}_{xzt}$
profiles from the raw LES data (figure 17a) with the dealiased
$\left<u^{\prime\prime 2}\right>^{*}_{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
($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 22\times 31$) is suitable for
studying $DR$, mean velocity profiles $U^{*}$ and
$\left<\tilde{w}^{2}\right>^{*}_{xt}$ (for the Stokes layer dynamics). The
fine LES resolution ($\Delta^{+}_{x}\times\Delta^{+}_{z}\simeq 14\times 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_{\tau}=4000$
with $A^{+}=12,\kappa^{+}_{x}=0.0014,\omega^{+}=-0.044$. Here, we further
study the domain size effect for some of our production cases at
$Re_{\tau}=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
$(\kappa^{+}_{x},\omega^{+})=(0.021,-0.1),(0.021,+0.1),(0.007,+0.05)$. The
cases with $\kappa^{+}_{x}=0.021$ fall at the upper bound of our range of
interest for $\kappa^{+}_{x}$, and the case with $\kappa^{+}_{x}=0.007$ falls
within this range. Also, we consider cases with upstream travelling waves
($\omega^{+}<0$) and downstream travelling waves ($\omega^{+}>0$). We
deliberately choose the case with $(\kappa^{+}_{x},\omega^{+})=(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\times 0.6h,y^{+}_{\mathrm{res}}\simeq 1000$) and the large domain size
(figure 1b, $4.0h\times 1.2h,y^{+}_{\mathrm{res}}\simeq 2000$).

domain | $Re_{\tau}$ | $y^{+}_{\mathrm{res}}$ | $\kappa^{+}_{x}$ | $\omega^{+}$ | $L_{x}/h,,L_{z}/h$ | $N_{x},N_{y},N_{z}$ | $\Delta^{+}_{x},\Delta^{+}_{z}$ | $DR$
---|---|---|---|---|---|---|---|---
medium | $4000$ | $1000$ | $0.021$ | $-0.1$ | $2.04,0.63$ | $384,192,80$ | $21,31$ | $17.1$
large | $4000$ | $2000$ | $0.021$ | $-0.1$ | $4.08,1.25$ | $768,192,160$ | $21,31$ | $18.7$
medium | $4000$ | $1000$ | $0.021$ | $+0.1$ | $2.04,0.63$ | $384,192,80$ | $21,31$ | $21.1$
large | $4000$ | $2000$ | $0.021$ | $+0.1$ | $4.08,1.25$ | $768,192,160$ | $21,31$ | $21.6$
medium | $4000$ | $1000$ | $0.007$ | $+0.05$ | $2.04,0.63$ | $384,192,80$ | $21,31$ | $11.3$
large | $4000$ | $2000$ | $0.007$ | $+0.05$ | $4.08,1.25$ | $768,192,160$ | $21,31$ | $12.7$
Table 4: Summary of the LES cases for the domain size study. For all cases,
$A^{+}=12$. We consider three cases with
$(\kappa^{+}_{x},\omega^{+})=(0.021,-0.1),(0.021,+0.1),(0.007,+0.05)$. For
each case, we perform a medium domain simulation (figure 1a, $L_{x}\times
L_{z}\simeq 2.0h\times 0.6h,y^{+}_{\mathrm{res}}\simeq 1000$), and a large
domain simulation (figure 1b, $L_{x}\times L_{z}\simeq 4.0h\times
1.2h,y^{+}_{\mathrm{res}}\simeq 2000$). Figure 18: Comparison of the mean
velocity profiles $U^{*}$ between the medium domain simulation ( $L_{x}\times
L_{z}\simeq 2.0h\times 0.6h,y^{+}_{\mathrm{res}}\simeq 1000$) and the large
domain simulation ( $L_{x}\times L_{z}\simeq 4.0h\times
1.2h,y^{+}_{\mathrm{res}}\simeq 2000$) for two actuated cases from table 4:
(a) $Re_{\tau}=4000,A^{+}=12,\kappa^{+}_{x}=0.021,\omega^{+}=+0.1$, and (b)
$Re_{\tau}=4000,A^{+}=12,\kappa^{+}_{x}=0.007,\omega^{+}=+0.05$. In both
(a,b), we also plot the non-actuated case at $Re_{\tau}=4000$ with the medium
domain size ( ). The profiles are presented in viscous units (scaled by their
actual $u_{\tau}$ and kinematic viscosity $\nu$). The bullet points
($\bullet$, $\bullet$, $\bullet$) locate the resolved height
$y^{*}_{\mathrm{res}}$. The profiles ( ,  ,  ) beyond $y^{*}_{\mathrm{res}}$
are reconstructed using the composite profile of Nagib & Chauhan (2008) (see §
2.3). The insets plot the velocity difference $\Delta
U^{*}=U^{*}_{\mathrm{act}}-U^{*}_{\mathrm{non{\text{-}}act}}$.

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^{*}_{\mathrm{res}}$ using
Nagib & Chauhan (2008)’s composite profile, indicated with dashed line in
figure 18. Then, we obtain $C_{f}\equiv 2/{U^{*}_{b}}^{2}$ by integrating the
resolved portion of the $U^{*}$ profile up to $y^{*}_{\mathrm{res}}$ (solid
line in figure 18) and its reconstructed portion beyond $y^{*}_{\mathrm{res}}$
(dashed line in figure 18). Therefore, for the medium domain, $C_{f}$ is
obtained by integrating the resolved $U^{*}$ profile up to
$y^{*}_{\mathrm{res}}\simeq 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^{*}_{\mathrm{res}}\simeq 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^{*}_{\mathrm{res}}\simeq 1000$). Beyond $y^{*}_{\mathrm{res}}\simeq 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 $\Delta U^{*}$ between the medium domain size (red solid
line) and the large domain size (blue solid line) for two cases from table 4;
$\kappa^{+}_{x}=0.021,\omega^{+}=+0.1$ (figure 18a) and
$\kappa^{+}_{x}=0.007,\omega^{+}=+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 $\Delta U^{*}$
profiles (the insets). For both cases, the logarithmic $U^{*}$ profile appears
by $y^{*}\simeq 200$. This allows to use the composite profile beyond
$y^{*}\simeq 200$. Overall, we conclude that the medium domain size
($L_{x}\times L_{z}\simeq 2.0h\times 0.6h$, figure 1a) is suitable for our
production simulations at $Re_{\tau}=4000$ (table 1).

## References

  * Abbassi et al. (2017) Abbassi, M. R., Baars, W. J., Hutchins, N. & Marusic, I. 2017 Skin-friction drag reduction in a high-Reynolds-number turbulent boundary layer via real-time control of large-scale structures. Int. J. Heat Fluid Flow 67, 30–41.
  * Aider et al. (2010) Aider, J. L., Beaudoin, J. F. & Wesfreid, J. E. 2010 Drag and lift reduction of a 3D bluff-body using active vortex generators. Exp. Fluids 48, 771–789.
  * Auteri et al. (2010) Auteri, F., Baron, A., Belan, M., Campanardi, G. & Quadrio, M. 2010 Experimental assessment of drag reduction by traveling waves in a turbulent pipe flow. Phys. Fluids 22, 115103\.
  * Baidya et al. (2021) Baidya, R., Philip, J., Hutchins, N., Monty, J. P. & Marusic, I. 2021 Spanwise velocity statistics in high-Reynolds-number turbulent boundary layers. J. Fluid Mech. 913, A35.
  * Banyassady & Piomelli (2014) Banyassady, R. & Piomelli, U. 2014 Turbulent plane wall jets over smooth and rough surfaces. J. Turbul. 15, 186–207.
  * Baron & Quadrio (1995) Baron, A. & Quadrio, M. 1995 Turbulent drag reduction by spanwise wall oscillations. Appl. Sci. Res. 55, 311–326.
  * Batchelor (2000) Batchelor, G. K. 2000 An Introduction to Fluid Dynamics. Cambridge university press.
  * Bird et al. (2018) Bird, J., Santer, M. & Morrison, J. F. 2018 Experimental control of turbulent boundary layers with in-plane travelling waves. Flow Turbul. Combust. 100, 1015–1035.
  * Chan et al. (2015) Chan, L., MacDonald, M., Chung, D., Hutchins, N. & Ooi, A. 2015 A systematic investigation of roughness height and wavelength in turbulent pipe flow in the transitionally rough regime. J. Fluid Mech. 771, 743–777.
  * Chandran et al. (2020) Chandran, D., Monty, J.P. & Marusic, I. 2020 Spectral-scaling-based extension to the attached eddy model of wall turbulence. Physical Review Fluids 5 (10), 104606.
  * Chandran et al. (2022) Chandran, D., Zampiron, A., Rouhi, A., Fu, M. K, , Wine, D., Holloway, B., Smits, A. J. & Marusic, I. 2022 Turbulent drag reduction by spanwise wall forcing. Part 2: High-Reynolds-number experiments. J. Fluid Mech. (under consideration) .
  * Chin et al. (2014) Chin, C., Monty, J. P. & Ooi, A. 2014 Reynolds number effects in dns of pipe flow and comparison with channels and boundary layers. Int. J. Heat Fluid Flow 45, 33–40.
  * Choi & Clayton (2001) Choi, K. S. & Clayton, B. R. 2001 The mechanism of turbulent drag reduction with wall oscillation. Int. J. Heat Fluid Flow 22, 1–9.
  * Choi et al. (1998) Choi, K. S., DeBisschop, J. R. & Clayton, B. R. 1998 Turbulent boundary-layer control by means of spanwise-wall oscillation. AIAA J 36, 1157–1163.
  * Choi & Graham (1998) Choi, K. S. & Graham, M. 1998 Drag reduction of turbulent pipe flows by circular-wall oscillation. Phys. Fluids 10, 7–9.
  * Chung et al. (2015) Chung, D., Chan, L., MacDonald, M., Hutchins, N. & Ooi, A. 2015 A fast direct numerical simulation method for characterising hydraulic roughness. J. Fluid Mech. 773, 418–431.
  * Chung et al. (2014) Chung, D., Monty, J. P. & Ooi, A. 2014 An idealised assessment of Townsend’s outer-layer similarity hypothesis for wall turbulence. J. Fluid Mech. 742.
  * Dean & Bhushan (2010) Dean, B. & Bhushan, B. 2010 Shark-skin surfaces for fluid-drag reduction in turbulent flow: a review. Philo. Trans. R. Soc. A 368, 4775–4806.
  * Dean (1978) Dean, R. B. 1978 Reynolds number dependence of skin friction and other bulk flow variables in two-dimensional rectangular duct flow. J. Fluids Eng. 100, 215–223.
  * Di Cicca et al. (2002) Di Cicca, G. M., Iuso, G., Spazzini, P. G. & Onorato, M. 2002 Particle image velocimetry investigation of a turbulent boundary layer manipulated by spanwise wall oscillations. J. Fluid Mech. 467, 41–56.
  * Endrikat et al. (2021a) Endrikat, S., Modesti, D., García-Mayoral, R., Hutchins, N. & Chung, D. 2021a Influence of riblet shapes on the occurrence of Kelvin–Helmholtz rollers. J. Fluid Mech. 913, A37.
  * Endrikat et al. (2021b) Endrikat, S., Modesti, D., MacDonald, M., García-Mayoral, R., Hutchins, N. & Chung, D. 2021b Direct numerical simulations of turbulent flow over various riblet shapes in minimal-span channels. Flow Turbul. Combust. 107, 1–29.
  * Endrikat et al. (2022) Endrikat, S., Newton, R., Modesti, D., García-Mayoral, R., Hutchins, N. & Chung, D. 2022 Reorganisation of turbulence by large and spanwise-varying riblets. J. Fluid Mech. 952, A27.
  * García-Mayoral & Jiménez (2011) García-Mayoral, R. & Jiménez, J. 2011 Drag reduction by riblets. Phil. Trans. R. Soc. A 369 (1940), 1412–1427.
  * García-Mayoral & Jiménez (2012) García-Mayoral, R. & Jiménez, J. 2012 Scaling of turbulent structures in riblet channels up to ${R}e_{\tau}=550$. Phys. Fluids 24, 105101.
  * Gatti & Quadrio (2012) Gatti, D. & Quadrio, M. 2012 What happens to turbulent skin-friction drag reduction at high Re? arXiv preprint arXiv:1212.6052 .
  * Gatti & Quadrio (2013) Gatti, D. & Quadrio, M. 2013 Performance losses of drag-reducing spanwise forcing at moderate values of the Reynolds number. Phys. Fluids 25, 125109.
  * Gatti & Quadrio (2016) Gatti, D. & Quadrio, M. 2016 Reynolds-number dependence of turbulent skin-friction drag reduction induced by spanwise forcing. J. Fluid Mech. 802, 553–582.
  * Germano et al. (1991) Germano, M., Piomelli, U., Moin, P. & Cabot, W. H. 1991 A dynamic subgrid-scale eddy viscosity model. Phys. Fluids 3, 1760–1765.
  * Hasegawa & Kasagi (2011) Hasegawa, Y. & Kasagi, N. 2011 Dissimilar control of momentum and heat transfer in a fully developed turbulent channel flow. J. Fluid Mech. 683, 57–93.
  * Hurst et al. (2014) Hurst, E., Yang, Q. & Chung, Y. M. 2014 The effect of Reynolds number on turbulent drag reduction by streamwise travelling waves. J. Fluid Mech 759, 11.
  * Jiménez & Moin (1991) Jiménez, J. & Moin, P. 1991 The minimal flow unit in near-wall turbulence. J. Fluid Mech. 225, 213–240.
  * Jung et al. (1992) Jung, W.-J., Mangiavacchi, N. & Akhavan, R. 1992 Suppression of turbulence in wall-bounded flows by high-frequency spanwise oscillations. Phys. Fluids 4 (8), 1605–1607.
  * Kametani et al. (2015) Kametani, Y., Fukagata, K., Örlü, R. & Schlatter, P. 2015 Effect of uniform blowing/suction in a turbulent boundary layer at moderate Reynolds number. Int. J. Heat Fluid Flow 55, 132–142.
  * Koike et al. (2004) Koike, M., Nagayoshi, T. & Hamamoto, N. 2004 Research on aerodynamic drag reduction by vortex generators. Mitsubishi Motors Tech. Rev. 16, 11–16.
  * Kravchenko & Moin (1997) Kravchenko, A. G. & Moin, P. 1997 On the effect of numerical errors in large eddy simulations of turbulent flows. J. Comput.l Phys. 131, 310–322.
  * Lee & Moser (2015) Lee, M. & Moser, R. D. 2015 Direct numerical simulation of turbulent channel flow up to. J. Fluid Mech. 774, 395–415.
  * Lilly (1992) Lilly, D. K. 1992 A proposed modification of the Germano subgrid-scale closure method. Phys. Fluids 4, 633–635.
  * Lin (2002) Lin, J. C. 2002 Review of research on low-profile vortex generators to control boundary-layer separation. Prog. Aerosp. Sci. 38, 389–420.
  * Lozano-Durán & Jiménez (2014) Lozano-Durán, A. & Jiménez, J. 2014 Effect of the computational domain on direct simulations of turbulent channels up to ${R}e_{\tau}$= 4200. Phys. Fluids 26, 011702.
  * Luchini & Quadrio (2022) Luchini, P. & Quadrio, M. 2022 Wall turbulence and turbulent drag reduction. In 50+ Years of AIMETA, pp. 349–364. Springer.
  * MacDonald et al. (2017) MacDonald, M., Chung, D., Hutchins, N., Chan, L., Ooi, A. & García-Mayoral, R. 2017 The minimal-span channel for rough-wall turbulent flows. J. Fluid Mech. 816, 5–42.
  * MacDonald et al. (2019) MacDonald, M., Hutchins, N. & Chung, D. 2019 Roughness effects in turbulent forced convection. J. Fluid Mech. 861, 138–162.
  * MacDonald et al. (2018) MacDonald, M., Ooi, A., García-Mayoral, R., Hutchins, N. & Chung, D. 2018 Direct numerical simulation of high aspect ratio spanwise-aligned bars. J. Fluid Mech. 843, 126–155.
  * Marusic et al. (2021) Marusic, I., Chandran, D., Rouhi, A., Fu, M. K., Wine, D., Holloway, B., Chung, D. & Smits, A. J. 2021 An energy-efficient pathway to turbulent drag reduction. Nat. Commun. 12, 1–8.
  * Marusic et al. (2010) Marusic, I., Mathis, R. & Hutchins, N. 2010 Predictive model for wall-bounded turbulent flow. Science 329, 193–196.
  * Mathis et al. (2013) Mathis, R., Marusic, I., Chernyshenko, S. I & Hutchins, N. 2013 Estimating wall-shear-stress fluctuations given an outer region input. J. Fluid Mech. 715, 163–180.
  * Mathis et al. (2009) Mathis, R., Monty, J. P., Hutchins, N. & Marusic, I. 2009 Comparison of large-scale amplitude modulation in turbulent boundary layers, pipes, and channel flows. Phys. Fluids 21, 111703\.
  * Modesti et al. (2021) Modesti, D., Endrikat, S., Hutchins, N. & Chung, D. 2021 Dispersive stresses in turbulent flow over riblets. J. Fluid Mech. 917, A55.
  * Monty et al. (2009) Monty, J. P., Hutchins, N., Ng, H. C. H., Marusic, I. & Chong, M. S. 2009 A comparison of turbulent pipe, channel and boundary layer flows. J. Fluid Mech. 632, 431–442.
  * Musgrave & Tarazaga (2019) Musgrave, P. F. & Tarazaga, P. A. 2019 Turbulent boundary layer over a piezoelectrically excited traveling wave surface. In AIAA Scitech, p. 1354.
  * Nagib & Chauhan (2008) Nagib, H. M. & Chauhan, K. A. 2008 Variations of von Kármán coefficient in canonical flows. Phys. Fluids 20, 101518.
  * Park et al. (2004) Park, N., Yoo, J. Y. & Choi, H. 2004 Discretization errors in large eddy simulation: on the suitability of centered and upwind-biased compact difference schemes. J. Comput. Phys. 198, 580–616.
  * Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge University Press.
  * Ptasinski et al. (2003) Ptasinski, P. K., Boersma, B. J., Nieuwstadt, F. T. M., Hulsen, M. A., Van den Brule, B. H. A. A. & Hunt, J. C. R. 2003 Turbulent channel flow near maximum drag reduction: simulations, experiments and mechanisms. J. Fluid Mech. 490, 251–291.
  * Quadrio (2011) Quadrio, M. 2011 Drag reduction in turbulent boundary layers by in-plane wall motion. Phil. Trans. R. Soc. A 369, 1428–1442.
  * Quadrio & Ricco (2011) Quadrio, M. & Ricco, P. 2011 The laminar generalized Stokes layer and turbulent drag reduction. J. Fluid Mech. 667, 135–157.
  * Quadrio et al. (2009) Quadrio, M., Ricco, P. & Viotti, C. 2009 Streamwise-travelling waves of spanwise wall velocity for turbulent drag reduction. J. Fluid Mech. 627, 161–178.
  * Quadrio & Sibilla (2000) Quadrio, M. & Sibilla, S. 2000 Numerical simulation of turbulent flow in a pipe oscillating around its axis. J. Fluid Mech. 424, 217–241.
  * Ricco et al. (2012) Ricco, P., Ottonelli, C., Hasegawa, Y. & Quadrio, M. 2012 Changes in turbulent dissipation in a channel flow with oscillating walls. J. Fluid Mech. 700, 77.
  * Ricco et al. (2021) Ricco, P., Skote, M. & Leschziner, M. A. 2021 A review of turbulent skin-friction drag reduction by near-wall transverse forcing. Prog. Aero. Sci. 123, 100713.
  * Ricco & Wu (2004) Ricco, P. & Wu, S. 2004 On the effects of lateral wall oscillations on a turbulent boundary layer. Exp. Therm Fluid Sci. 29, 41–52.
  * Rouhi et al. (2022) Rouhi, A., Endrikat, S., Modesti, D., Sandberg, R. D., Oda, T., Tanimoto, K., Hutchins, N. & Chung, D. 2022 Riblet-generated flow mechanisms that lead to local breaking of reynolds analogy. J. Fluid Mech. 951, A45.
  * Schatzman et al. (2014) Schatzman, D., Wilson, J., Arad, E., Seifert, A. & Shtendel, T. 2014 Drag-reduction mechanisms of suction-and-oscillatory-blowing flow control. AIAA J. 52, 2491–2505.
  * Segawa et al. (2007) Segawa, T., Mizunuma, H., Murakami, K., Li, F. C. & Yoshida, H. 2007 Turbulent drag reduction by means of alternating suction and blowing jets. Fluid Dyn. Res. 39, 552\.
  * Skote (2022) Skote, M. 2022 Drag reduction of turbulent boundary layers by travelling and non-travelling waves of spanwise wall oscillations. Fluids 7, 65.
  * Smits et al. (2011) Smits, A. J., McKeon, B. J. & Marusic, I. 2011 High–Reynolds number wall turbulence. Annu. Rev. Fluid Mech. 43, 353–375.
  * Squire et al. (2016) Squire, D. T., Morrill-Winter, C., Hutchins, N., Schultz, M. P., Klewicki, J. C. & Marusic, I. 2016 Comparison of turbulent boundary layers over smooth and rough surfaces up to high Reynolds numbers. J. Fluid Mech. 795, 210–240.
  * Taylor (1938) Taylor, G. I. 1938 The spectrum of turbulence. Proc. R. Soc. London, Ser. A 164, 476–490.
  * Touber & Leschziner (2012) Touber, E. & Leschziner, M. A. 2012 Near-wall streak modification by spanwise oscillatory wall motion and drag-reduction mechanisms. J. Fluid Mech. 693, 150–200.
  * Viotti et al. (2009) Viotti, C., Quadrio, M. & Luchini, P. 2009 Streamwise oscillation of spanwise velocity at the wall of a channel for turbulent drag reduction. Phys. Fluids 21, 115109.
  * White & Mungal (2008) White, C. M. & Mungal, M. G. 2008 Mechanics and prediction of turbulent drag reduction with polymer additives. Annu. Rev. Fluid Mech. 40, 235–256.
  * Wu & Piomelli (2018) Wu, W. & Piomelli, U. 2018 Effects of surface roughness on a separating turbulent boundary layer. J. Fluid Mech. 841, 552–580.
  * Yamamoto et al. (2013) Yamamoto, A., Hasegawa, Y. & Kasagi, N. 2013 Optimal control of dissimilar heat and momentum transfer in a fully developed turbulent channel flow. J. Fluid Mech. 733, 189–220.
  * Yao et al. (2022) Yao, J., Chen, X. & Hussain, F. 2022 Direct numerical simulation of turbulent open channel flows at moderately high reynolds numbers. J. Fluid Mech. 953, A19.