1. Introduction
Flow separation occurs when the flow streamlines are not able to follow the contour of a solid surface and depart significantly from it. The problem is of interest as flows around large vehicles and buildings often experience separation (Oliveira & Younis Reference Oliveira and Younis2000; Choi, Lee & Park Reference Choi, Lee and Park2014).
In many engineering systems, separation causes a large increase in pressure drag that has to be overcome by an undesirable power expenditure, thus contributing negatively to pressing environmental issues, such as excessive fuel consumption, noise and pollutant emissions. The development of flow control techniques that aim to reduce the pressure drag due to flow separation is an extremely active area of academic and industrial research. Control methods can be classified into active and passive, depending on whether energy is supplied or not to the physical system. A clear advantage of active controllers is the large beneficial flow modifications that can be achieved. Active controllers can be deactivated when their action is no longer required, but a weight penalty is introduced due to the actuators. Active methods are also prone to damage and often require expensive maintenance. Passive methods often involve a geometrical modification of the surface, offer the benefit of not requiring external power supply, although the small wall alterations may necessitate regular repairment. The development of systems which, at the same time, guarantee a sizeable drag-reduction margin and low costs for production, maintenance and energy consumption is therefore of paramount importance to achieve economic advantages.
1.1. Active methods
1.1.1. Suction and blowing
The suction of fluid from a slit on the surface was the first active technique used to control flow separation (Chang Reference Chang2014). The amount of control is typically quantified via a mass flow rate coefficient 
$C_m$ (Nuber & Needham Reference Nuber and Needham1948) or a momentum coefficient 
$C_\mu$ (Chng et al. Reference Chng, Rachman, Tsai and Zha2009). Nuber & Needham (Reference Nuber and Needham1948) found a 
$25\,\%$ increase of the maximum lift coefficient for a NACA 
$64_1A212$ airfoil by applying suction with a coefficient 
$C_m=1.8\times 10^{-3}$ along the upper surface near the leading edge. Seifert & Pack (Reference Seifert and Pack2002) applied steady and oscillating suction and blowing near the separation point of a Glauert–Goldschmied airfoil. Oscillating suction and blowing effectively reduced the pressure drag for 
$C_\mu <0.1\,\%$. Separation drag was suppressed by steady suction or blowing using 
$C_\mu =0.8\,\%$ and 
$C_\mu =2\,\%$ respectively, with suction leading to the best performance. Greenblatt et al. (Reference Greenblatt, Paschal, Yao, Harris, Schaeffler and Washburn2006) removed the fluid around the upper surface of a modified Glauert profile and the form drag was reduced by increasing 
$C_\mu$, showing a strong dependence on the Reynolds number. In the same study, the separation bubble was suppressed when 
$C_\mu =2\,\%$.
Some authors applied blowing over the separation line of supersonic airfoils (Bradley & Wray Reference Bradley and Wray1974; Meyer & Seginer Reference Meyer and Seginer1994; Wong & Kontis Reference Wong and Kontis2007) or backward-facing steps (Chun, Lee & Sung Reference Chun, Lee and Sung1999), effectively reducing the reattachment length. More recent studies, such as Brackston, Wynn & Morrison (Reference Brackston, Wynn and Morrison2016b), have employed extremum-seeking feedback control, implemented via pulsed-jet actuators, to minimize the separation drag given by a bluff-body wake. Cho, Choi & Choi (Reference Cho, Choi and Choi2016) reduced the extent of the separation bubble in a turbulent boundary layer by means of oscillating suction and blowing. The control, applied near the separation line, gave the best performance for a frequency 
$f=0.5$, normalized via the free-stream velocity and the uncontrolled bubble length. The creation of intense spanwise vortices energized the fluid near the wall, thereby reducing the bubble.
Critical aspects of these active methods are the minimization of the pumping power required by the actuation, the skin-friction penalty (Gad-el Hak, Pollard & Bonnet Reference Gad-el Hak, Pollard and Bonnet1998) due to the enhanced fluid mixing, and the decrease of thrust due to the air drawn from the aircraft engines (Chng et al. Reference Chng, Rachman, Tsai and Zha2009).
1.1.2. Plasma actuators
Plasma actuators are devices composed of two metal electrodes separated by a dielectric material and subjected to a high voltage. The fluid is ionized around the actuators, by creating a volume force equivalent to a pressure gradient (Corke et al. Reference Corke, Jumper, Post, Orlov and McLaughlin2002).
Post & Corke (Reference Post and Corke2004) investigated the ability of plasma actuators to reattach the flow over a NACA
$66_3-018$ airfoil at the post-stall incidence. The actuators, applied at the leading edge of the airfoil and at the maximum chamber location, created a two-dimensional steady flow field that energized the boundary layer around the body. The lift increased with the forcing magnitude up to a peak value corresponding to the fluid reattachment, while the drag reduced proportionally to the actuation intensity, causing a maximum 
$400\,\%$ increase of the aerodynamic efficiency. Post & Corke (Reference Post and Corke2006) controlled the dynamic stall and separation on a NACA
$0015$ airfoil oscillating around its axis by employing three types of plasma actuators that all reattached the flow. The steady actuator increased the lift, except at high angles of attack, while the unsteady actuator raised the lift only during the pitch-down phases. Huang, Corke & Thomas (Reference Huang, Corke and Thomas2006) employed a plasma actuator to control separation on a blade cascade. The actuator was located at two positions upstream of the separation point and led to a reduction of the reattachment length, similar to that obtained with vortex generators in the same configuration. The main drawback of plasma actuators is the energy used for their activation, causing a low efficiency at high Reynolds numbers (Neretti Reference Neretti2016).
1.2. Passive methods
Vortex generators are the most common passive devices for separation control. They are in the form of geometrical modifications, usually placed upstream of the region of separation. They create streamwise vortices that enhance the turbulent mixing and allow the flow to sustain adverse pressure gradients, thereby reducing the tendency of the flow to separate. Calarese, Crisler & Gustafson (Reference Calarese, Crisler and Gustafson1985) conducted experiments on a replica of a C-130 aircraft with arrays of vortex generators placed along the fuselage circumference. Vortex generators at the most upstream position were most effective at an incidence of 
$4^{\circ }$, leading a 
$7\,\%$ reduction of the total aircraft drag. Bragg & Gregorek (Reference Bragg and Gregorek1987) applied three kind of vortex generators on a canard wing to alleviate the effect of boundary-layer separation caused by dirt deposition. All vortex generators were able to increase the lift.
Classical vortex generators create parasite drag even when no separation occurs, and therefore research has focused on developing devices that are small enough to be completely submerged in the boundary layer. Lin, Selby & Howard (Reference Lin, Selby and Howard1991) proved that wishbone vortex generators with a height of 
$20\,\%$ of the boundary-layer thickness were effective in reducing the reattachment length in a descending ramp. They created a pair of counter-rotating vortices that enhanced mixing in the area with the lowest momentum. Lin et al. (Reference Lin, Robinson, McGhee and Valarezo1994) experimentally studied how delta and trapezoidal submerged vortex generators influenced the flow on a three-section airfoil with a flap. At low incidence (
$8^{\circ }$) vortex generators were placed at a quarter of the flap chord and the separation around the flap was controlled effectively as the wake thickness was reduced by 
$60\,\%$. At the same angle of attack they increased the lift by 
$11\,\%$ and generated a drag reduction of 
$38\,\%$ with an improvement of 
$80\,\%$ on the aerodynamic efficiency.
The main shortcoming of submerged vortex generators is their tendency to break because of their small size.
1.3. Objectives and structure of the paper
The aim of this numerical work is to study the effectiveness of a spanwise oscillating pressure gradient in controlling flow separation. The objectives are to reduce the pressure drag caused by the flow separation around square bars located on the opposite walls of a channel and to understand the physics behind the flow alterations induced by the spanwise oscillating pressure gradient.
Spanwise oscillating walls and pressure gradients have been used extensively for turbulent skin-friction drag reduction (Jung, Mangiavacchi & Akhavan Reference Jung, Mangiavacchi and Akhavan1992; Trujillo, Bogard & Ball Reference Trujillo, Bogard and Ball1997; Quadrio & Ricco Reference Quadrio and Ricco2004), but only a very few studies have focused on the effect of spanwise forcing on large flow separation. Jukes & Choi (Reference Jukes and Choi2012) were able to reattach the flow on a 
$20^{\circ }$ inclined ramp by using spanwise jets, while Brackston et al. (Reference Brackston, García De La Cruz, Wynn, Rigas and Morrison2016a) and Brackston, Wynn & Morrison (Reference Brackston, Wynn and Morrison2018) used oscillating flaps positioned at the back sides of a bluff body in a wind tunnel to control the separation wake via reduced-order feedback control. The recent numerical study of Banchetti, Luchini & Quadrio (Reference Banchetti, Luchini and Quadrio2020) showed that streamwise travelling waves of spanwise wall velocity can reduce the pressure drag of flows separating marginally behind smooth obstacles.
In the present study, square bars are chosen because the pressure drag and the skin-friction drag are completely distinct around these obstacles. The pressure drag is produced by the difference in the surface integrated pressures in front and behind the bars and the skin-friction drag is the result of the surface integrated wall-shear stress on the bar crests and over the channel walls between bars. The contribution of the oscillating pressure gradient to the global drag reduction can therefore be precisely quantified by the separate reductions of the pressure drag and the skin-friction drag. For the chosen geometry, the pressure drag contributes significantly to the total drag as the ratio between the pressure drag and the skin-friction drag is 4.3, that is, the pressure drag is 81 % of the total drag. The selection of an uncontrolled flow with such a feature is central to our study because we are particularly interested in the effect of spanwise forcing on the separated flow behind the square bars.
Section 2 illustrates the flow configuration and discusses the numerical method and the parameters used to evaluate the control performance. Section 3 describes the derivation and the numerical solution of the laminar auxiliary problem employed to compute the power spent for the actuation. In § 4 the mean pressure and the mean streamlines of the uncontrolled flow are first discussed and the drag-reduction margin, the power spent for the control actuation and the net power are quantified. In § 5, the effect of the actuation on the flow statistics is further investigated and the role of the Reynolds stresses between the bars is outlined. Section 6 presents the conclusions of our work.
2. Numerical procedures
We study the turbulent channel flow depicted in figure 1 where, within the computational domain, two identical square bars are located on the opposite walls and oriented perpendicularly to the mean streamwise direction.
Figure 1. Sketch of the system configuration. The bar is indicated by the grey square and the drag coefficients, assigned to each surface, are defined in § 2.3. The pressure coefficient 
$C_p$ is computed by considering both sides of the bar. Only the bottom half of the channel is shown. The dash-dotted line represents the symmetry plane at 
$y^*=\delta ^*+h^*$ and the height of the bar is 
$h^*=0.2\delta ^*$. The distance between the bar and the domain sides is not drawn to scale.
The boundary conditions are periodic over the streamwise and spanwise directions, while the no-slip condition is imposed on the solid walls. Dimensional quantities are indicated with an asterisk. The Cartesian coordinates are 
$x^*$, 
$y^*$ and 
$z^*$ along the streamwise, vertical and spanwise directions respectively, where 
$x^*$ is defined from the upstream start of the computational domain and 
$y^*$ is defined from the bottom channel wall. The side of the square section of a bar is of length 
$h^*$ and 
$\delta ^*=5 h^*$ is half of the distance between the bar crests. The centre of each bar is located at 
$x^*=4\delta ^*$.
The flow used to initiate the computations is of zero velocity in the two strips between the bars, 
$0 < y^* < h^*$ and 
$2\delta ^*+h^* < y^* < 2(\delta ^*+h^*)$, and a laminar Poiseuille flow plus a random three-dimensional perturbation velocity field of zero mean in the range 
$h^* < y^* < 2\delta ^*+h^*$. The centreline velocity of this Poiseuille flow is 
$U_p^*$. The channel studied in this work differs from the one of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003), where multiple bars were used, because in our case only one bar is placed on each wall within the computational domain and symmetry between the top and bottom walls applies. The distance between bars on the same wall is 
$39 h^*$, therefore higher than those used by Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) (
$0.33h^*-19 h^*$). We have chosen this large distance because Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) showed that there is no influence between the recirculation regions for distances larger than 
$7 h^*$. All the variables are scaled using 
$\delta ^*$ and 
$U_p^*$ and the Reynolds number is 
$Re_p=U_p^* \delta ^*/\nu ^*=4200$, where 
$\nu ^*$ is the kinematic viscosity of the fluid. The simulations are performed at a constant mass flow rate. The components of the velocity vector along 
$x^*$, 
$y^*$ and 
$z^*$ are indicated by 
$u^*$, 
$v^*$ and 
$w^*$, respectively. The dimensions of the computational domain are 
$L_x=8$, 
$L_y=2.4$ and 
$L_z={\rm \pi}$, where the subscripts denote the Cartesian coordinates. The pressure 
$p^*$, scaled using the density of the fluid 
$\rho ^*$ as 
$p=p^*/(\rho ^*{U_p^*}^2)$, is defined as 
$p(x,y,z,t)=\phi (x,y,z,t)+\varPi _x(t)\, (x-L_x)+\varPi _z(t)\, z$, where 
$\phi$ is a periodic function of 
$x$ and 
$z$, 
$\varPi _x <0$ is the spanwise-averaged streamwise pressure gradient and 
$\varPi _z$ is the spanwise pressure gradient used for controlling the flow. The latter has the form
\begin{equation} \varPi^*_z=A^* \cos\left(\frac{2{\rm \pi} t^*}{T^*}\right), \end{equation}
where 
$A^*$ and 
$T^*$ are the amplitude and the period, scaled as 
$A=A^*\delta ^*/ (\rho ^* {U^*_p}^2)$ and 
$T=T^* U^*_p/ \delta ^*$. Analogous to the studies on the wall-oscillation technique where the peak-to-peak displacement is used as a third forcing parameter (Ricco & Quadrio Reference Ricco and Quadrio2008), we define the peak-to-peak impulse 
$J=J^*/\rho ^* U_p^*=2\max _t(\int _0^t \varPi _z \,\mathrm {d}{\hat t})=AT/{\rm \pi}$ as an additional parameter related to 
$A$ and 
$T$. The impulse 
$J^*$ has units 
$\textrm {kg}\,\textrm {m}^{-2}\,\textrm {s}$ and therefore it represents the change of momentum along 
$z$ per unit volume imposed by the pressure gradient.
2.1. Direct numerical simulations
The fluid flow is simulated by the spectral element code Nek5000 (Fischer Reference Fischer2017). The code employs the spectral finite element method (SEM), based on the Galerkin discretization of the weak form of the Navier–Stokes equations. The fluid domain is split into hexahedral elements where each fluid quantity 
$q$ is expressed in terms of Lagrange polynomials with nodes following the distribution of Gauss–Lobatto–Legendre (GLL) points. Each element is discretized with 
$l_x l_y l_z=l_x^3$ GLL points, where 
$l_x=l_y=l_z$ are the number of points used for each edge of the element. The Navier–Stokes equations are discretized in time using an operator integration factor splitting scheme (Maday, Patera & Rønquist Reference Maday, Patera and Rønquist1990) that permits large time steps 
$\Delta t$ due to its semi-implicit nature. The time step is set to the constant value 
$\Delta t= 1.69\times 10^{-4}$ and the related Courant–Friedrichs–Lewy number is in the range 1.5–2.0. Figure 2 depicts the mesh around the lower bar. The mesh radiates from the bar to guarantee a high clustering of points especially around the left and the right top corners, where singular behaviours of the local pressure and the local wall-shear stress occur. A total number of 
$76\,496$ hexahedral elements are employed. The mesh is identical on each 
$x$–
$y$ plane and 14 identical elements are used along 
$z$. A number of 
$6$ GLL points is adequate to compute the total drag accurately, as shown by the convergence study in appendix A. The first point above the cavity is at a distance 
$\Delta y=1.1\times 10^{-3}$ from the wall, while the last point along 
$y$ is at a distance 
$\Delta y=1.3\times 10^{-2}$ from the centreline. The maximum streamwise resolution is 
$\Delta x=1.5\times 10^{-2}$. The minimal 
$\Delta x=\Delta y =5.7\times 10^{-5}$ refer to the grid elements at the side edges and the crests of the bars, respectively. Even though the spectral elements along 
$z$ are identical, the spanwise resolution varies from 
$\Delta z_{min}=2.6\times 10^{-2}$ to 
$\Delta z_{max}=6.4\times 10^{-2}$. All the simulations are performed on the high performance computing facility of the A*STAR Computational Resource Centre in Singapore, employing 2000–8000 processors. The flow statistics and the drag coefficients are computed after the end of an initial transient time, evaluated by the inspection of the time history of total drag. This is crucial to assess that the flow has reached a statistically steady state.
Figure 2. Mesh for the channel with bars. The grid is symmetrical around the axes 
$x^*=4\delta ^*$ and 
$y^*=1.2\delta ^*$ and only the quarter of the computational domain in the range 
$\left [ 0,4\delta ^* \right ] \times \left [ 0,1.2\delta ^* \right ]$ is shown. Only the hexahedral elements are shown.
2.2. Averaging operators and flow decomposition
The average of a quantity 
$q$ along the spanwise direction is
\begin{equation} \left \langle q \right\rangle(x,y,t)=\frac{1}{L_z}\int_0^{L_z} q(x,y,\hat{z},t) \, \mathrm{d} \hat{z}, \end{equation}where the integral is approximated via the Gauss–Legendre quadrature that takes into account the distribution of the GLL points. The time-ensemble average is
\begin{equation} \hat{q}\ (x,y,z,\tau)=\frac{1}{N}\sum_{n=0}^{N-1} q \left( x,y,z,\tau+n \, T \right), \quad 0\leq \tau< T, \end{equation}
where 
$N$ is the number of periods, and the time average between times 
$t_i$ and 
$t_f$ is
\begin{equation} \bar{q}(x,y,z)=\frac{1}{t_f-t_i}\int_{t_i}^{t_f} q\left(x,y,z,\hat{t}\ \right) \,\mathrm{d} \hat{t}. \end{equation}
A quantity 
$q$ of the flow field is decomposed as
\begin{equation} q=q_m\left( x,y \right) +q_{osc}\left( x,y,\tau \right) +q_t\left( x,y,z,t \right), \end{equation}
where 
$q_m\left ( x,y \right ) =\left \langle \bar {q} \right \rangle$ is the mean flow, 
$q_{osc}\left ( x,y,\tau \right ) =\left \langle \hat {q} \right \rangle -q_m$ is the time-periodic flow and 
$q_t=q-\left \langle \hat {q} \right \rangle$ is the turbulent fluctuation.
2.3. Parameters defining the actuation performance
For a separated flow the total drag force is composed of a pressure-drag force 
$F_p^*$ and a skin-friction drag force 
$F_f^*$, defined as
\begin{equation} F_p^*=\int_{S^*} p^*\, \boldsymbol{i} \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S^*, \quad F_f^*=-\int_{S^*} \left.\mu^* \frac{\partial u^*}{\partial \boldsymbol{n}}\right\vert_w \, \mathrm{d} S^*, \end{equation}
where 
$\boldsymbol {n}$ is the normal unit vector pointing out of the fluid domain, 
$\boldsymbol {i}$ is the unit vector along the 
$x$ direction, 
$S^*$ is the solid surface and the 
$w$ denotes the wall. For the vertical surfaces of the bar 
$\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {i}\neq 0$, and therefore the pressure drag force is
\begin{equation} F_p^*=\int_{S^*_l} \overline{p^*} \, \mathrm{d} S^* -\int_{S^*_r} \overline{p^*} \, \mathrm{d} S^*, \end{equation}
where 
$S^*_l$ and 
$S^*_r$ are the left and right vertical surfaces of the bar of area 
$h^*L_z^*$. The skin-friction drag is
\begin{equation} F_f^*=F^*_{ca}+F^*_{cr}=\int_{S^*_{ca}} \left.\mu^* \frac{\partial \overline{u^*}}{\partial y^*}\right|_w \, \mathrm{d} S^* +\int_{S^*_{cr}} \left.\mu^* \frac{\partial \overline{u^*}}{\partial y^*}\right|_w \, \mathrm{d} S^*, \end{equation}
where 
$S^*_{ca}$ and 
$S^*_{cr}$ are the cavity and crest surfaces, respectively. Note that the minus sign in (2.6a,b) is absent in (2.8) because, for the bottom horizontal wall, 
$\boldsymbol {n}=-\boldsymbol {j}$, where 
$\boldsymbol {j}$ is the unit vector along the 
$y$ axis. The forces 
$F^*_{ca}$, 
$F^*_{cr}$ and 
$F_p^*$ include contributions from the bottom and top walls. The drag coefficients of the skin friction along the cavity (
$C_{ca}$), the skin friction along the crest (
$C_{cr}$), the pressure (
$C_p$) are defined by dividing 
$F^*_{ca}$, 
$F^*_{cr}$ and 
$F_p^*$ by 
$\rho ^*{U^*_p}^2 L^*_xL^*_z$, respectively. The skin-friction coefficient is composed of the crest and the cavity drag coefficients, i.e. 
$C_f=C_{cr}+C_{ca}$. We do not define the drag coefficients on the left and the right bar sides and their reductions because, as remarked by Banchetti et al. (Reference Banchetti, Luchini and Quadrio2020), it is not meaningful to define a local pressure coefficient for an incompressible flow as the pressure values can be shifted by any constant.
We define a percentage reduction for each drag coefficient with respect to its value in the uncontrolled case, denoting these reductions with the letter 
$\mathcal {R}$ and the same subscript of the drag coefficients. For instance, the percentage reduction of the skin-friction coefficient along the cavity is
\begin{equation} \mathcal{R}_{ca}(\%)=100 \left( 1-\frac{{C_{ca}}_o}{C_{ca}} \right), \end{equation}
where 
$C_{ca}$ and 
${C_{ca}}_o$ are the uncontrolled and controlled values, respectively (the subscript 
$o$ stands for oscillating). The total drag reduction 
$\mathcal {R}(\%)$ is defined as
\begin{equation} \mathcal{R}(\%) = 100 \left( 1- \frac{{F_p}_{o}^*+{F_f}_{o}^*}{F_p^*+F_f^*}\right)= 100 \left( 1- \frac{{C_f}_{o}+{C_p}_{o}}{{C_f}+{C_p}} \right), \end{equation}
where the pressure-drag force 
$F_p^*$ is defined in (2.7) and the skin-friction drag force 
$F_f^*$ is defined in (2.8). Banchetti et al. (Reference Banchetti, Luchini and Quadrio2020) fixed the total pressure at the inflow location of their computational domain in order to quantify the effect of the pressure drop given by their control method. In all our cases, we instead maintain the pressure constant at the further end of our computational domain.
As we study an active technique, it is crucial to account for the power 
$\mathcal {P}_z(\%)$ supplied to the fluid for controlling the flow, defined as
\begin{equation} \mathcal{P}_z(\%)=\frac{100}{\bar{\varPi}_x U_b (t_f-t_i)} \int_{t_i}^{t_f} \varPi_z\left( \hat{t} \right)\, W_b\left( \hat{t} \right) \, {\mathrm{d} \hat{t}}, \end{equation}
expressed as a percentage of the power used for driving the fluid along the streamwise direction. In (2.11) 
$\bar {\varPi }_x$ is the mean pressure gradient driving the flow along 
$x$ in the uncontrolled case, 
$U_b$ is the constant volume-averaged velocity, 
$U_b=V^{-1} \int _{V} \bar {u} \, \mathrm {d} V$, 
$W_b(t)$ is the spanwise volume-averaged velocity due to the control, 
$W_b=V^{-1} \int _{V} w(x,y,z,t) \, \mathrm {d} V$ and 
$V$ is the fluid volume. Appendix B presents the derivation for the power used to drive the fluid along 
$x$, 
${P}_x=-\bar {\varPi }_x U_b>0$.
It is verified numerically that the force exerted by the pressure gradient to drive the fluid along the streamwise direction matches the sum of the drag forces, i.e. 
$-\bar {\varPi }_x L_y=C_f+C_p$ and 
$-{\bar {\varPi }_x}_o L_y={C_f}_o+{C_p}_o$. By using this force balance, the percentage of power saved is
\begin{equation} \mathcal{P}_{sav}(\%) = 100 \left( 1 - \frac{\varPi_{x_o} U_b}{\bar{\varPi}_{x} U_b}\right) = 100 \left( 1- \frac{\bar{\varPi}_{x_o}}{\bar{\varPi}_x} \right) = 100 \left( 1 - \frac{C_{f_o} + C_{p_o}}{C_f + C_p} \right) =\mathcal{R}(\%), \tag{2.12}\end{equation}
because the volume-averaged velocity 
$U_b$ is constant.
The difference between the power 
$\mathcal {P}_{sav}(\%)$ saved through the oscillating pressure gradient and the power 
$\mathcal {P}_z(\%)$ spent for controlling the flow is the net power saved
\begin{equation} \mathcal{P}_{net}(\%)=\mathcal{R}(\%)-\mathcal{P}_z(\%). \end{equation}3. Laminar-flow results
Previous studies (Ricco & Quadrio Reference Ricco and Quadrio2008; Quadrio & Ricco Reference Quadrio and Ricco2011; Ricco & Hahn Reference Ricco and Hahn2013) showed that useful information about control in the turbulent case can be extracted from the corresponding laminar-flow solutions. Here, for the first time, we apply the same strategy to the study of a separated flow. We consider the time-ensemble and 
$z$-averaged spanwise momentum equation
\begin{align} \frac{\partial \left \langle \widehat{w} \right\rangle}{\partial \tau}&+\frac{\partial \left \langle \widehat{u w} \right\rangle}{\partial x}+ \frac{\partial \left \langle \widehat{v w} \right\rangle}{\partial y}=-\frac{1}{L_z}\left[ \hat{p}(x,y,L_z,\tau)-\hat{p}(x,y,0,\tau)\right]\nonumber\\ &+\frac{1}{Re_p}\left( \frac{\partial^2 \left \langle \widehat{w} \right\rangle}{\partial x^2}+\frac{\partial^2 \left \langle \widehat{w} \right\rangle}{\partial y^2} \right), \end{align}where
\begin{equation} -\frac{1}{L_z}\left[ \hat{p}(x,y,L_z,\tau)-\hat{p}(x,y,0,\tau)\right] = -\left\langle \frac{\partial \hat{p}}{\partial z} \right\rangle. \end{equation}
Substitution of the pressure 
$p=\phi +\bar {\varPi }_x\, (L_x - x) + \varPi _z\, z$, the Reynolds decomposition 
$w=\left \langle \widehat {w} \right \rangle +w_t$ and the 
$z$-averaged continuity equation
\begin{equation} \frac{\partial \left \langle \widehat{u} \right\rangle}{\partial x}+\frac{\partial \left \langle \widehat{v} \right\rangle}{\partial y}=0 \end{equation}into (3.1) leads to
\begin{equation} \frac{\partial \left \langle \widehat{w} \right\rangle}{\partial \tau}=\frac{1}{Re_p}\left( \frac{\partial^2 \left \langle \widehat{w} \right\rangle}{\partial x^2}+\frac{\partial^2 \left \langle \widehat{w} \right\rangle}{\partial y^2} \right) -\varPi_z-\left \langle \widehat{u} \right\rangle\frac{\partial \left \langle \widehat{w} \right\rangle}{\partial x} -\left \langle \widehat{v} \right\rangle\frac{\partial \left \langle \widehat{w} \right\rangle}{\partial y} -\frac{\partial \left \langle \widehat{u_t w_t} \right\rangle}{\partial x}- \frac{\partial \left \langle \widehat{v_t w_t} \right\rangle}{\partial y}. \end{equation} Ricco & Quadrio (Reference Ricco and Quadrio2008) showed that, for a smooth channel with oscillating walls, there is a very good agreement between the 
$x$- and 
$z$-averaged spanwise velocity profile and the laminar Stokes layer profile because the term 
$\partial \left \langle \widehat {v_t w_t} \right \rangle /\partial y$ is negligible in the spanwise momentum equation. Following Ricco & Quadrio (Reference Ricco and Quadrio2008), we assume that 
$\partial \left \langle \widehat {v_t w_t} \right \rangle / \partial y$ and 
$\partial \left \langle \widehat {u_t w_t} \right \rangle /\partial x$ are negligible and obtain
\begin{equation} \frac{\partial w_l}{\partial \tau}=\frac{1}{Re_p} \left( \frac{\partial^2 w_l}{\partial x^2}+\frac{\partial^2 w_l}{\partial y^2} \right) -\varPi_z- \left( \left \langle \widehat{u} \right\rangle \frac{\partial w_l}{\partial x} +\left \langle \widehat{v} \right\rangle \frac{\partial w_l}{\partial y} \right), \end{equation}
where 
$w_l$ is a laminar approximation to 
$\left \langle \widehat {w} \right \rangle$. We further note that the velocities 
$\left \langle \widehat {u} \right \rangle$ and 
$\left \langle \widehat {v} \right \rangle$ are independent of 
$w_l$, as shown by the time-ensemble and 
$z$-averaged momentum equations,
\begin{gather} \frac{\partial \left \langle \widehat{u} \right\rangle}{\partial \tau}=\frac{1}{Re_p} \left( \frac{\partial^2 \left \langle \widehat{u} \right\rangle}{\partial x^2}+\frac{\partial^2 \left \langle \widehat{u} \right\rangle}{\partial y^2} \right) -\left( \frac{\partial \left \langle \widehat{uu} \right\rangle}{\partial x}+\frac{\partial \left \langle \widehat{uv} \right\rangle}{\partial y} \right) -\frac{\partial \left \langle \widehat{p} \right\rangle}{\partial x}, \end{gather}
\begin{gather} \frac{\partial \left \langle \widehat{v} \right\rangle}{\partial \tau}=\frac{1}{Re_p} \left( \frac{\partial^2 \left \langle \widehat{v} \right\rangle}{\partial x^2}+\frac{\partial^2 \left \langle \widehat{v} \right\rangle}{\partial y^2} \right) -\left( \frac{\partial \left \langle \widehat{uv} \right\rangle}{\partial x}+\frac{\partial \left \langle \widehat{vv} \right\rangle}{\partial y} \right) -\frac{\partial \left \langle \widehat{p} \right\rangle}{\partial y}. \end{gather}
We assume that 
$w_l$ is a periodic function of 
$\tau$ and thus Fourier series expansions can be employed. The quantities 
$\varPi _z$ and 
$w_l$ are expressed as
\begin{equation} w_l=\sum_{n=-\infty}^{+\infty}\frac{\hat{w}_n}{2} \textrm{e}^{\textrm{i} n \omega \tau}, \quad \varPi_z=A \cos\left( \omega \tau \right) =\frac{A}{2} \left( \textrm{e}^{\textrm{i}\omega \tau}+\textrm{e}^{-\textrm{i}\omega \tau} \right), \end{equation}
where 
$\omega =2{\rm \pi} /T$. Assuming that the mean velocities 
$\left \langle \widehat {u} \right \rangle$ and 
$\left \langle \widehat {v} \right \rangle$ do not depend on 
$\tau$, substitution of (3.7a,b) into (3.5) leads to
\begin{align} &\sum_{n=-\infty}^{+\infty} \left[ \frac{\textrm{i}\,n\omega}{2} \hat{w}_n -\frac{1}{Re_p}\left( \frac{\partial^2 \hat{w}_n}{\partial x^2}+\frac{\partial^2 \hat{w}_n}{\partial y^2} \right) +\left \langle \widehat{u} \right\rangle \frac{\partial \hat{w}_n}{\partial x}+\left \langle \widehat{v} \right\rangle\frac{\partial \hat{w}_n}{\partial y} \right] \textrm{e}^{\textrm{i}\,n\,\omega\,\tau}\nonumber\\ &\qquad\qquad \qquad\qquad= -\frac{A}{2} \left( \textrm{e}^{\textrm{i}\,\omega\,\tau}+\textrm{e}^{-\textrm{i}\,\omega\,\tau} \right). \end{align}
The assumptions of the neglect of the Reynolds stresses, the steadiness of 
$\left \langle \widehat {u} \right \rangle$ and 
$\left \langle \widehat {v} \right \rangle$ and the periodicity of 
$w_l$ on 
$\tau$ are confirmed by the turbulent-flow numerical calculations, discussed in § 5.3. The modes 
$\hat {w}_n$ satisfy
\begin{equation} \textrm{i}\,n\omega \hat{w}_n =\frac{1}{Re_p}\left( \frac{\partial^2 \hat{w}_n}{\partial x^2} + \frac{\partial^2 \hat{w}_n}{\partial y^2}\right) -\left \langle \widehat{u} \right\rangle\frac{\partial \hat{w}_n}{\partial x}-\left \langle \widehat{v} \right\rangle\frac{\partial \hat{w}_n}{\partial y}, \quad n\neq \pm 1, \end{equation}and
\begin{equation} \left( \pm \textrm{i}\,\omega \right) \hat{w}_{\pm 1} =\frac{1}{Re_p}\left( \frac{\partial^2 \hat{w}_{\pm 1}}{\partial x^2}+ \frac{\partial^2 \hat{w}_{\pm 1}}{\partial y^2} \right) -\left \langle \widehat{u} \right\rangle\frac{\partial \hat{w}_{\pm 1}}{\partial x}-\left \langle \widehat{v} \right\rangle\frac{\partial \hat{w}_{\pm 1}}{\partial y} -A. \end{equation}The boundary conditions for (3.9b) are
\begin{gather} \hat{w}_{\pm 1}\left( 0,y \right) =\hat{w}_{\pm 1}\left( h,y \right) =0, \quad 0\leq y\leq h, \end{gather}
\begin{gather} \hat{w}_{\pm 1}\left( x,h \right) =\hat{w}_{\pm 1}\left( x,0\right) =0, \quad h\leq x\leq L_x, \end{gather}
\begin{gather} \hat{w}_{\pm 1}\left( h,y\right) =0, \quad 0\leq y\leq h, \end{gather}
\begin{gather} \frac{\partial \hat{w}_{\pm 1}}{\partial y}\left( x,1+h \right) =0, \quad 0\leq x\leq L_x, \end{gather}
\begin{gather} \hat{w}_{\pm 1} \left( 0,y \right) =\hat{w}_{\pm 1} \left( L_x,y \right), \quad 0\leq y\leq 1+h. \end{gather} Conditions (3.10a), (3.10b), (3.10c) are the no-slip conditions at the walls, condition (3.10d) expresses the symmetry along the centreline and (3.10e) denotes the periodicity along 
$x$. The only non-null modes are those for which 
$n={\pm }1$ and from (3.7a,b) it follows that
\begin{equation} w_l\left( x,y,\tau \right) =\tfrac{1}{2}\left( \hat{w}_{1}\,\textrm{e}^{\textrm{i}\,\omega\tau}+\hat{w}_{-1}\,\textrm{e}^{-\textrm{i}\,\omega\tau} \right). \end{equation}
The spanwise velocity 
$w_l\left ( x,y,\tau \right )$ is real and thus 
$\hat {w}_{1}=\hat {w}_{-1}^{cc}$ because of the Hermitian property (the superscript 
$cc$ indicates the complex conjugate).
In order to eliminate the non-homogeneous term 
$A$ from (3.9b), we introduce 
$\tilde {w}_1\left ( x,y \right ) =\hat {w}_1\left ( x,y \right ) - A \textrm {i}/\omega$, which transforms (3.9b) into
\begin{equation} \frac{\partial^2 \tilde{w}}{\partial x^2}+ \frac{\partial^2 \tilde{w}}{\partial y^2}-\textrm{i}\,\omega\,Re_p \tilde{w}-Re_p\left( \left \langle \widehat{u} \right\rangle\frac{\partial \tilde{w}}{\partial x}+\left \langle \widehat{v} \right\rangle\frac{\partial \tilde{w}}{\partial y} \right) =0, \end{equation}
where the subscript of 
$\tilde {w}_1$ is removed for clarity. The boundary conditions for (3.12) are
\begin{gather} \tilde{w}\left( 0,y\right) =\tilde{w}\left( h,y\right) =-\frac{A}{\omega}\, \textrm{i}, \quad 0\leq y\leq h, \end{gather}
\begin{gather} \tilde{w}\left( x,h \right) =\tilde{w}\left( x,0 \right) =-\frac{A}{\omega}\, \textrm{i}, \quad h\leq x\leq L_x, \end{gather}
\begin{gather} \tilde{w}\left( h,y \right) =-\frac{A}{\omega}\, \textrm{i}, \quad 0\leq y\leq h, \end{gather}
\begin{gather} \frac{\partial \tilde{w}}{\partial y}\left( x,1+h\right) =0, \quad 0\leq x\leq L_x, \end{gather}
\begin{gather} \tilde{w} \left( 0,y\right) =\tilde{w} \left( L_x,y\right), \quad 0\leq y\leq 1+h. \end{gather}The laminar power spent is
\begin{equation} \mathcal{P}_{z,l}=-\frac{1}{T}\int_0^T \varPi_z(\tau)\; {W_b}_l(\tau) \, {\mathrm{d} \tau}, \end{equation}
where the laminar volume-averaged velocity 
${W_b}_l=V^{-1} \int _{V} w_l(x,y,\tau ) \, \mathrm {d} V$ can be simplified by using (3.11),
\begin{equation} {W_b}_l=\frac{1}{2 V}\int_{V} \left( \hat{w}_{1}\,\textrm{e}^{\textrm{i}\,\omega\tau}+\hat{w}^{cc}_{1}\,\textrm{e}^{-\textrm{i}\,\omega\tau} \right) \, \mathrm{d} V = \left( \hat{w}_b \textrm{e}^{\textrm{i}\,\omega\tau}+\hat{w}^{cc}_b \textrm{e}^{-\textrm{i}\,\omega\tau} \right) /2, \end{equation}
where 
$\hat {w}_b=V^{-1} \int _{V} \hat {w}_1 \, \mathrm {d} V$. Substitution of (3.15) and the second of (3.7a,b) into (3.14) leads to
\begin{equation} \mathcal{P}_{z,l}=-\frac{A}{4}\left( \hat{w}_b+\hat{w}^{cc}_b \right). \end{equation}
By introducing 
$\hat {w}_b=A\hat {w}_{u,b}$, where 
$\hat {w}_{u,b}$ is the value of 
$\hat {w}_b$ obtained when 
$A=1$, expression (3.16) can be written in terms of 
$\hat {w}_u$ as
\begin{equation} \mathcal{P}_{z,l}=-\frac{A^2}{4} \left( \hat{w}_{u,b} +\hat{w}^{cc}_{u,b} \right) =A^2 P_{z,lu}, \end{equation}
where 
$P_{z,lu}$ is the laminar power obtained for 
$A=1$. The laminar power (3.17) can also be expressed as the percentage of the power employed for driving the uncontrolled turbulent flow along 
$x$,
\begin{equation} \mathcal{P}_{z,l}(\%)=- \frac{100 \mathcal{P}_{z,l}}{\bar{\varPi}_x U_b}=\frac{25 A^2}{\bar{\varPi}_x U_b} \left( \hat{w}_{u,b} +\hat{w}^{cc}_{u,b} \right). \end{equation}
Equation (3.12) is solved numerically via a finite-difference method, as described in appendix C.. The laminar and the averaged turbulent-flow solutions are compared in § 5.3. We also solve equation (3.12) by neglecting the last two terms in parenthesis on the left-hand side. We compare the solutions in these two cases to quantify the effect of the advective terms on 
$w_l$ to find the locations where these terms are important and to determine whether they have an impact on the laminar power spent 
$\mathcal {P}_{z,l}$.
Figure 3 shows the absolute values of the differences between the real parts (3a) and the imaginary parts (3b) of the solutions obtained including and excluding the advective terms. The contours are more intense in front of the cavity and on the crest, while a mild disagreement occurs around the separated region.
Figure 3. Absolute values of the differences between the solution 
$\hat {w}_{\pm 1}$ of (3.9b) including or neglecting the advective terms for 
$A=1$ and 
$T=2.3$. The contour of the real parts is depicted in (a), while the contour of the imaginary parts is shown in (b).
For a smooth channel, it is possible to find an analytical expression relating the power 
$\mathcal {P}_{z,l}$ to the forcing parameters 
$A$ and 
$T$ (Ricco & Quadrio Reference Ricco and Quadrio2008). We herein use the analytical expression for 
$\mathcal {P}_{z,l}$ in the smooth channel case without bars to obtain an empirical relationship linking the power spent 
$\mathcal {P}_{z,l}$ to the period 
$T$. In the smooth channel case, the streamwise viscous term and the spatial advection terms are null. The streamwise viscous term is zero because the boundary conditions are independent of the 
$x$ direction. The laminar flow thus satisfies
\begin{equation} \left. \begin{gathered} \frac{\mathrm{d}^2 \tilde{w}}{{\textrm{d}y}^2}-\textrm{i}\,\omega Re_p\tilde{w}=0, \\ \tilde{w}(0)=\tilde{w}(2)=-\frac{A\,\textrm{i}}{\omega}. \end{gathered} \right\} \end{equation}
By imposing the boundary conditions on the general solution 
$\tilde {w}=C_1\textrm {e}^{\xi y}+C_2 \textrm {e}^{-\xi y}$ and using 
$\hat {w}_1=\tilde {w}+A \textrm {i}/\omega$, one finds
\begin{equation} \hat{w}=\frac{A \textrm{i}}{\omega} \left[ 1-\frac{\textrm{e}^{\xi (y-2)}+\textrm{e}^{-\xi y}}{1+\textrm{e}^{-2\xi}} \right], \end{equation}
where 
$\xi =\sqrt {\textrm {i}\omega Re_p}$. The integrated velocity 
$\hat {w}_b$ for the smooth channel is therefore
\begin{equation} \hat{w}_b=\frac{1}{2} \int_0^{2} \hat{w} \,{\textrm{d}y}=\int_0^{1} \hat{w} \,{\textrm{d}y}=\frac{A \textrm{i}}{\omega} \left[ 1- \frac{1-\textrm{e}^{-2 \xi}}{\xi(1+\textrm{e}^{-2 \xi})} \right]. \end{equation}
By noting that 
$\hat {w}_b+\hat {w}^{cc}_b$ is twice the real part of 
$\hat {w}_b$, we obtain
\begin{equation} \hat{w}_b+\hat{w}^{cc}_b=-\frac{2A}{\omega} {\rm Re} \left[ \frac{\textrm{i}\left(1-\textrm{e}^{-2 \xi}\right)}{\xi \left(1+\textrm{e}^{-2 \xi}\right)}\right] =-\frac{\sqrt{2}A}{\omega^{3/2} \sqrt{Re_p}}{\rm Re}\left[ (\textrm{i}+1) \frac{1-\textrm{e}^{-2 \xi}}{1+\textrm{e}^{-2 \xi}}\right], \end{equation}
where 
${\rm Re}$ indicates the real part. When 
$2\sqrt {\omega Re_p}\gg 1$, which is always satisfied in the studied cases, (3.22) can be approximated by
\begin{equation} \hat{w}_b+\hat{w}^{cc}_b\approx-\frac{\sqrt{2}A}{\omega^{3/2} \sqrt{Re_p}}. \end{equation}The laminar power spent for a smooth channel can be obtained by substituting (3.23) into (3.16) to find
\begin{equation} \mathcal{P}_{sc}=\frac{\sqrt{2} A^2\,\omega^{-3/2}}{4 \sqrt{Re_p}}=\frac{A^2 \, T^{3/2}}{8 {\rm \pi}^{3/2}\sqrt{Re_p}}, \end{equation}
where 
$sc$ stands for smooth channel. Equation (3.24) suggests that an expression proportional to 
$A^2 T^\gamma$ could be an approximation to 
$\mathcal {P}_{l}$. Best-line fitting of the laminar data points gives
\begin{equation} \mathcal{P}_{z,l}(\%) \approx -\frac{3.03\times 10^{-2} A^2 T^{3/2}}{\bar{\varPi}_x U_b}. \end{equation}
The same exponent 
$\gamma =3/2$ as in the smooth channel case is found. Expression (3.25) is plotted in figure 4 together with the values obtained via the laminar solution.The agreement is excellent, proving that the mean advection and the mean viscous diffusion along 
$x$ have a small impact on 
$\mathcal {P}_{z,l}$, as assumed when deriving (3.24). Ricco & Quadrio (Reference Ricco and Quadrio2008) found that the laminar power spent to oscillating the wall is proportional to 
$T^{-1/2}$, while we find that 
$\mathcal {P}_{z,l}$ is proportional to 
$T^{3/2}$ in our case. The reason of this apparent discrepancy is explained by considering a transformation relating the two laminar problems, as discussed in appendix D.
Figure 4. Power spent 
$\mathcal {P}_{l}/A^2$ as a function of 
$T$ for (3.12) with (squares) and without (circles) the mean advective terms. The solid line represents the function 
$\mathcal {P}_{l}/A^2 =3.03 \times 10^{-4} T^{1.5}$.
4. Turbulent-flow performance quantities
In this section, we present results on the uncontrolled turbulent flow and the performance parameters in the controlled cases.
4.1. Uncontrolled flow
The mean-flow streamlines in the region surrounding the bar in the uncontrolled flow case are depicted in figure 5. The separation zone is occupied by three recirculating vortices: two of them, A and B, are adjacent to the vertical surfaces of the bar and have a size comparable with the bar height. Vortex C is located downstream and on top of vortex B, extending for a length of more than five bar heights along the streamwise direction and covering the crest where a mild separation occurs. Vortices A and C rotate clockwise, while vortex B rotates anticlockwise. Sufficiently far from the bar (
$x<3$ and 
$x>6.5$) the flow is not influenced by the obstacle and the mean-flow streamlines become straight and aligned with the streamwise direction. We define the reattachment points at the positions where the time and spanwise-averaged wall-shear stress is zero. They separate the forward and backward mean flows near the wall and are found at 
$x=3.56, 4.39, 5.38$. The mean streamlines agree visually with those of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) and Ikeda & Durbin (Reference Ikeda and Durbin2007), although their Reynolds number is slightly larger.
Figure 5. Spanwise- and time-averaged streamlines for the flow without control. The recirculation vortices are indicated by the letters A, B and C.
Figure 6 shows the mean pressure of the uncontrolled flow. The pressure reaches high values right upstream of the bar where vortex A occurs, with the maximum pressure occurring along the left vertical surface of the bar. The core of vortex C, located approximately at 
$x=4.6$, is characterized by a low pressure. The effect of the bar on the wall-shear stress is clearly shown in figure 7 by the distribution of the local friction Reynolds number 
$Re_\tau =u_\tau ^* (h^*+\delta ^*)/\nu ^*$ as a function of the streamwise direction, where 
$u_\tau ^*=\sqrt {\tau _w^*/\rho ^*}$ is the local wall-friction velocity and 
$\tau _w^*$ is the time- and spanwise-averaged wall-shear stress. The Reynolds number 
$Re_\tau$ is only computed where the wall-shear stress is positive. It varies significantly in the proximity of the bar, while it is almost constant at 
$Re_{\tau }=218$ away from the bar, for 
$x<2$ and 
$x>6$. It must be remarked that scaling quantities with the local wall-friction velocity is only meaningful in the region where the flow is fully attached and relatively undisturbed from the separation regions so that the skin friction approaches a positive constant. It is instead not appropriate to employ wall-friction scaling in the recirculating regions where the skin friction is either vanishingly small or negative. The constant wall-friction velocity in the fully attached region, 
$u_\tau =0.044$, is used to scale quantities in wall units, as follows
\begin{gather} A^+ = \frac{A^* \nu^*}{\rho^* u_{\tau}^{*3}} = \frac{A(h+1)^3 Re_p^2}{h Re_\tau^3}, \end{gather}
\begin{gather}T^+ = \frac{T^* u_{\tau}^{*2}}{\nu^*} = \frac{T Re_\tau^2}{(h+1)^2 Re_p}, \end{gather}
\begin{gather}J^+ = \frac{J^*}{\rho^* u^*_{\tau}} = \frac{J (h+1) Re_p}{Re_\tau}. \end{gather}
Scaling the forcing parameters in wall units is useful to compare the performance of the actuation in reducing the skin-friction drag of the attached flow along the cavity with the widely studied spanwise wall-oscillation technique (Quadrio & Ricco Reference Quadrio and Ricco2004). We also find that the uncontrolled scaled wall-shear stress of the attached flow away from the bars, 
$Re_p^{-1}\mathrm {d}U/\mathrm {d}y\vert _{y=0}=0.0034$, is less than 1 % different from the coefficient found from the empirical correlation 
$Re_p^{-1}\mathrm {d}U/\mathrm {d}y\vert _{y=0}=0.0336 U_b^2 Re_\tau ^{-0.273}$ (Pope Reference Pope2000), valid for fully developed turbulent channel flows over flat walls.
Figure 6. Colour map and contour lines of the mean pressure in the uncontrolled case. The isoline values are equispaced at intervals of 
$0.02$ in the range 
$[-0.25,0.1]$.
Figure 7. Local friction Reynolds number along the channel in the uncontrolled case. The location of the bar is marked by the shaded region. The thin solid line denotes the friction Reynolds number 
$Re_{\tau }=218$, averaged in the constant-friction region, 
$x<2$ and 
$x>6.5$.
4.2. Skin-friction and pressure-drag coefficients
As the current control strategy aims at reducing both the viscous drag along the streamwise-parallel surfaces and the significant form drag experienced by the bar because of the recirculation regions, it is central to precisely quantify the drag coefficients and their changes in the controlled cases.
The effect of the control on the drag coefficients for a large drag-reduction case with 
$T=6.33$ and 
$A=1.08$ (
$J^+=50$, 
$T^+=50$) is presented in figure 8. The flow responds quickly to the control activation at 
$t=250$ and the coefficients adjust to the new value monotonically. This transient response is shorter and smoother than in the case of spanwise wall oscillations, where the wall-shear stress may evolve non-monotonically and reach a very low instantaneous value before adjusting to the new fully developed value (Quadrio & Ricco Reference Quadrio and Ricco2003). The time-averaged coefficients are computed after the end of the transient phase, estimated at around 
$t=270$, i.e. when all the coefficients reach a new statistically steady state. The control guarantees a substantial reduction of total drag 
$C_p+C_f$ due to a reduction of both 
$C_p$ and 
$C_f$. This result proves that a spanwise oscillating pressure gradient can reduce the skin-friction drag of a cavity flow and the pressure drag given by the difference on the bar sides in the integrated pressures of a fully separated flow. The drop of total pressure drag may not be achieved by the traditional separation control methods, such as vortex generators and jets, which may introduce a viscous drag penalty due to the enhanced flow mixing. The time histories of the pressure coefficient is much more oscillatory than that of the skin-friction coefficient, which is expected because the pressure drag is affected by the intense fluctuations in the separated region downstream of the bar. The skin-friction reduction is mostly due to the decrease of the cavity drag. The friction drag on the crest is very small and negative because of the reversed flow in the small recirculation region occurring due to separation. The spanwise forcing has a marginal effect on the crest friction.
Figure 8. (a) Time evolutions of the total drag coefficient 
$C_p+C_f$. (b) Time evolutions of the pressure-drag coefficients 
$C_{p}$ and the skin-friction coefficients 
$C_{cr}$ and 
$C_{ca}$. The control parameters are 
$T^+=50$ and 
$J^+=50$. The control is activated at 
$t=250$ and the circles denote the start of averaging. The horizontal dashed and solid lines indicate the mean uncontrolled and controlled coefficients, respectively.
Table 1 lists the values of the drag coefficients for different amplitudes 
$A^+$, periods of oscillations 
$T^+$ and impulses in the range 
$J^+=1\text {--}50$. In uncontrolled conditions, the contribution of the pressure drag to the total drag is more than four times the contribution of the skin-friction drag.
Table 1. Control parameters and mean drag coefficients. In this table and in tables 2 and 3, the case highlighted in bold is studied in § 5.
Table 2 presents the computed percentage variations of the drag coefficients. Note that, while the partial skin-friction and pressure-drag coefficients in table 1 add to give the corresponding total coefficients, the partial reductions of the coefficients in table 2 do not because each reduction is defined with respect to the corresponding value in the uncontrolled case, as clear from the definitions (2.9) and (2.10).
Table 2. Percentage variations of the drag coefficients. Positive values indicate drag reduction and negative values denote drag increase.
The control with small period 
$T=0.2$ (
$T^+=1.57$) reduces the total drag marginally as the impulse 
$J$ is too small. The small drag reduction is due to the small reduction of 
$C_p$, while both the cavity and the total skin friction increase slightly. For such a low period, a variation of 
$A$ of two orders of magnitude in the range 
$0.1\text {--}10$ (
$A^+=0.29\text {--}29$) does not cause any improvement. The drag increases slightly in some cases for 
$J^+=3.2$ and 
$4.4$. Both the skin-friction drag and the pressure drag decrease at higher impulse values, 
$J^+=15, 50$. The case with 
$T=10.1$ and 
$A=0.67$ (
$T^+=80$ and 
$J^+=50$) shows the maximum total drag-reduction margin, 
$\mathcal {R}=25\,\%$, the maximum skin-friction drag-reduction margin, 
$\mathcal {R}_f=37\,\%$, and the maximum pressure-drag reduction margin, 
$\mathcal {R}_p=22\,\%$.
Figure 9 shows the total drag-reduction margin as a function of the actuation period at the impulse values 
$J^+=15$ and 50. For a constant 
$T^+$, the drag-reduction margin is larger as the impulse increases, while, for a constant 
$J^+$, it increases with the period up to a maximum at 
$T^+=80$. The existence of an optimum period of oscillation is consistent with similar control strategy for skin-friction drag reduction, such as the oscillating-wall technique in smooth channels (Ricco & Quadrio Reference Ricco and Quadrio2008). The optimal period of wall oscillation in the flat-wall case is 
$T^+=100$ (Touber & Leschziner Reference Touber and Leschziner2012), which gives maximum drag-reduction margins of 35 %–40 %, thus in good agreement with the maximum found in the channel with bars altered by a spanwise pressure gradient.
Figure 9. Drag-reduction margin as a function of the oscillating period. The circles indicate the numerical results and the solid lines are the best fits at constant 
$J^+=A^+ T^+/{\rm \pi}$.
Figure 10 reveals that the reduction in viscous drag along the cavity contributes to most of the total viscous reduction because the skin friction on the bar crest is negligible. It is also shown that the total friction drag reduction 
$\mathcal {R}_f$ is almost double the total pressure-drag reduction 
$\mathcal {R}_p$. The viscous drag reduction reaches a maximum of approximately 
$37\,\%$ for 
$T^+=50\text {--}100$, while the pressure-drag reduction reaches its maximum value of approximately 
$22\,\%$ for 
$T^+=80\text {--}100$. The amount of skin-friction drag reduction is comparable with the maximum reduction brought about by spanwise wall oscillations in a fully developed channel flow for 
$T^+=100$, and similar wall-friction Reynolds numbers and maximum velocity amplitude of oscillation (Quadrio & Ricco Reference Quadrio and Ricco2004). While the scaled forces are additive, i.e. 
$-\bar {\varPi }_x L_y=C_f+C_p$, the reductions of viscous drag, 
$\mathcal {R}_f$, and of the pressure drag, 
$\mathcal {R}_p$, lead to 
$\mathcal {R}$ via
\begin{equation} \mathcal{R} = 100 - \frac{1}{\bar{C}+1}\left[ 100 - \mathcal{R}_f + \bar{C} \left(100-\mathcal{R}_p \right) \right], \end{equation}
where 
$\bar {C}=C_p/C_f=4.3$.
Figure 10. Drag-reduction margins as functions of the oscillating period for an impulse 
$J^+=50$.
4.3. Power spent and net power saved
Table 3 gives the values of the drag-reduction margin, 
$\mathcal {R}(\%)$ (which coincides with the power saved), the power spent to activate the control, 
$\mathcal {P}_z(\%)$, and the laminar power 
$\mathcal {P}_{z,l}(\%)$ for different parameters 
$T$ and 
$A$. Very marginal savings in terms of 
$\mathcal {P}_{net}(\%)$ are found when the drag-reduction effect is small and these values are within the uncertainty of the computed 
$\mathcal {R}(\%)$. For the cases with large drag reduction, the decrease in pressure gradient along the streamwise direction comes at the expense of a high consumption of energy, which leads to no savings in terms of 
$\mathcal {P}_{net}(\%)$. This finding opens new avenues of research to improve the net energy gain with a view to rendering this technique competitive with vortex generators which, at this stage, are more effective and guarantee estimated net drag reductions of the order of 
$10\,\%$ (Calarese et al. Reference Calarese, Crisler and Gustafson1985).
Table 3. Power saved 
$\mathcal {P}_z(\%)$, net power 
$\mathcal {P}_{net}(\%)$ and laminar power 
$\mathcal {P}_{z,l}(\%)$ as functions of the control parameters 
$A$ and 
$T$.
Figure 11 shows the values of 
$\mathcal {P}_z$ obtained via the direct numerical simulations versus the laminar prediction 
$\mathcal {P}_{z,l}$. The laminar power 
$\mathcal {P}_{z,l}$ gives an excellent approximation of 
$\mathcal {P}_z$ for all the values of 
$A$ and 
$T$ considered. The laminar solution found in § 3 is therefore able to predict the turbulent power spent very well.
Figure 11. Turbulent power spent 
$\mathcal {P}_z(\%)$ versus its laminar prediction 
$\mathcal {P}_{z,l}(\%)$. The solid line indicates the line 
$\mathcal {P}_z(\%)=\mathcal {P}_{z,l}(\%)$.
5. Turbulent-flow physics
In this section, instantaneous flow fields and flow statistics are presented for the controlled case with 
$A=1.08$ and 
$T=6.33$ (
$J^+=50$ and 
$T^+=50$), which features a large total drag-reduction margin, 
$\mathcal {R}=23\,\%$.
5.1. Controlled instantaneous fields and average flow
Instantaneous flow fields are shown in figures 12–14 to illustrate the effects of the square bar and of the spanwise-forcing control on the flow dynamics. Figure 12 shows isosurfaces of the invariant 
$\lambda _2$, proposed by Jeong & Hussain (Reference Jeong and Hussain1995) to visualize vortical structures. Vortices are formed on top and downstream of the bar, especially where the flow fully separates, and are mostly attenuated upstream of the bar. Figures 13 and 14 depict instantaneous streamwise and spanwise velocity fields in the central 
$x$–
$y$ plane at 
$z={\rm \pi} /2$. The controlled streamwise velocity is much more uniform than in the reference case, both upstream and downstream of the bar. In the controlled case, the spanwise velocity is positive everywhere in the visualization plane and much larger than in the reference case because of the large spanwise pressure gradient at this phase of the oscillation.
Figure 12. Isosurfaces of 
$\lambda _2=-1.0$ in the channel for the uncontrolled case (a) and the controlled case (b). The isosurface contours are coloured by the streamwise velocity.
Figure 13. Contour plots of instantaneous streamline velocity 
$u$ in the channel for the uncontrolled case (a) and the controlled case (b).
Figure 14. Contour plots of instantaneous spanwise velocity, 
$w$, in the channel for uncontrolled (a) and controlled case (b).
The contours of the mean streamwise velocity, shown in figure 15, confirm that the main flow features around the bar are qualitatively unvaried in the controlled case. The three recirculation regions are still present; the large recirculation area behind the bar is slightly compressed towards the wall and the one upstream of the bar is weakened. Figure 16 compares the mean streamlines in the uncontrolled and controlled cases. The forcing reduces the size of vortices B and C and brings them closer to the wall, thus rendering the flow more attached downstream of the bar. The reattachment point upstream of the vortex A displaces slightly and a small vortex forms upstream of vortex A with its core at 
$x=3.72$. After the bar, although the control shifts the second reattachment point downstream, this change only marginally affects the pressure-drag balance, which is instead mainly influenced by the dynamics in the whole separated region and by the flow modifications upstream of the bar. Similar changes of the separation area as an effect of the control were also found in previous separation drag studies (Minelli et al. Reference Minelli, Tokarev, Zhang, Liu, Chernoray, Basara and Krajnović2019; Banchetti et al. Reference Banchetti, Luchini and Quadrio2020).
Figure 15. Contour plots of time- and spanwise-averaged velocity 
$\langle \bar {u} \rangle$ in the channel for uncontrolled (a) and controlled case (b). The contour lines are for the negative velocity range of 
$\langle \bar {u} \rangle$ in 
$(-0.2,-0.15,-0.1,-0.05,-0.025)$.
Figure 16. Effect of the control on the separation areas. Dashed and solid lines indicate the time- and spanwise-averaged streamlines for the uncontrolled and controlled flow, respectively.
5.2. Distribution of pressure and skin-friction drag coefficients
The difference of the time- and spanwise-averaged pressure between the controlled case and the uncontrolled case is shown in the contour plot of figure 17. The pressure is unvaried in a region that extends from the bar upward toward the channel centreline. The pressure upstream of the bar decreases, while the pressure instead increases in the region downstream the bar with the maximum increase inside the area delimited by the isoline around the point 
$(x,y)=(4.75,0.1)$. This dual effect on the pressure distribution is the main cause of the reduced pressure difference on the two sides of the bar, which gives a lower overall pressure drag.
Figure 17. Colour map and contours of the difference of the mean pressure between the uncontrolled case and the controlled case. The positive values indicate a pressure increase and the negative values denote a pressure decrease. The isoline values are equispaced at intervals of 
$0.002$ in the range 
$[-0.02,0.02]$.
Figure 18 depicts the mean pressures in the controlled and uncontrolled cases (18a) and their difference (18b). The change of pressure is much smaller away from the bar, especially downstream of the large separation region, in the fully attached region for 
$s>6.5$, where 
$s$ is the streamwise coordinate running along the solid wall. The pressures are kept equal to zero at the further end of the domain, as discussed in § 2.3. Excellent agreement is found between our data and the data behind the bar from the direct numerical simulations of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) (red circles), obtained for a flow over a squared bar with the same height as ours, a slightly different Reynolds number and a distance between bars which is approximately half of ours (
$h/L_x=19$). This match confirms that the recirculation area is largely unaffected once the bars are sufficiently spaced apart. The distribution of the pressure difference in the separated region is qualitatively similar to that found by Banchetti et al. (Reference Banchetti, Luchini and Quadrio2020) for a channel flow modified by a smooth bump and altered by wall travelling waves. In their flow, the total pressure sufficiently downstream of the obstacle recovers the negative constant slope behaviour proper of channel flows over smooth walls. In our case, figure 18(a) instead shows that the pressure further downstream of the bar tends to a constant value because the effect of the intense separation is still present. The same trend was found by Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003). This finding is expected because Banchetti et al.'s (Reference Banchetti, Luchini and Quadrio2020) flow is characterized by a mild separation being their obstacles smoother and thinner than our bars and more spaced apart along the streamwise direction.
Figure 18. (a) Time- and spanwise-averaged total pressure in the uncontrolled case (dashed line) and the controlled case (solid line) as a function of the coordinate 
$s$ running along the wall surface. The sides of the bar are shown in grey and the crest of the bar is marked in light blue. The red circles are data from the direct numerical simulations of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003). (b) Difference between the time- and spanwise-averaged pressure in the controlled and uncontrolled cases.
The graphs of figure 19 show that the wall-shear stress decreases along the major part of the cavity area where the flow is attached (
$x<3.5$ and 
$x>5.4$), while there is no significant change of skin friction upstream of the bar (
$3.5<x<3.9$). The wall-shear stress is mostly affected by the forcing when it is positive, i.e. where the flow is attached along the cavity, reaching a reduction of 
$31\,\%$. This value is consistent with the drag-reduction margin of 
$\mathcal {R}=35\,\%$ obtained by forcing a flat-wall turbulent channel flow by spanwise wall oscillations at a similar friction Reynolds number, period 
$T^+=50$ and maximum spanwise velocity of approximately 
$27u_\tau ^*$ (Quadrio & Ricco Reference Quadrio and Ricco2004). As for the pressure distribution, the red circles show that the skin-friction distribution in the recirculation area agrees very well with that of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003).
Figure 19. (a) Time- and spanwise-averaged wall-shear stress in the uncontrolled case (dashed line) and the controlled case (solid line). The red circles are data from the direct numerical simulations of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003). (b) Difference between the time- and spanwise-averaged wall-shear stress caused by the control. The prime indicates derivation with respect to 
$y$.
The negative skin friction in the region behind the bar is instead largely unperturbed. The wall-shear stress in the separation area B between 
$x=4.1$ and 
$x=4.4$, defined in figure 16, is unaffected by the forcing, although the near-wall flow is along the streamwise direction. This finding is likely to be due to the stress being small and the flow not being in the fully developed turbulent regime there. This result is consistent with the spanwise forcing not having an effect on a streamwise laminar flow (Ricco & Hicks Reference Ricco and Hicks2018). Figure 19(a) also shows that the reattachment points upstream of the bar and immediately downstream of the bar are unaffected by the control, while the reattachment point further downstream is shifted only marginally in the streamwise direction.
5.3. Turbulent-flow statistics
Figures 20 and 21 show the time- and spanwise-averaged profiles of the streamwise and wall-normal mean velocities along the channel for the uncontrolled and controlled cases. Close to the wall, the streamwise velocity is smaller in the controlled case resulting in the lower skin friction. The conservation of mass flow rate manifests itself in the wall-normal location increasing from the wall. Away from the wall, the streamwise velocity in the controlled case becomes first larger and then smaller than in the uncontrolled case as the channel centreline is approached when the flow is streamwise only at every wall-normal location. The opposite occurs when the wall-shear stress is instead negative. The decrease of pressure on the left-hand side of the bar could be due to the reduced dynamic pressure of the flow encountering the bar, due to the lower mean velocity, shown in the first three profiles of figure 20(a). The wall-normal velocity is unaffected by the control when the flow is fully attached, while it slightly increases in the recirculation region.
Figure 20. Profiles of streamwise mean velocity 
$\langle \bar {u} \rangle$ at different 
$x$ locations along the channel. The velocity profiles are shown up to the height of the bar (a) and to the channel centreline (b). In this figure and in figures 21 and 22, solid and dashed lines show profiles in the controlled and uncontrolled cases, respectively. The top horizontal axis denotes the velocity magnitude. The bars are not drawn to scale.
Figure 21. Profiles of the wall-normal mean velocity 
$\langle \bar {v} \rangle$ at different 
$x$ locations along the channel. The bar is not drawn to scale.
Figure 22. Mean streamwise velocity (a) and wall-normal velocity (b) profiles scaled in viscous units at different locations along the channel. The wall-friction velocity for each case, defined in § 4, is used for scaling.
Figure 22 shows the velocity profiles scaled in viscous units and at locations where the streamwise flow is attached. The wall-friction velocity of each case is used. The streamwise velocity profile shifts upwards in the controlled case, thus showing a thicker viscous layer and buffer region, a well-known feature of drag-reduced attached flows. The wall-normal velocity is instead only marginally affected.
Figure 23 depicts the spanwise- and phase-averaged velocity components at different streamwise locations and phase of oscillation 
$\tau$. The velocities 
$\left \langle \hat {u} \right \rangle$ and 
$\left \langle \hat {v} \right \rangle$ do not depend on 
$\tau$ and overlap with the time-averaged velocities 
$\left \langle \bar {u} \right \rangle$ and 
$\left \langle \bar {v} \right \rangle$ when the flow is fully attached and show only a marginal modulation with the phase in the recirculation region downstream of the bar. This result is in accordance with the behaviour of turbulent channel flows with spanwise oscillating walls (Jung et al. Reference Jung, Mangiavacchi and Akhavan1992) and validates the assumption of the steadiness of 
$\left \langle \hat {u} \right \rangle$ and 
$\left \langle \hat {v} \right \rangle$, adopted for the laminar-flow analysis of § 3. The velocity 
$\left \langle \hat {v} \right \rangle$ is smaller than 
$\left \langle \hat {u} \right \rangle$ for every 
$x$ and 
$\tau$, and is negligible outside the recirculation areas (
$x=7$). The spanwise profiles 
$\left \langle \hat {w} \right \rangle$ at 
$\tau =0$ and 
$\tau =T/4$ are 
$90^{\circ }$ out of phase with the 
$\left \langle \hat {w} \right \rangle$ profiles at 
$\tau =T/2$ and 
$\tau =3/4 T$. We also note that the spanwise velocity 
$\left \langle \hat {w} \right \rangle$ is comparable with 
$\left \langle \bar {u} \right \rangle$ in the bulk of the channel.
Figure 23. Profiles of spanwise- and phase-averaged velocities 
$\left \langle \hat {u}\right \rangle$, 
$\left \langle \hat {v}\right \rangle$ and 
$\left \langle \hat {w}\right \rangle$ at different 
$x$ locations. Thin solid lines with symbols are spanwise- and phase-averaged profiles at different phases (circle: 
$\tau =0$; square: 
$\tau =T/4$; diamond: 
$\tau = T/2$ and triangle: 
$\tau =3T/4$). The lines indicate the time- and spanwise-averaged velocities 
$\left \langle \bar {u}\right \rangle$, 
$\left \langle \bar {v}\right \rangle$ and 
$\left \langle \bar {w}\right \rangle$ for the controlled case (solid line) and the uncontrolled case (dashed line). The bars are not drawn to scale.
Figure 24 displays the profiles of the viscous-scaled phase-averaged velocities 
$\left \langle \hat {u} \right \rangle ^+$ and 
$\left \langle \hat {v} \right \rangle ^+$ at two locations where the flow is attached. Both velocity components fluctuate more downstream than upstream of the bar. The velocity 
$\left \langle \hat {v} \right \rangle$ increases up to an absolute maximum and then decreases monotonically to zero as the channel centreline is approached.
Figure 24. Profiles of spanwise- and phase-averaged velocities 
$\left \langle \hat {u}\right \rangle$, 
$\left \langle \hat {v}\right \rangle$ scaled in viscous units at two 
$x$ locations. Symbols and lines are the same in figure 23.
Figure 25 depicts the profiles of 
$\left \langle \hat {w} \right \rangle$ and 
$w_l$ for 
$\tau =0, T/4, T/2,3/4 T$ at different 
$x$ stations. The laminar solution 
$w_l$ agrees well with 
$\left \langle \hat {w} \right \rangle$ outside the separation area (
$x=2$ and 
$x=7$) for 
$\tau =T/4, 3/4 T$. At 
$\tau =0, T/2$, a large discrepancy occurs between the laminar and the turbulent profiles. The laminar 
$w_l$ tends to be null in the outer region of 
$y^+>10$, contrary to the finite values of the turbulent 
$\left \langle \hat {w} \right \rangle$. The differences between the laminar and turbulent flows are due to the transverse Reynolds stresses 
$-\partial \left \langle \widehat {u_t w_t} \right \rangle /\partial x$ and 
$-\partial \left \langle \widehat {v_t w_t} \right \rangle /\partial y$ in (3.4). Overall, the turbulent velocity 
$\left \langle \hat {w} \right \rangle$ shows a satisfactory agreement with the laminar 
$w_l$, proving that the analysis of § 3 gives a relevant physical model for the transverse oscillating mean flow obtained from the direct numerical simulations. This result is confirmed by the excellent laminar prediction of the spanwise power 
$\mathcal {P}_z$, shown in figure 11.
Figure 25. Comparisons between the profiles of 
$\left \langle \hat {w} \right \rangle$ (solid lines) and the laminar 
$w_l$ obtained via (3.5) the last term on the right-hand side (dashed lines). Profiles are at 
$\tau =0$ (circles), 
$\tau =T/4$ (squares), 
$\tau = T/2$ (diamonds) and 
$\tau =3T/4$ (triangles); (a) 
$x = 2$, (b) 
$x = 5$, (c) 
$x=6$, (d) 
$x=7$.
The turbulence intensities, 
$\left \langle \widehat {u_t u_t}\right \rangle$, 
$\left \langle \widehat {v_t v_t}\right \rangle$, 
$\left \langle \widehat {w_t w_t}\right \rangle$ and 
$\left \langle \overline {u_t u_t}\right \rangle$, 
$\left \langle \overline {v_t v_t}\right \rangle$, 
$\left \langle \overline {w_t w_t}\right \rangle$, are shown in figure 26. They are highest at the 
$x$ locations corresponding to the bar sides where peaks of 
$\left \langle \widehat {u_t u_t}\right \rangle$ and 
$\left \langle \overline {u_t u_t}\right \rangle$ reach values two and four times higher than the values far from the separation regions (
$x=3$ and 
$7$). The profiles of 
$\left \langle \widehat {u_t u_t}\right \rangle$ and 
$\left \langle \overline {u_t u_t}\right \rangle$ show two local maxima far from the bar with the global maxima occurring at 
$y=0.1$ for 
$x=3$ and 
$y=0.5$ for 
$x=7$. Similar to the mean-flow quantities, the phase-averaged quantities fluctuate more downstream (
$4 < x < 7$) than upstream of the bar. The streamwise and spanwise correlations are more dominant than the wall-normal correlation in both the controlled and uncontrolled cases. Both 
$\left \langle \widehat {w_t w_t}\right \rangle$ and 
$\left \langle \overline {w_t w_t}\right \rangle$ show a peak around the bar height 
$y=0.2$ for every 
$\tau$. The quantities are all reduced for any 
$\tau$ and at any location along the channel. Figure 27 confirms that the turbulent kinetic energy 
$0.5(\left \langle \overline {u_t u_t}\right \rangle +\left \langle \overline {v_t v_t}\right \rangle +\left \langle \overline {w_t w_t}\right \rangle )$ decreases in the detached region downstream of the bar.
Figure 26. Profiles of 
$\left \langle \widehat {u_t u_t}\right \rangle$, 
$\left \langle \widehat {v_t v_t}\right \rangle$, 
$\left \langle \widehat {w_t w_t}\right \rangle$ at different 
$x$ stations and 
$\tau$ (thin solid lines with symbols). Profiles are at 
$\tau =0$ (circle), 
$\tau =T/4$ (square), 
$\tau = T/2$ (diamond) and 
$\tau =3T/4$ (triangle). The phase-averaged quantities are compared with time-averaged flow statistics 
$\left \langle \overline {u_t u_t}\right \rangle$, 
$\left \langle \overline {v_t v_t}\right \rangle$, 
$\left \langle \overline {w_t w_t}\right \rangle$ for the control (solid lines) and uncontrolled (dashed lines) case. The bars are not drawn to scale.
Figure 27. Contour plots of turbulent kinetic energy in the region around the bar for uncontrolled (a) and controlled (b) cases.
Figure 28 compares the Reynolds stresses of the controlled and uncontrolled flows. The Reynolds stress 
$\left \langle \widehat {u_t v_t}\right \rangle$ is dominant over 
$\left \langle \widehat {u_t w_t}\right \rangle$ and 
$\left \langle \widehat {v_t w_t}\right \rangle$ and decreases with respect to the uncontrolled case for almost every 
$y$ at any 
$x$ and 
$\tau$, similar to Jung et al. (Reference Jung, Mangiavacchi and Akhavan1992). The attenuating effect on 
$\left \langle \widehat {u_t v_t}\right \rangle$ is most intense downstream of the bar. The stresses 
$\left \langle \overline {v_t w_t}\right \rangle$ and 
$\left \langle \overline {u_t w_t}\right \rangle$ remain smaller than 
$\left \langle \widehat {u_t v_t}\right \rangle$ in the controlled case. (Note that 
$\left \langle \overline {u_t w_t}\right \rangle$ vanishes in the averaged 
$z$-momentum equation of a flat-wall turbulent channel flow because of streamwise averaging, while in the channel with bars the streamwise derivative of this term is present in the balance, as shown in (3.4).) The small amplitudes of these Reynolds stresses give support to the assumption adopted in the laminar analysis of (3.4), where these quantities are neglected.
Figure 28. Profiles of 
$\left \langle \widehat {u_t v_t}\right \rangle$, 
$\left \langle \widehat {v_t w_t}\right \rangle$, 
$\left \langle \widehat {w_t u_t}\right \rangle$ at different 
$x$ stations and 
$\tau$ (thin solid lines with symbols). Lines and symbols are the same as in figure 26. The bars are not drawn to scale.
5.4. Role of the Reynolds stresses between bars
In order to gain further insight into the role of the Reynolds stresses in the controlled channel flow, we consider a momentum balance relating the total drag to the Reynolds stresses 
$-\left \langle \widehat {u_tv_t} \right \rangle$ along the horizontal line 
$l_m$ connecting the crests of two consecutive bars. The sum of the drag coefficients 
$C_p(\tau )+C_f(\tau )$ for the controlled flow is related to the product 
$\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$, the Reynolds stress 
$-\left \langle \widehat {u_tv_t} \right \rangle$ and the friction term 
$Re_p^{-1}\, \partial \left \langle \widehat {u} \right \rangle /\partial y|_h$ at 
$y=h$ via
\begin{equation} C_p(\tau)+C_f(\tau) \approx -\frac{L_y}{L_x}\int_{l_m} \left \langle \widehat{u} \right\rangle\left \langle \widehat{v} \right\rangle \,{\textrm{d}x}-\frac{L_y}{L_x}\int_{l_m} \left \langle \widehat{u_tv_t} \right\rangle \,{\textrm{d}x} +\frac{L_y}{L_x Re_p}\int_{l_m} \frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \,{\textrm{d}x}, \end{equation}where the phase dependence of the quantities is kept herein.
Appendix E presents the derivation of (5.1). Note that in (5.1) the total drag 
$C_p(\tau )+C_f(\tau )$ is only approximated by the terms on the right-hand side because 
$C_{cr}$, which is much smaller than 
$C_{ca}$, is neglected. Equation (5.1) differs from (5.3) in Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) as they did not include the term due to 
$\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$ and we use the time-ensemble average instead of the time average. Figure 29 depicts the integrands of the terms in (5.1) for the controlled and uncontrolled flows. The Reynolds stresses 
$-\left \langle \widehat {u_tv_t} \right \rangle$ and shear stresses 
$\partial \left \langle \widehat {u} \right \rangle /\partial y$ are positive at every 
$x$ along the cavity mouth, giving a positive contribution to the overall drag for both flows. The term involving 
$\partial \left \langle \widehat {u} \right \rangle /\partial y$ is, however, negligible compared to the other terms. The term 
$-\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$ is positive and streamwise oscillatory downstream of the bar in the controlled case, while it is negative upstream of the bar (
$2<x<4$) because the wall-normal velocity 
$\left \langle \widehat {v} \right \rangle$ is large and positive. Due to conservation of mass, this upward velocity balances the decreasing streamwise mean velocity caused by the presence of the bar. Figure 29 shows that the reduction of the total drag mainly occurs because the control attenuates the Reynolds stresses 
$-\left \langle \widehat {u_tv_t} \right \rangle$, while the control effect of the time average of 
$-\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$ and 
$\partial \left \langle \widehat {u} \right \rangle /\partial y$ is marginal. The attenuation of 
$-\left \langle \widehat {u_tv_t} \right \rangle$ occurs not only in the region where the flow is attached (
$x<3.5$ and 
$x>5.5$), but also along 
$l_m$ above the vortex A (
$3.5<x<4$). Good agreement is obtained behind the bar between our 
$\partial \left \langle \widehat {u} \right \rangle /\partial y$ data along the channel mouth and the data from the direct numerical simulations of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003) (red circles), while the comparison between the 
$-\left \langle \widehat {u_tv_t} \right \rangle$ data is only satisfactory downstream of the peak. The mismatch of the maximum values could be due to the large uncertainty caused by the large velocity fluctuations in the recirculation region.
Figure 29. Comparisons between the profiles of 
$-\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$, 
$-\left \langle \widehat {u_tv_t} \right \rangle$, 
$Re_p^{-1} \partial \left \langle \widehat {u} \right \rangle /\partial y |_h$ of the controlled flow at phases 
$\tau =0, T/4, T/2, 3/4 T$ (thin solid lines) and the uncontrolled terms 
$-\left \langle \bar {u} \right \rangle \left \langle \bar {v} \right \rangle$, 
$-\left \langle {\overline {u_tv_t}}\right \rangle$, 
$Re_p^{-1} \partial \left \langle \bar {u} \right \rangle /\partial y |_h$ (thick solid lines). The red circles are data from the direct numerical simulations of Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003).
We can interpret the reduction of the Reynolds stresses as the effect of the oscillation on the bulk turbulent flow that in turn manifests itself at the wall surfaces as the attenuation of the pressure difference on the vertical surfaces of the bars and of the wall-shear stress over the cavity wall. Further research work is needed to establish the link between the bulk Reynolds-stress reduction and the pressure changes at the bar sides.
6. Conclusions
A drag-reduction study on a channel flow with square bars located on each wall has been conducted using direct numerical simulations. The flow is controlled via a spanwise oscillating pressure gradient and a maximum total drag-reduction margin of 
$25\,\%$ is achieved. We have shown that a reduction of both pressure and skin-friction drag is possible via an active method, contrary to vortex generators and jets that may energize the fluid and increase the skin-friction drag to reduce the form drag. The drag reduction has been linked to a reduction of the Reynolds stress along the dividing streamline connecting the bar crest and the reattachment point along the cavity. The skin-friction drag is mostly reduced because the wall-shear stress of the attached flow along the cavity wall is decreased under the effect of the control, whereas the control has no significant effect on the friction drag of the separated flow and of the flow over the bar crest because the wall-shear stress is low or negative there. The pressure drag decreases because the pressure decreases in front of the bar and increases behind the bar. The lower pressure on the front of the bar is caused by the lower dynamical pressure of the mean flow approaching the bar, which is due to the reduced wall-shear stress along the cavity upstream of the bar.
The power to control the flow has been quantified and compared with the power saved thanks to drag reduction to compute the net energy saved. For the amplitudes and periods considered, the solution of the related spanwise laminar flow gives an excellent prediction of the turbulent power spent. No net power saved has been computed for the cases that feature the largest drag-reduction margins. The small net power savings, obtained for small drag-reduction margins, are within the uncertainty error of the computations, as discussed in appendix A.
Further studies are necessary to study the mechanism causing drag reduction by, for example, an analysis analogous to that of Ricco et al. (Reference Ricco, Ottonelli, Hasegawa and Quadrio2012) based on the turbulent enstrophy and dissipation or by extending the useful identity of Fukagata, Iwamoto & Kasagi (Reference Fukagata, Iwamoto and Kasagi2002) to a channel with bars to investigate the role of the Reynolds stresses on the drag-reduction effect. It is also interesting to determine the effect of the Reynolds number on the drag-reduction performance and whether some combinations of parameters and of the forcing actuation can further reduce drag. This technique is promising because it is shown that reduction of both pressure and skin-friction drag is possible, paving the way for research into more effective types of separation control.
Future work should certainly be directed at improving the power balance to obtain a measurable net power saved. It is known that a travelling wave forcing is more efficient than spanwise oscillation (Quadrio, Ricco & Viotti Reference Quadrio, Ricco and Viotti2009) in reducing the skin-friction drag. A possibility to control separated flows could thus be a spanwise pressure gradient of the form 
$\varPi _z(x,t)=A\sin (k_x x-\omega t)$ and giving an effect similar to the one of Quadrio et al. (Reference Quadrio, Ricco and Viotti2009) in the attached-flow region. For improved efficiency, such an active technique could be confined in the proximity of the bars and combined with a wall-based passive control method along the cavity where the flow separates. An idea could be the method proposed by Heenan & Morrison (Reference Heenan and Morrison1998), who were able to reduce the separation drag by 
$9\,\%$ in a forward-facing step by making permeable the cavity wall in front of the step and obtaining a decrease of the root mean square of the pressure fluctuations.
Acknowledgements
We would like to thank the Department of Mechanical Engineering at the University of Sheffield and the Agency for Science, Technology and Research (A*STAR) in Singapore for funding this research. We would like to express our gratitude towards the A*STAR Computational Resource Centre for their kind assistance during the use of their high performance computing resource. G.P. chose the references, conceived and wrote §§ 2, 3, 5.4 and the appendixes, modified Nek5000 for this problem and generated half of the data at low 
$T$ in tables 1–3 and figures 1–5. Part of this work was presented by G.P. at the European Drag Reduction and Flow Control Meeting, Bad Herrenalb, Germany in March 2019. We would like to thank Professors S. Chernyshenko, S. He and P. Luchini for useful discussions. We are also indebted to Professor P. Fischer for the advice concerning the mesh creation and for the useful information on running Nek5000.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Uncertainty of the drag components for the uncontrolled flow
The sensitivity of the drag coefficients on the mesh resolution and the time step, for the uncontrolled case, is examined in table 4. Time accuracy is checked by comparing the two results obtained using the same grid with 
$l_x=6$, but with two different time steps, 
$\Delta t=1.69 \times 10^{-4}$ and 
$\Delta t=8.45\times 10^{-5}$. A percentage variation of 
$0.4\,\%$ in the total drag 
$C_p+C_f$ is obtained and variations of all the other coefficients are less than 
$1\,\%$, indicating that the time step used is adequate. The accuracy of the drag with space resolution is evaluated comparing results with three different grids. Each grid has the same number of spectral elements, but a different number 
$l_x$ of GLL points is employed along each direction because the SEM order of convergence is spectral with 
$l_x$. The computation with 
$l_x=8$ gives variations of less than 
$1\,\%$ if compared to the one with 
$l_x=6$ except for 
$C_{cr}$ which, however, represents a small contribution to the total drag. The most resolved case (
$l_x=10$) gives the highest percentage variation of the total drag (
$1\,\%$) with respect to the less resolved run (
$l_x=6$). These tests prove that the computation with 
$l_x=6$ and a time step 
$\Delta t=1.69 \times 10^{-4}$ guarantees accurate values of the drag coefficients and therefore this configuration is adopted for all the simulations with control.
Table 4. Space and time resolution checks for the drag components.
Appendix B. Drag balance and power for driving the flow along 
$x$
In this appendix, we find the expression for the power spent to drive the flow along the streamwise direction. We first consider the volume integral of the streamwise momentum equation
\begin{equation} \frac{\partial }{\partial t}\int_{V} u \, \mathrm{d} V +\int_{S} u \boldsymbol{u} \cdot \boldsymbol{n}\, \mathrm{d} S =-\int_{S} p\, \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S +\frac{1}{Re_p}\int_{S} \, \cdot \boldsymbol{n}\, \mathrm{d} S \boldsymbol{\nabla} \boldsymbol{u}, \end{equation}
where 
$S$ is the boundary surface of 
$V$, 
$\boldsymbol {u}$ is the velocity vector and 
$\boldsymbol {n}$ the unit vector pointing out of the control volume. The first integral in (B1) is null because the simulations are at constant mass flow rate. The boundary surface 
$S$ is composed of the opposite wall surfaces 
$S_w$, the plane surfaces 
${S_x}_i$ and 
${S_x}_f$ with normal vectors 
$\boldsymbol {i}$ and 
$-\boldsymbol {i}$ and the plane surfaces 
${S_z}_i$ and 
${S_z}_f$ of unit vectors 
$\boldsymbol {k}$ and 
$-\boldsymbol {k}$, as shown in figure 30. The second integral in (B1) is null because 
$\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {n}=0$ on 
$S_w$ and because the sum of the integrals over surfaces of opposite unit vectors is zero due to the periodic boundary conditions.
Figure 30. Schematic of the fluid volume 
$V$. The unit normal vectors to each surface are represented by arrows with full tip. The bars are not drawn to scale.
The pressure term in (B1) is obtained by defining the vector 
$\boldsymbol {p}=[p;0;0]$ and applying the divergence theorem, i.e.
\begin{equation} -\int_{V} \frac{\partial p}{\partial x} \, \mathrm{d} V =-\int_{V} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{p} \, \mathrm{d} V =-\int_{S} \boldsymbol{p} \cdot \boldsymbol{n}\, \mathrm{d} S =-\int_{S} p \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S. \end{equation}
The surfaces 
$S_l$, 
$S_r$, 
${S_x}_f$ and 
${S_x}_i$ are the only ones with 
$\boldsymbol {i}\boldsymbol {\cdot } \boldsymbol {n}\ne 0$ and the pressure integral in (B1) is
\begin{align} -\int_{S} p\, \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S &=-\left( \int_{S_l} p \, \mathrm{d} S -\int_{S_r} p \, \mathrm{d} S +\int_{{S_x}_f} p \, \mathrm{d} S -\int_{{S_x}_i} p \, \mathrm{d} S \right) \nonumber\\ &=-\left( C_pL_x L_z+\int_{0}^{L_y} \int_{0}^{L_z} p|_{L_x}\, \mathrm{d} y\, \mathrm{d} z-\int_{0}^{L_y} \int_{0}^{L_z} p|_0\, \mathrm{d} y\, \mathrm{d} z\right) \nonumber\\ &=-C_p(t)L_x L_z-\varPi_x(t)L_x L_y L_z, \end{align}
where the definitions 
$p|_{L_x}=\phi |_{L_x}+\varPi _z z$, 
$p|_0=\phi |_0+\varPi _x L_x+\varPi _z z$ and 
$\phi |_{L_x}=\phi |_0$ have been used. The contributions of 
${S_x}_i$, 
${S_x}_f$, 
${S_z}_i$, 
${S_z}_f$ to the last term in (B1) is zero, similar to the term 
$\int _{S} u \boldsymbol {u} \cdot \boldsymbol {n}\, \mathrm {d} S$ and
\begin{align} \frac{1}{Re_p}\int_{S} \boldsymbol{\nabla} \boldsymbol{u} \cdot \boldsymbol{n}\, \mathrm{d} S &=\frac{1}{Re_p}\int_{S_{ca}} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S +\frac{1}{Re_p}\int_{S_{cr}} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S \nonumber\\ &\quad + \frac{1}{Re_p}\int_{S_l} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S +\frac{1}{Re_p}\int_{S_r} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S \nonumber\\ &= \frac{1}{Re_p}\int_{S_{ca}} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S +\frac{1}{Re_p}\int_{S_{cr}} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S \nonumber\\ &= -C_{ca}(t)L_x L_z-C_{cr}(t)L_x L_z=-C_f(t) L_x L_z, \end{align}
where the minus sign is due to the orientation of 
$\boldsymbol {n}$, 
$S_{ca}$ and 
$S_{cr}$ are the cavity and crest surfaces, respectively. The integrals over 
$S_l$ and 
$S_r$ in (B4) are null because of continuity and the no-slip condition for 
$v$ and 
$w$, i.e. over the surface 
$S_l$
\begin{equation} \frac{1}{Re_p}\int_{S_l} \boldsymbol{\nabla} u \boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} S =\frac{1}{Re_p}\int_{0}^{h} \int_{0}^{L_z} \frac{\partial u}{\partial x}\, \mathrm{d} y\, \mathrm{d} z=-\frac{1}{Re_p}\int_{0}^{h} \int_{0}^{L_z} \left( \frac{\partial v}{\partial y}+\frac{\partial w}{\partial z} \right) \, \mathrm{d} y\, \mathrm{d} z=0. \end{equation}By substituting (B4) and (B3) into (B1) we find
\begin{equation} C_p(t)+C_f(t)=-\varPi_x(t) L_y \end{equation}
or its time-averaged version 
$C_p+C_f=-\bar {\varPi }_x L_y$.
We consider the first term of 1.108 in Hinze (Reference Hinze1975) averaged over 
$V$
\begin{equation} -\frac{1}{V}\int_{V} \frac{\partial up}{\partial x} \, \mathrm{d} V =-\frac{1}{V} \int_{S} up \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S. \end{equation}
As the surfaces 
$S_l$, 
$S_r$, 
${S_x}_f$ and 
${S_x}_i$ are the only ones with 
$\boldsymbol {i}\boldsymbol {\cdot } \boldsymbol {n}\ne 0$, we obtain
\begin{align} \int_{S} up \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S &=\int_{S_l} up \, \mathrm{d} S -\int_{S_r} up \, \mathrm{d} S +\int_{{S_x}_f} up \, \mathrm{d} S -\int_{{S_x}_i} up \, \mathrm{d} S \nonumber\\ &= \int_{{S_x}_f} up \, \mathrm{d} S -\int_{{S_x}_i} up \, \mathrm{d} S =\varPi_x(t) L_x \int_{0}^{L_y} \int_{0}^{L_z} u|_{L_x}\, \mathrm{d} y\, \mathrm{d} z, \end{align}
because of the no-penetration condition on 
$S_l$ and 
$S_r$, 
$p|_{L_x}=\phi |_{L_x}+\varPi _z z$, 
$p|_0=\phi |_0+\varPi _x L_x + \varPi _z z$ and 
$\phi |_{L_x}=\phi |_0$. The continuity equation in integral form states that the mass flow rate is the same at any cross section 
$S_x(x)$ of normal 
$\pm \boldsymbol {i}$ and independent of 
$x$
\begin{equation} \int_{{S_x}(x)} u \, \mathrm{d} S =\int_{0}^{L_y} \int_{0}^{L_z} u|_{L_x}\, \mathrm{d} y\, \mathrm{d} z. \end{equation}
The volume-averaged streamwise velocity 
$U_b(t)$ can be written using (B9) as
\begin{align} U_b(t)&=\frac{1}{V} \int_{V} u \, \mathrm{d} V =\frac{1}{V} \int_0^{L_x} \mathrm{d} x \int_{{S_x}(x)} u \, \mathrm{d} S \nonumber\\ &= \frac{1}{V} \int_0^{L_x} \mathrm{d} x \int_{0}^{L_y} \int_{0}^{L_z} u|_{L_x}\, \mathrm{d} y\, \mathrm{d} z=\frac{L_x}{V} \int_{0}^{L_y} \int_{0}^{L_z} u|_{L_x}\, \mathrm{d} y\, \mathrm{d} z, \end{align}
where 
$U_b$ does not depend on 
$t$ because the mass flow rate is constant.
By using (B10) and (B8), the power 
${P}_x$ per unit of volume used to drive the flow along 
$x$ is the time average of (B7), i.e.
\begin{equation} {P}_x=-\frac{1}{t_f-t_i}\int_{t_i}^{t_f} \frac{1}{V} \int_{S} up \boldsymbol{i} \cdot \boldsymbol{n}\, \mathrm{d} S \,\mathrm{d} t =-\bar{\varPi}_x U_b. \end{equation}Appendix C. Numerical solution of the laminar-flow equation
Equation (3.12) is discretized by using a centred finite-difference scheme
\begin{equation} \tilde{w}_{a,b}=C_{a,b}\,\tilde{w}_{a,b-1}+D_{a,b}\,\tilde{w}_{a-1,b}+F_{a,b}\,\tilde{w}_{a+1,b}+G_{a,b}\,\tilde{w}_{a,b+1}, \end{equation}
where 
$a$ and 
$b$ are the indices to the grid points and 
$C_{a,b}$, 
$D_{a,b}$, 
$F_{a,b}$ and 
$G_{a,b}$ become
\begin{gather} C_{a,b}=\dfrac{\dfrac{1}{\Delta y^2}+\dfrac{Re_p}{2\Delta y}\, \left \langle \hat{v}\right\rangle_{a,b} }{\textrm{i}\,\omega\,Re_p+\dfrac{2}{\Delta x^2}+\dfrac{2}{\Delta y^2}}, \quad D_{a,b}=\dfrac{\dfrac{1}{\Delta x^2}+\dfrac{Re_p}{2\Delta x}\, \left \langle \hat{u}\right\rangle_{a,b} }{\textrm{i}\,\omega\,Re_p+\dfrac{2}{\Delta x^2}+\dfrac{2}{\Delta y^2}}, \end{gather}
\begin{gather} F_{a,b}=\dfrac{\dfrac{1}{\Delta x^2}-\dfrac{Re_p}{2\Delta x}\, \left \langle \hat{u}\right\rangle_{a,b} }{\textrm{i}\,\omega\,Re_p+\dfrac{2}{\Delta x^2}+\dfrac{2}{\Delta y^2}}, \quad G_{a,b}=\dfrac{\dfrac{1}{\Delta y^2}-\dfrac{Re_p}{2\Delta y}\, \left \langle \hat{v}\right\rangle_{a,b} }{\textrm{i}\,\omega\,Re_p+\dfrac{2}{\Delta x^2}+\dfrac{2}{\Delta y^2}}. \end{gather}The boundary conditions (3.13a)–(3.13e) read
\begin{gather} \tilde{w}_{1,b}=\tilde{w}_{l,b}=-\frac{A \textrm{i}}{\omega}, \quad 1\leq b\leq m, \end{gather}
\begin{gather} \tilde{w}_{a,m}=\tilde{w}_{a,1}=-\frac{A \textrm{i}}{\omega}, \quad l\leq a\leq n_x, \end{gather}
\begin{gather} \tilde{w}_{l,b}=-\frac{A \textrm{i}}{\omega}, \quad 1\leq b\leq m, \end{gather}
\begin{gather} \tilde{w}_{a,n_y+1}=\tilde{w}_{a,n_y-1}, \quad 1\le a\le n_x, \end{gather}
\begin{gather} \tilde{w}_{0,b}=\tilde{w}_{n_x-1,b}, \quad 1\le b\le n_y, \quad \tilde{w}_{n_x,b}=\tilde{w}_{1,b}, \quad 1\le b\le n_y, \end{gather}
where 
$n_x$ and 
$n_y$ are the indices of the last points along 
$x$ and 
$y$. The scheme (C1) can be directly applied to the interior points. The (C1) is modified near the boundaries to take boundary conditions (C3a)–(C3e) into account
(i) ‘square’ scheme (
$\tilde {w}_{0,b}=\tilde {w}_{n_x-1,b}$)
(C4a)\begin{gather} \tilde{w}_{1,b}=C_{1,b} \tilde{w}_{1,b-1}+D_{1,b}\tilde{w}_{n_x-1,b}+F_{1,b} \tilde{w}_{2,b}+G_{1,b}\tilde{w}_{1,b+1}, \end{gather}(C4b)\begin{gather} m+1\le\ b \le n_y-1, \end{gather}(ii) ‘triangle’ scheme (
$\tilde {w}_{a,n_y+1}=\tilde {w}_{a,ny-1}$)
(C5a)\begin{gather} \tilde{w}_{a,n_y}=C_{a,n_y}\tilde{w}_{a,n_y-1}+D_{a,n_y}\tilde{w}_{a-1,n_y}+F_{a,n_y} \tilde{w}_{a+1,n_y}+G_{a,n_y}\tilde{w}_{a,ny-1}, \end{gather}(C5b)\begin{gather} 2 \le a\le n_x-2, \end{gather}(iii) ‘empty circle’ scheme (
$\tilde {w}_{n_x,b}=\tilde {w}_{1,b}$)
(C6a)\begin{gather} \tilde{w}_{n_x-1,b}=C_{n_x-1,b}\tilde{w}_{n_x-1,b-1}+D_{n_x-1,b}\tilde{w}_{n_x-2,b}+F_{n_x-1,b} \tilde{w}_{1,b}+G_{n_x-1,b}\tilde{w}_{n_x-1,b+1}, \end{gather}(C6b)\begin{gather} 3\le\ b \le n_y-1. \end{gather}
Points with indices 
$1,n_y$, 
$n_x-1,2$ and 
$n_x-1,n_y$ are not marked with any symbol because they satisfy these two boundary conditions
(i) point
$1,n_y$, ($\tilde {w}_{1,n_y+1}=\tilde {w}_{1,n_y-1}$, $\tilde {w}_{0,n_y}=\tilde {w}_{n_x-1,n_y}$),
(C7)\begin{equation} \tilde{w}_{1,n_y}=C_{1,n_y}\tilde{w}_{1,n_y-1}+D_{1,n_y}\tilde{w}_{n_x-1,n_y}+F_{1,n_y} \tilde{w}_{2,n_y}+G_{1,n_y}\tilde{w}_{1,n_y-1}, \end{equation}(ii) point
$n_x-1,2$, $\tilde {w}_{n_x-1,1}=-({A}/{\omega })\,\textrm {i}$, $\tilde {w}_{n_x,2}=\tilde {w}_{1,2}$,
(C8)\begin{equation} \tilde{w}_{n_x-1,2}=-C_{n_x-1,2}\frac{A}{\omega}+D_{n_x-1,2}\tilde{w}_{n_x-2,2}+F_{n_x-1,2} \tilde{w}_{1,2}+G_{n_x-1,2}\tilde{w}_{n_x-1,3} \end{equation}(iii) point
$n_x-1,n_y$, $\tilde {w}_{n_x-1,n_y+1}=\tilde {w}_{n_x-1,n_y-1}$, $\tilde {w}_{n_x,n_y}=\tilde {w}_{1,n_y}$,
(C9)\begin{align} \tilde{w}_{n_x-1,n_y}&=C_{n_x-1,n_y}\tilde{w}_{n_x-1,n_y-1}+D_{n_x-1,n_y} \tilde{w}_{n_x-2,n_y}+F_{n_x-1,n_y}\tilde{w}_{1,n_y}\nonumber\\ &\quad +G_{n_x-1,n_y}\,\tilde{w}_{n_x-1,n_y-1}. \end{align}
The equations are solved via a successive over-relaxation method
\begin{equation} \tilde{w}^{new}_{a,b}=\tilde{w}_{a,b}+\alpha r_{a,b}, \end{equation}where the residue is defined as
\begin{equation} r_{a,b}=C_{a,b}\tilde{w}_{a,b-1}+D_{a,b}\tilde{w}_{a-1,b}+F_{a,b} \tilde{w}_{a+1,b}+G_{a,b}\tilde{w}_{a,b+1}-\tilde{w}_{a,b}. \end{equation}
The convergence of the algorithm is guaranteed by monitoring 
$\| r \|= \sqrt {\sum _{a=1}^{n_x-1} \sum _{b=1}^{n_y} r_{a,b}\,r^{cc}_{a,b}}$, where 
$cc$ indicates the complex conjugate.
Appendix D. Relation with the laminar power of the oscillating wall
We consider a smooth laminar channel with spanwise oscillation of the walls
\begin{equation} \left. \begin{gathered} \frac{\partial w_s}{\partial \tau}=\frac{1}{Re_p} \frac{\partial^2 w_s}{\partial y^2}, \\ w_s(0,\tau)=w_s(2,\tau)= A \sin{(\omega \tau)}=\omega \int_0^{\tau} A \cos{(\omega \hat{\tau})} \mathrm{d} \hat{\tau}=\omega \int_0^{\tau} \varPi_z \mathrm{d} \hat{\tau}, \end{gathered} \right\} \end{equation}with the laminar power
\begin{equation} \mathcal{P}_s=\frac{1}{T Re_p}\int_0^T w_s \left.\frac{\partial w_s}{\partial y} \right\vert_{y=0} \,\mathrm{d} \tau. \end{equation}The system (D1) is transformed into the problem
\begin{equation} \left. \begin{gathered} \frac{\partial w_l}{\partial \tau}=\frac{1}{Re_p} \frac{\partial^2 w_l}{\partial y^2}-\varPi_z, \\ w_l(0,\tau)=w_l(2,\tau)=0, \end{gathered} \right\} \end{equation}by using the equation
\begin{equation} w_s=\omega \left( w_l+\int_0^{\tau} \varPi_z\, \mathrm{d} \hat{\tau} \right). \end{equation}
An expression for the spanwise volume-averaged velocity 
$W_{bs}=0.5\int _0^2 w_s \,\mathrm {d} y$ is obtained by averaging along the height of the channel (
$0 < y < 2$),
\begin{equation} \frac{\mathrm{d} W_{bs}}{\mathrm{d} \tau}=\frac{1}{2Re_p} \left( \left. \frac{\partial w_s}{\partial y} \right\vert_{y=2}-\left.\frac{\partial w_s}{\partial y} \right\vert_{y=0} \right) =-\frac{1}{Re_p}\left.\frac{\partial w_s}{\partial y} \right\vert_{y=0}. \end{equation}
The relationship between 
$W_{bs}$ and 
$W_{bl}$, defined in § 3, is obtained by averaging (D4) in the range 
$0 < y <2$,
\begin{equation} W_{bs}=\omega \left( W_{bl}+\int_0^{\tau} \varPi_z\, \mathrm{d} \hat{\tau} \right). \end{equation}By substituting (D6) into (D5) one finds
\begin{equation} \frac{1}{Re_p}\left.\frac{\partial w_s}{\partial y} \right\vert_{y=0}=-\omega\left( \frac{\mathrm{d} W_{bl}}{\mathrm{d} \tau}+ \varPi_z\right). \end{equation}
\begin{align} \mathcal{P}_s &= -\frac{\omega^2}{T}\int_0^T \left( \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \right) \left( \frac{\mathrm{d} W_{bl}}{\mathrm{d} \hat{\tau}}+ \varPi_z\right) \mathrm{d} \hat{\tau}\nonumber\\ &= - \frac{\omega^2}{T} \left( \int_0^T \frac{\mathrm{d} W_{bl}}{\mathrm{d} \hat{\tau}} \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \,\mathrm{d} \hat{\tau} + \int_0^T \varPi_z \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \,\mathrm{d} \hat{\tau} \right). \end{align}By integrating by parts the first integral in the square bracket of (D8) one finds
\begin{equation} \int_0^T \frac{\partial W_{bl}}{\partial \hat{\tau}} \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \,\mathrm{d} \hat{\tau}=\left[ W_{bl} \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \right] _0^T-\int_0^T W_{bl} \varPi_z \,\mathrm{d} \hat{\tau}=-\int_0^T W_{bl} \varPi_z \,\mathrm{d} \tau, \end{equation}
because 
$\int _0^T \varPi _z \,\mathrm {d} \tau =0$. The second integral in (D8) is null because
\begin{equation} \int_0^T \varPi_z \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau\,\mathrm{d}\hat{\tau} = \left[ \left( \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \right) ^2 \right] _0^T-\int_0^T \varPi_z \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \,\mathrm{d}\hat{\tau}, \end{equation}and
\begin{equation} \int_0^T \varPi_z\int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \,\mathrm{d}\hat{\tau} = \frac{1}{2}\left[ \left( \int_0^{\hat{\tau}} \varPi_z \,\mathrm{d} \tau \right) ^2 \right] _0^T=0. \end{equation}The substitution of (D9) and (D11) into (D8) gives
\begin{equation} \mathcal{P}_s=-\omega^2\left( -\frac{1}{T}\int_0^T W_{bl} \varPi_z \,\mathrm{d} \tau \right) =-\omega^2 \mathcal{P}_{z,l}=-4{\rm \pi}^2 T^{-2} \mathcal{P}_{z,l}. \end{equation}
Equation (D12) shows that the laminar power spent of a smooth oscillating wall 
$\mathcal {P}_s$ can be obtained by multiplying the laminar power spent of the oscillating pressure gradient 
$\mathcal {P}_{z,l}$ by a term proportional to 
$T^{-2}$.
Appendix E. Drag dependence on the Reynolds stresses
In this appendix, we study the mechanism by which the control decreases the time-ensemble-averaged total drag 
$C_p(\tau )+C_f(\tau )$ at different phases of the oscillation. The dependence of the quantities on the phase is kept herein. We consider a momentum balance relating the total drag to the Reynolds stress 
$-\left \langle \widehat {u_tv_t} \right \rangle$ along the horizontal line 
$l_m$ connecting the crests of two consecutive bars, i.e. the mouth of the cavity, as shown in figure 31. Inspired by the works of Gharib & Roshko (Reference Gharib and Roshko1987) and Leonardi et al. (Reference Leonardi, Orlandi, Smalley, Djenidi and Antonia2003), we consider the time-ensemble and 
$z$ averaged streamwise momentum equation
\begin{equation} \frac{\partial \left \langle \widehat{u u} \right\rangle}{\partial x}+ \frac{\partial \left \langle \widehat{u v} \right\rangle}{\partial y}=-\frac{\partial \left \langle \hat{p} \right\rangle}{\partial x}+\frac{1}{Re_p}\left( \frac{\partial^2 \left \langle \widehat{u} \right\rangle}{\partial x^2}+\frac{\partial^2 \left \langle \widehat{u} \right\rangle}{\partial y^2} \right). \end{equation}
Equation (E1) is steady because 
$\partial \left \langle \widehat {u} \right \rangle /\partial \tau =0$, as shown by figure 23. By defining the vector 
$\left \langle \widehat {uu} \right \rangle \boldsymbol {i}+\left \langle \widehat {uv} \right \rangle \boldsymbol {j}$ and noticing that 
$\partial \left \langle \hat {p} \right \rangle /\partial x=\boldsymbol {\nabla } \boldsymbol {\cdot }{(\left \langle \hat {p} \right \rangle \boldsymbol {i})}$, we write (E1) as
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot}{(\left \langle \widehat{uu} \right\rangle\boldsymbol{i}+\left \langle \widehat{uv} \right\rangle\boldsymbol{j})} = \boldsymbol{\nabla} \boldsymbol{\cdot}{\left( -\left \langle \hat{p} \right\rangle \boldsymbol{i}+\frac{1}{Re_p}\boldsymbol{\nabla} {\left \langle \widehat{u} \right\rangle}\right) }. \end{equation}
Integrating (E2) on a generic surface 
$S_z$ with boundary 
$\partial S_z$ lying in the plane 
$x-y$ and applying the divergence theorem leads to
\begin{equation} \int_{\partial S_z} (\left \langle \widehat{uu} \right\rangle\boldsymbol{i}+\left \langle \widehat{uv} \right\rangle\boldsymbol{j})\boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} l =-\int_{\partial S_z} \left \langle \hat{p} \right\rangle \boldsymbol{i}\boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} l +\frac{1}{Re_p}\int_{\partial S_z} \boldsymbol{\nabla} \left \langle \widehat{u} \right\rangle\boldsymbol{\cdot} \boldsymbol{n} \, \mathrm{d} l. \end{equation}
As the flow is periodic along 
$x$, it is possible to define 
$S_z$ as the rectangular area between two bars depicted in figure 31. Equation (E3) becomes
\begin{align} \left.\int_{x_b}^{x_a+L_x} \left \langle \widehat{uv} \right\rangle\right|_{h} \,{\textrm{d}x}&=\int_0^{h}\left( \left. \left \langle \widehat{p} \right\rangle\right|_{x_b}-\left.\left \langle \widehat{p} \right\rangle\right|_{x_a+L_x}\right) \,{\textrm{d}y}\nonumber\\ &\quad -\frac{1}{Re_p}\int_0^{h} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial x}\right|_{x_b} \,{\textrm{d}y}-\frac{1}{Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_0 \,{\textrm{d}x}\nonumber\\ &\quad + \frac{1}{Re_p}\int_0^{h} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial x}\right|_{x_a+L_x} \,{\textrm{d}y}+\frac{1}{Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_{h} \,{\textrm{d}x}. \end{align}
The terms in (E4) involving 
$\partial \left \langle \widehat {u} \right \rangle /\partial x$ are null because of continuity and the no-slip condition. Equation (E4) simplifies to
\begin{align} \left.\int_{x_b}^{x_a+L_x} \left \langle \widehat{uv} \right\rangle\right|_{h} \,{\textrm{d}x} &=\int_0^{h}\left(\left.\left \langle \widehat{p} \right\rangle\right|_{x_b}-\left.\left \langle \widehat{p} \right\rangle\right|_{x_a+L_x}\right) \,{\textrm{d}y}-\frac{1}{Re_p} \int_{x_b}^{x_a+L_x}\left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_0 \,{\textrm{d}x}\nonumber\\ &\quad + \frac{1}{Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_{h} \,{\textrm{d}x}. \end{align}
As the time-ensemble- and spanwise-averaged pressure is 
$\left \langle \widehat {p} \right \rangle =\left \langle \widehat {\phi } \right \rangle +\hat {\varPi }_x (x-L_x)+0.5\hat {\varPi }_z L_z$ and 
$\left \langle \widehat {\phi } \right \rangle =0$, we obtain
\begin{equation} \left.\left \langle \widehat{p} \right\rangle\right|_{x_b}-\left.\left \langle \widehat{p} \right\rangle\right|_{x_a+L_x} =\left.\left \langle \widehat{p} \right\rangle\right|_{x_b}-\left.\left \langle \widehat{p} \right\rangle\right|_{x_a}-\hat{\varPi}_x L_x. \end{equation}
Substituting (E6) into (E5) and considering the Reynolds decomposition 
$\left \langle \widehat {uv} \right \rangle =\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle +\left \langle \widehat {u_tv_t} \right \rangle$ leads to
\begin{align} &\frac{1}{L_x}\int_0^{h}\left(\left.\left \langle \widehat{p} \right\rangle\right|_{x_a}-\left.\left \langle \widehat{p} \right\rangle\right|_{x_b}\right) \,{\textrm{d}y}+\frac{1}{L_x Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y}\right|_0 \,{\textrm{d}x}+\hat{\varPi}_x h\nonumber\\ &\quad = \left.-\frac{1}{L_x}\int_{x_b}^{x_a+L_x} \left \langle \widehat{u} \right\rangle\left \langle \widehat{v} \right\rangle\right|_{h} \,{\textrm{d}x}-\left.\frac{1}{L_x}\int_{x_b}^{x_a+L_x} \left \langle \widehat{u_tv_t} \right\rangle\right|_{h} \,{\textrm{d}x}+\frac{1}{L_x Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_{h} \,{\textrm{d}x}. \end{align}
The first two terms of (E7) are half of the ensemble time-averaged coefficients 
$C_p(\tau )$ and 
$C_{ca}(\tau )$ as we only consider the bottom channel wall. Therefore (E7) reads
\begin{align} &\frac{C_p(\tau)+C_{ca}(\tau)}{2}+\hat{\varPi}_x h=-\frac{1}{L_x}\left.\int_{x_b}^{x_a+L_x} \left \langle \widehat{u} \right\rangle\left \langle \widehat{v} \right\rangle\right|_{h} \,{\textrm{d}x}-\frac{1}{L_x}\left.\int_{x_b}^{x_a+L_x} \left \langle \widehat{u_tv_t} \right\rangle\right|_{h} \,{\textrm{d}x}\nonumber\\ &\quad +\frac{1}{L_x Re_p}\int_{x_b}^{x_a+L_x} \left.\frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \right|_{h} \,{\textrm{d}x}. \end{align}
Figure 31. Sketch of the rectangle 
$S_z$ of boundary 
$\partial S_z$ used for (E4). The tangent and normal unit vectors to 
$\partial S_z$ are depicted as arrows with full tip. Only the normal unit vector is indicated. The locations 
$x_a$ and 
$x_b$ are the abscissae of the left and right vertical edges of the bar. The dashed line indicates the mouth of the cavity located at 
$y=h$. Only the bottom wall is shown. The distance between bars is not drawn to scale.
As in Gharib & Roshko (Reference Gharib and Roshko1987), the crest skin-friction drag is two orders of magnitude smaller than the total drag 
$-\varPi _x(\tau ) L_y=C_p(\tau )+C_{f,ca}(\tau )+C_{f,cr}(\tau )\approx C_p(\tau )+C_{f,ca}(\tau )$. Equation (E8) becomes
\begin{equation} C_p(\tau)+C_f(\tau) \approx -\frac{L_y}{L_x}\int_{l_m} \left \langle \widehat{u} \right\rangle\left \langle \widehat{v} \right\rangle \,{\textrm{d}x}-\frac{L_y}{L_x}\int_{l_m} \left \langle \widehat{u_tv_t} \right\rangle \,{\textrm{d}x} +\frac{L_y}{L_x Re_p}\int_{l_m} \frac{\partial \left \langle \widehat{u} \right\rangle}{\partial y} \,{\textrm{d}x}. \end{equation}
Equation (E9) shows that the time-ensemble average of the total drag is due to the product 
$\left \langle \widehat {u} \right \rangle \left \langle \widehat {v} \right \rangle$, the Reynolds stresses 
$\left \langle \widehat {u_tv_t} \right \rangle$ and the friction term 
$Re_p^{-1}\partial \left \langle \widehat {u} \right \rangle /\partial y$ along the mouth of the cavity.
