Min et al. found that an upstream traveling wave of blowing and suction can achieve a sustained sub-laminar drag in a turbulent channel flow [38]. However, some studies have proved that this kind of traveling wave can generate a flux opposed to the wave propagation direction [40], and it is more like a pumping than drag reduction effect [39]. Inspired by the streamwise traveling wave, we make the wave reverse in the spanwise direction periodically to mimic the oscillatory wall motion in WOS, shown in Eq. 2. v ^ w = V ^ 0 ⁡ cos k ^ z z ^ − ω ^ t ^ , 0 < t ^ < T ^ a V ^ 0 ⁡ cos k ^ z z ^ + ω ^ t ^ + ϕ ^ 0 , T ^ a < t ^ < 2 T ^ a (2) where ϕ ^ 0 = − 2 ω ^ T ^ a is a constant phase as a result of the continuity at t ^ = T ^ a In the first half period 0 < t ^ < T ^ a , the wave travels in +z direction and generates a−z direction spanwise motion, while in the second half period T ^ a < t ^ < 2 T ^ a , it travels in −z direction and generates a+z direction spanwise motion, as shown in Figure 1. With a non-dimensionalization in length scale by z = k ^ z z ^ , time scale by t = ω ^ t ^ , and velocity by v = v ^ / V ^ 0 , respectively, the actuation is rewritten as: v w = cos z − t , 0 < t < ω ^ T ^ a cos z + t + ϕ 0 , ω ^ T ^ a < t < 2 ω ^ T ^ a (3)

To resolve the WBS induced flow in a quiescent background, an asymptotic expansion method is applied as shown in Eq. 4a–c: w = w 0 + α w 1 + α 2 w 2 + O α 3 (4a) v = v 0 + α v 1 + α 2 v 2 + O α 3 (4b) p = p 0 + α p 1 + α 2 p 2 + O α 3 (4c) where α = V ^ 0 / c ^ , the ratio of the blowing and suction amplitude over the wave speed, is assumed as a small parameter, w 0 , v 0 , p 0 are the zero-order terms, w 1 , v 1 , p 1 are the first order terms, and w 2 , v 2 , p 2 and higher order terms are not considered here. Without losing generality, the zero order terms are solved in the first half period t < ω ^ T ^ a as: w 0 y , z , t = i H 2 1 − H e − y − e − H y e i z − t + c . c . (5a) v 0 y , z , t = 1 2 1 − H − H e − y + e − H y e i z − t + c . c . (5b) H = 1 − i β = a − b i (5c) a = 1 2 β 2 + 1 + 1 (5d) b = 1 2 β 2 + 1 − 1 (5e) where i = − 1 is the imaginary unit, and c.c. represents the complex conjugate. The Reynolds number β is based on the frequency ω ^ and wave number k ^ z as: β = ρ ^ ω ^ μ ^ k ^ z 2 (6)

The zero-order terms are a traveling wave in the same direction as actuation, which is also called the harmonic part in the manuscript. The harmonic part decays as the distance to the wall y increases, and the protrusion height is defined as δ ^ = 0.73 λ ^ z , where the amplitude of the velocity is only 1% of that at the wall.

The average of the first order terms is of much significance in this study: w ¯ 1 y , t = lim L → ∞ ∫ − L L w 1 y , z , t d z (7)

And it can be solved with known zero order terms (Eq. 5a, b): w ¯ 1 y , t = ∑ n = 1 ∞ f n y exp i n ω 0 t + c . c . f n y = C n ⁡ exp − i n β ω 0 y + f n , 1 y + f n , 2 y + f n , 3 y (8a) f n , 1 y = 2 1 + a − i n γ cos b y − 2 1 + a b ⁡ sin b y 2 1 + a − i n γ 2 + 4 1 + a 2 b 2 × b 2 a − a 2 + 1 β 2 a a − 1 i n π − 1 n − 1 e − 1 + a y (8b) f n , 2 y = 2 1 + a − i n γ sin b y + 2 1 + a b ⁡ cos b y 2 1 + a − i n γ 2 + 4 1 + a 2 b 2 × b 2 a + 1 β 2 a a − 1 i n π − 1 n − 1 e − 1 + a y (8c) f n , 3 y = b 2 a − 1 4 a 2 − i n γ i n π β − 1 n − 1 e − 2 a y (8d) γ = β ω 0 , ω 0 = π ω T a (8e) where the constant C n = − f n , 1 0 + f n , 2 0 + f n , 3 0 determined by the w ¯ 1 y , t = 0 boundary condition at y = 0. This asymptotic solution matches the numerical results well with a slightly smaller amplitude, as shown in Figure 1A. The expression exp − i n β ω 0 y   exp i n ω 0 t in Eqs 8a-e is the same as the Stokes layer with frequency n ω 0 .In other words, the induced oscillatory spanwise motion w ¯ 1 y , t comprises a series of the Stokes layers with a frequency n π / ω ^ T ^ a . Obviously, a large divergence exists between w ¯ 1 y , t and the Stokes layers within a very thin region above the wall as a result of different boundary conditions. The thickness of the first Stokes layer in the series (Eqs 8a-e), δ s = 4.6 2 / ω 0 β or δ ^ s ∝ T ^ a with dimensions, is also regarded as the penetration height of the spanwise motion. The higher order terms in Eqs 8a-e, or the Stokes layers with multiple frequency n ω 0 , are always neglected for analysis as their thickness is 1 / n less and the amplitude is 1 / n less. In the limit of ω ^ → ∞ , the amplitude of w ¯ 1 y , t can be approximated by W ^ m ≈ V ^ 0 2 / 8 ν ^ ω ^ . When the frequency is zero, ω ^ = 0 , there is no average spanwise motion since the traveling wave degenerates to a constant blowing and suction. It is noticed that the relative error of asymptotic solutions is approximately less than 10% when the parameter α = V ^ 0 k ^ z / ω ^ satisfies α < 0.5 . As a result, the oscillatory spanwise motion in Eqs 8a-e is called the Stokes part in the rest of the manuscript.

In summary, the motion induced by the spanwise traveling wave can be approximately decomposed into a harmonic part (Eq. 5a, b) and a Stokes part (Eqs 8a-e) with the asymptotic solution.

Incompressible fully developed turbulent channel flows with R e τ = u τ h / ν ≈ 180 ( R e b = ρ U b h / μ = 2800 ) and R e τ = u τ h / ν ≈ 550 ( R e b = ρ U b h / μ = 10000 ) are chosen as the base flow for control, where u τ = τ w /   ρ is the friction velocity and U b is the bulk velocity of the channel. For the R e τ = 180 channel flow, the computation domain is set as 0,2 π h in streamwise (x) direction, − 1,1 in wall-normal (y) direction, and 0 , π h in spanwise (z) direction, with h the half channel height. In this study, u , v , w represents the velocities in x , y , z direction, respectively. The grid points are N x × N y × N y = 113 × 169 × 113 with a corresponding resolution ∆ x + ≈ 10.1 , ∆ z + ≈ 5.0 . Here, the “+” superscript represents inner-scale variables normalized by u τ and δ v = ν / u τ . A hyperbolic tangent function y = ⁡ tanh r 2 y i − 1 / ⁡ tanh r , y i = i − 1 / N y − 1 decides the grid distribution in wall-normal direction. The resolution 0.2 < Δ y + < 5.2 is achieved with a stretch ratio r = 2.4 . And for the channel flow with R e τ = 550 , the computation domain is 0,3 π h × − 1,1 × 0,2 π h to resolve the large-scale structures [41, 42], with the resolution ∆ x + ≈ 7.6 , ∆ z + ≈ 7.7 , and 0.04 < Δ y + < 2.9 .

High order methods based on CPR(correction procedure via reconstruction) [43] are used here and seventh order polynomials are used for space discretization. For the viscous term, it is discretized with BR2 (the second approach of Bassi and Rebay) [44] and interior penalty method. The time advancement strategy is a third order explicity Euler scheme with a time step d t = 1 × 10 − 4 ( d t + = d t ∙ u τ / δ v ≈ 0.00116 ). A constant flow rate (CFR) method is chosen in the computations of the channel flows. All the computations in this study are conducted with this type of direct numerical simulations (DNS).

The length and velocity scales are normalized by h and U b in the turbulent channel flows without particular specification, respectively. Comparisons of our numerical results with those from Kim et al. [45]( R e τ = 180 , and with those from Jiménez et al. [46]( R e τ = 550 ) in Reynolds stress of the uncontrolled channel flows is done firstly for verification of our computation code, which show a satisfying match in Figure 2. The vortex structures detected by Q-value also show a good match with previous work, see Figures 3A, B.

The variables in WBS control are decomposed into three parts: a x , y , z , t = a ¯ y + a ∼ y , ϕ + a ″ x , y , z , t = a y , ϕ + a ″ x , y , z , t (10) where a ¯ is the ensemble average, a y , ϕ the phase average with ϕ representing the phase, and a ″ x , y , z , t the turbulent fluctuations a ∼ = a ¯ − a is the difference between the phase average and ensemble average, which is also called the periodic part. Fukagata, Iwamoto and Kasagi proposed an identity (known as FIK identity) relating the skin friction at the wall and the turbulent shear stress in the interior field. While in WBS control, it is rewritten as: 2 C f b − C f u = 6 R e b ︸ C f L + 3 ∫ − 1 1 1 − y − u ″ v ″ ¯ ︸ C f T d y + 3 ∫ − 1 1 1 − y − u ∼ v ∼ ¯ d y ︸ C f P (11) where C f L is the skin friction in uncontrolled flow, C f T the contribution from turbulent shear stress, and C f P the contribution from the shear stress formed by the periodic parts, respectively.

The turbulent skin friction part can be further decomposed into two parts [23, 24], one from the inner region ( y + < y P + ≈ 30 ) and the other from the outer region ( y + > 30 ), as shown in Eq. 12. 3 ∫ − 1 1 1 − y − u ″ v ″ ¯ d y = ∫ − 1 y P 1 − y − u ″ v ″ ¯ d y + ∫ y P 1 1 − y − u ″ v ″ ¯ d y (12)

Therefore the drag reduction rate can be decomposed as shown in Eqs 13a-d: 2 ∙ D R = D R I + D R O + D R P (13a) D R I = 1 C f 0 3 ∫ − 1 y P 1 − y u ″ v ″ ¯ 0 − u ″ v ″ ¯ d y (13b) D R O = 1 C f 0 3 ∫ y P 1 1 − y u ″ v ″ ¯ 0 − u ″ v ″ ¯ d y (13c) D R P = 1 C f 0 3 ∫ − 1 1 1 − y − u ″ v ″ ¯ d y (13d) where D R I , D R O , and D R P represents the drag reduction rate from the inner region, outer region, and periodic part, respectively. The factor “2” on the left-hand side of Eqs 13a-d comes from the one-side control in the channel flows. The drag reduction decomposition of the cases is shown in Table 2. The D R is mainly from D R I and D R O . D R P is always negative, which means the periodic blowing and suction has a negative effects to the drag reduction. On the other hand, compared to D R I and D R O , D R P is much smaller and even can be neglected when R e τ = 180 . While R e τ = 550 , the D R P becomes larger and even dominates in Case WBS-3. The strength of the oscillatory spanwise motion W + is larger in Case WBS-3 compared to other cases with R e τ , but D R decreases a lot to a negative value. This means the negative effects of the periodic part on D R in this case become significant. The wavenumber increases from k z = 8 (Case WBS-6) to k z = 16 (Case WBS-5), the increases in D R are mainly from D R O , as a result of a smaller protrusion height δ = 0.73 λ z . In addition, the D R gets smaller mainly from the decreases of D R O when R e τ increases, similar to the results in WOS control [23, 24].

The effects of D R P can be described by the lift-up mechanism from the harmonic parts of blowing and suction. Based on the asymptotic expansion (Eq. 4a–c), the streamwise velocity can be written as u = u 0 + α u 1 + α 2 u 2 + O α 3 , and its equation is: ∂ u ∂ x + α v ∂ u ∂ y + w ∂ w ∂ z = − ∂ p ∂ x + 1 R e ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 (14)

Since it is uniform in streamwise direction, there is no ∂ / ∂ x terms in the equation, and u can be solved with a known v , w . In addition, u 0 is the prescribed baseflow and not dependent on v , w . Here, it is assumed that u 0 = K y without losing generality, since the linear relation is always true in the linear law region [48]. Thus, the first order term u 1 can be solved with v 0 , w 0 : u 1 = K i H 2 1 − H e − y − e − H y − R e 4 H 1 − H y e − H y exp i z − t + c . c . (15)

It is obvious that u 1 < 0 above the blowing regions and u 1 > 0 above the suction region, which is thought as a lift-up mechanism. Thus, there forms a negative shear stress u v , which is closely related to the C f P . Further, the spanwise average of the shear stress, noted as S y , can be deduced: S y = lim L z → ∞ 1 L z ∫ 0 L z u v d z = S 0 y + α S 1 y + O α 2 (16a) S 0 y = 0 (16b) S 1 y = lim L z → ∞ 1 L z ∫ 0 L z u 1 v 0 d z = K 1 − a 2 + b 2 S y 1 + S y 2 (16c) S y 1 = − a 2 + b 2 sin ⁡ b y + R e cos ⁡ b y − 2 a b ⁡ sin ⁡ b y 2 a 2 + b 2 y + b ⁡ cos ⁡ b y + a ⁡ sin ⁡ b y ] e − 1 + a y (16d) S y 2 = − b + a R e 2 a 2 + b 2 y e − 2 a y (16e)

In addition, the skin friction from the periodic part, C f P , can be evaluated as: C f P ≈ 3 ∫ − 1 1 1 − y − S y d y ≈ 3 K ∫ − 1 1 1 − y − α S y K V 0 2 d y (17)

As shown in Figure 4, the S y / K is always negative, leading to a positive C f P . As a result, the periodic parts always increase the skin friction, just as shown in Table 2. C f P decreases with an increasing frequency ω , which is also compliant with the results of Case WBS-2 and WBS-3. It is also found that C f P ∝ K , which means that the larger the shear K = ∂ U / ∂ y is in the background flow, the larger C f P is. When R e τ changes from 180 to 550, K increases and leads to a larger C f P , just as the results in Table 2.

As discussed in Control Strategy section, the induced flow by the traveling wave of blowing and suction with a periodically reversed wave speed can be decomposed into two parts, one is the harmonic part (zero order terms) and the other is the Stokes part (first order terms). The harmonic part will induce a negative shear stress, as shown in Eq. 16c, and lead to a positive contribution the skin friction through C f P . An effective way to reduce C f P is increasing ω , while an increased ω will lead to a smaller strength W = V 0 2 / 8 ν ω of the Stokes part. Therefore, there should be an optimal frequency in the WBS control, which is not discussed in details in this manuscript. In addition, C f P becomes larger with an increasing mean flow shear ∂ u / ∂ y in the viscous sublayer when R e τ changes from 180 to 550. This is thought to be the first reason that the drag reduction rate is smaller in R e τ = 550 cases.

In this section, a Fourier expansion based on spanwise direction is imposed to the fluctuation field variables, u ″ = ∑ κ = 1 ∞ u ^ κ exp i κ z + c . c . Thus the stress can be decomposed into different scales as: u i ″ u j ″ ¯ = ∑ κ = 1 ∞ u ^ i κ exp i κ z + c . c . ∑ κ = 1 ∞ u ^ j κ exp i κ z + c . c . ¯ = ∑ κ = 1 ∞ u ^ i κ u ^ j κ * ¯ (18) where the superscript “*” represents the complex conjugate. As a result, the skin friction coefficient from the turbulent shear stress, C f T can be also decomposed: C f T = 3 ∫ − 1 1 1 − y − u ″ v ″ ¯ d y = ∑ κ = 1 ∞ 3 ∫ − 1 1 1 − y − u ″ v ″ ¯ κ d y ︸ C f T κ (19)

Thus, the contribution from the scale λ = 2 π / κ to the turbulent skin friction is C f T κ . The pre-multiplied spectrum of turbulent shear stress, κ u ″ v ″ ¯ κ , and the skin friction coefficients, κ C f T κ , of Case WBS-5 with R e τ = 180 and Case WBS-2 with R e τ = 550 are shown in Figure 5. The skin friction of all scale structure decreases after control in R e τ = 180 case, while that of the medium scale, 100 < λ + < 700 , decreases in R e τ = 550 case. In addition, the spectrum κ u ″ v ″ ¯ κ gets reduced nearly in the whole region, y / h < 1 , while the spectrum gets reduced in the region y + < 30 . As the R e τ increases, the contribution to the skin friction in the outer region increases, while the spanwise motion mainly affects the inner region, so the drag reduction rate is smaller in R e τ = 550 than that in R e τ = 180 , as shown in Figure 5. To further analyze the control effects on the structures with different scales, we use a cut-off wavelength to divide the whole scale into large-scale (noted as LS) and small-scale (noted as SS) structures. The cut-off wave length is chosen as λ z c + = 300 , which separates the two local maximum in the pre-multiplied spectrum, that is: u ″ = ∑ κ = 1 ∞ u ″ κ = ∑ κ < 2 π / λ z c u ″ κ ︸ L S + ∑ κ > 2 π / λ z c u ″ κ ︸ S S (20)

It was found that the vortex structures are inclined by the Stokes layers in WOS, which has a close relation to the drag reduction mechanism [14, 15, 17–19]. In our previous analysis based on the R e τ = 180 channel flow, the inclination is thought to be a representation of the stretch of the streamwise vortex, ∂ u ″ / ∂ x . ∂ u ″ / ∂ x gets reduced after control, leading to a weakened pressure-strain term ϕ 11 = p ″ ∂ u ″ / ∂ x ¯ . Thus, the energy transported into the wall-normal stress v ″ v ″ ¯ from the streamwise stress u ″ u ″ ¯ is reduced as a result of the incompressibility, ϕ 11 + ϕ 22 + ϕ 33 = 0 [15]. Here we revisited the inclination of vortex structures in the controlled flow with R e τ = 550 . Figure 6 shows the streaks at y + ≈ 11 in Case WBS-2 in one period. The inclination of the streaks is similar as that in R e τ = 180 case. Also, the inclinations of streaks can be described by the inclinations of the two-point correlation, C ∆ x , ∆ z , which is defined as: C ∆ x , ∆ z = u ″ x , y , z , t u ″ x + ∆ x , y , z + ∆ z , t ¯ u ″ x , y , z , t u ″ x , y , z , t ¯ (21)

The contours of the streaks and C ∆ x , ∆ z of SS and LS structures are shown in Figure 7. Thus, an inclination angle, noted as θ in radians, is defined as the angle between the axis of the iso-value line of C ∆ x , ∆ z and the streamwise direction to give a quantitative description of the inclination of streaks. The angle θ y , t of whole scales and separated scales (SS and LS) in Case WBS-2 are shown in Figure 8. It is shown that the angle θ y , t has a phase difference at different wall distance y , which is induced by the spanwise motion [12, 14, 15]. The angle is largest at y + ≈ 11 , recalling that the spanwise motion affects the structures at this region most. Further for SS and LS structures, the angle of SS is nearly the same as all scale (AS) and much larger than that of LS. This tells us that the spanwise motion affects the SS structure more than the LS. In addition, the inclination of SS and LS structures has a phase delay.

In physical space, the oscillatory spanwise motion or the Stokes part mainly weakens the inner region ( y + < 30 ), as shown in Figure 5. When the R e τ gets larger with an increasing dominance of the outer region, the contribution from the inner region gets reduced, leading to a smaller drag reduction rate. In spectrum space, the Stokes part mainly affects the small-scale structures ( λ + < 300 ). While the R e τ gets larger, the large-scale structures become more significant, and the drag reduction is thus reduced.