Drag Reduction of Turbulent Boundary Layers by Travelling and Non-Travelling Waves of Spanwise Wall Oscillations

: Turbulence control in the form of a streamwise travelling wave of transverse wall motion was studied numerically by employing direct numerical simulations (DNS). Both total and phase averaging were utilised to examine the statistical behaviour of the turbulence affected by the wall forcing, with a focus on the skin friction. Comparison with results from pure temporal and spatial wall forcing are conducted, and a compilation of data is used to explore analogies with drag-reduced channel ﬂow.


Introduction
The first observations of wall oscillation as means for drag reduction (DR) were made by [1] through direct numerical simulations (DNS) of a channel flow. Since then, large research efforts were devoted to internal flows, such as in a channel or pipe flow. The boundary layer flows only recently started to be investigated numerically ( [2][3][4][5][6][7][8][9]), while experiments were previously performed by a number of researchers. In fact, the first experimental evidence that confirmed Jung's DNS results was provided by [10], who applied the oscillation technique to the boundary layer flow. Most of the experimental investigations since then focused on the boundary layer, see references in [11]. Extensive comparison between DNS, using the same numerical code as in the present work, and these experiments were made by [2]. A review of both internal and external flow controlled by wall oscillations is provided in [11].
The most commonly used form of control using spanwise wall-forcing is realised with temporal wall oscillations, which are imposed through a wall velocity (W) in the spanwise direction in the form of W = W m sin(ωt), where W m is the maximum wall velocity and ω is the angular frequency of the wall oscillation, which is related to the period (T) through ω = 2π/T. The amplitude of the wall velocity does usually not depend on the streamwise coordinate (x) in computational studies of channel flow, while the forcing (1) is complemented with a function describing the spatial distribution when applied to a limited section of the wall in the boundary layer experiments and simulations. Further description is provided below. The oscillation in time is relatively straightforward to implement in an experimental setting, however, a positive energy budget may not be easily obtained. Instead, researchers [3,5,[12][13][14][15] considered a steady variation in the streamwise direction along the plate instead of a time-dependent forcing. In this case, the wall velocity (W) is imposed in the form of W = W m sin(κx), where κ is the wavenumber of the spatial oscillation, which is related to the wavelength (λ x ) through κ = 2π/λ x .
Quadrio et al. [16] first studied (through DNS on a channel flow) the combination of spatial and temporal wall oscillation (streamwise travelling waves). The travelling waves as a wall forcing can be implemented in the form: The only experimental data to date for this type of wall forcing are provided by [17], who applied the streamwise travelling wave on the pipe flow, and by [18] for the high-Reynolds number boundary layer.
The appropriate scaling of the oscillation parameters are provided with the wallscaling, using the friction velocity u τ and kinematic viscosity ν. Hence, the parameters are discussed in the form of: The diminishing turbulence as a consequence of a moving wall underneath the flow may seem counter-intuitive, since the disturbance via the moving wall could be seen as a source for creating chaotic flow near the wall. However, the motion of the wall disturbs the intricate balance between low-speed streaks and streamwise vortices in the near-wall flow, thereby disrupting the self-generating mechanism of turbulence, which leads to diminished turbulent activities and decreased skin-friction. However, the precise explanation for the DR and how to best utilise the effects of a moving wall are still under debate. Please refer to the recent review paper [11] for an in-depth discussion of the various proposed physical mechanisms offering explanations for the reduced drag.
In the current work, previously produced DNS data of boundary layer flow controlled by all three forms of forcing above are presented. In addition, a complementary simulation (oscillating wall according to Equation (1) with T + = 100) was performed to clarify some questions regarding the dependence of DR on the oscillation parameters. The objective of the present work is to utilise the data available for the three variants of the wall-forcing at identical Reynolds number to clarify the differences in the spatial development of the DR in the downstream direction, starting from the onset of the wall-oscillations. The study is focused on the boundary layer forced by the travelling wave, since the results from the only numerical simulation that exists were not presented in full to date. Part of the results described herein were presented at the TSFP9 conference (Ref. [19]). Another aim of this paper is to collect results regarding the maximum DR obtained in various simulations that were never presented in a concise form previously. In this way, the present paper constitutes a complement to the recent review article [11].

Methodology
The numerical code and grid are the same as in the previous simulations of spatially and temporally oscillating turbulent boundary layers, see the references in Table 1. The code was developed at KTH, Stockholm [20] and was verified by comparison with experimental results for both transitional and turbulent flows.

Numerical Scheme
A pseudospectral method is employed, with Fourier discretization used in the streamwise and spanwise directions and Chebyshev polynomials in the wall-normal direction. The simulations start with a laminar boundary layer at the inflow that is triggered to transition by a random volume force near the wall. A fringe region is added at the end of the computational domain to enable simulations of spatially developing flows. In this region, the flow is forced from the outflow of the physical domain to the inflow. In this way, the physical domain and the fringe region together satisfy periodic boundary conditions. The time integration is performed using a third-order Runge-Kutta scheme for the nonlinear terms and a second-order Crank-Nicolson method for the linear terms. A 3/2 rule is applied to remove aliasing errors from the evaluation of the nonlinear terms when calculating fast Fourier transforms (FFTs) in the wall's parallel plane.

Numerical Parameters
All quantities are nondimensionalised by the freestream velocity (U) and the displacement thickness (δ * ) at the starting position of the simulation (x = 0), where the flow is laminar. The Reynolds number is set by specifying Re δ * = Uδ * /ν at x = 0. In all the simulations presented here, Re δ * = 450. The computational box is 600 in simulation length units (δ * ) long (including 100 units for the fringe), 30 units high, and 34 units wide. The resolution used for the simulations were 1000 modes in streamwise direction, 217 modes in wall-normal direction, and 200 modes in the spanwise direction. This grid size resulted in a spatial resolution of ∆X + × ∆Z + × ∆Y + mean = 13.6 × 3.9 × 3.1, with ∆Y + min close to the wall at 0.036. The + superscript indicates that the quantity is made nondimensional with the friction velocity of the reference case at x = 250, denoted u 0 τ , and the kinematic viscosity (ν). The number of modes is slightly lower for the simulations taken from [4] (800 streamwise; 201 wall-normal; 144 spanwise). However, this coarser resolution was shown to be sufficient to obtain even the higher-order turbulence statistics correctly represented [2].
The sampling time for the reference case was 10,000 in time units (δ * /U), started only after a stationary flow (in the statistical sense) was reached. In the case with wall forcing, the total sampling time was 23,000 after an initial simulation time of 4000 with oscillations.
As the fringe starts at x = 500, only results up to x = 470 are used to avoid any upstream influence of the fringe. This position is located upstream of the wall-forcing end-point; hence, the recovery region from the drag-reduced state is not included in these simulations. The transition region is roughly between x = 5 (where the trip is located) and x = 150. Thus, the region of a fully developed turbulent boundary layer, free from any influence of the numerical method, is x = 150 − 470. The Reynolds number based on the momentum thickness (Re θ ) varies between 390 and 750 in this region for the reference boundary layer. In inner-scaling (based on the friction velocity at x = 250), the region amounts to about 7000 wall units.
The wall forcing in the present simulation (oscillating wall with a period of T + = 100) and earlier simulations (temporal, spatial, and travelling wave) is applied in the spanwise direction over a particular region in streamwise direction. Therefore, the form of this boundary condition in the case of a travelling wave is given by where a profile function f (x) is utilised to select the domain where the oscillation takes place, and it is given by: with x start , x end , x rise and x f all set to 250, 487, 5, and 5, respectively. S(x) is a continuous step function that rises from zero for negative x to unity for x ≥ 1. The expression of S(x), which has the advantage of having continuous derivatives of all orders, is: In the case of pure temporal and spatial oscillations, the first and second part of Equation (7) (inside the square brackets), respectively, vanishes. The remaining parameters for all the cases presented in Section 3.1 are given in Table 1. The Reynolds number at the onset of oscillations (at x = 250) is Re θ = 505. The simulations are performed with constant W m , meaning that W + m is increasing downstream since the friction velocity u τ decreases downstream. An exception to this will be presented in Section 3.2. To distinguish the time averaged statistical mean from the phase averaged mean, a total of 36 individual statistics were created, i.e., a total of 36 bins were used to resolve one period of the oscillations.

Results
In the following, the streamwise, wall-normal, and spanwise coordinates are denoted x, y and z, respectively, while the corresponding mean velocity components are denoted as u, v, and w, with the mean taken over time and in the spanwise direction.

Compilation of Old and New Data
The friction coefficient is defined as: where u τ is the friction velocity: and ν is the kinematic viscosity. The DR is calculated from: where C 0 f = 2(u 0 τ /U) 2 is the skin friction of the reference boundary layer. The DR according to the above expression will vary downstream in the boundary layer. The C 0 f and C f are evaluated at the same downstream location (x). Other definitions could be based on, e.g., the skin friction coefficients at identical Re θ or boundary layer thickness.
In Figure 1a the DR is shown as a function of the streamwise coordinate (x). The DR reached with the travelling wave as forcing is around 42%, slightly less than what was obtained (45%) in channel flow, with identical T + , λ + x and W + m , by [16] at a similar Reynolds number. The Reynolds number for the channel flow simulations by [16] was Re τ = 200, which translates to Re θ = 464, as compared to Re θ = 500 − 730 in the boundary layer simulations presented here. To convert the Reynolds numbers, the expression proposed by [21], is used. The parameters T + = 176 and λ + x = 384 for the TW case were chosen for the boundary layer control since they are in the region of maximum drag reduction for the channel flow [16]. The region with maximum or near-maximum DR is larger and less clearly defined in the two-parameter space for the travelling wave forcing than the corresponding span of only frequency (or period) and wavenumber (or wavelength) for the pure temporal or spatial forcing, respectively.
In addition, previous and new simulations with pure temporal forcing and pure spatial forcing are shown for comparison Figure 1a. The previous boundary layer results are taken from [4,5,9] (see Table 1). All of these cases were performed with the same amplitude (W + m = 12) as in the travelling wave simulation, except the case with T + = 67, which has an amplitude of W + m = 11.3. We return to this special case in Section 3.2. In Figure 1a, the pure temporal forcing at T + = 176 does not produce a strong DR, while this value of the period is near optimum when employing the travelling wave. The optimal period is lower for the temporal forcing, around T + = 100 (the black line in Figure 1a), which agrees with the value obtained from channel flow at similar Reynolds numbers. However, the maximum DR obtained with pure temporal oscillations is lower than the values obtained from the travelling wave (the red line in Figure 1a). On the other hand, a pure spatial forcing (the blue dotted line in Figure 1a), with a wave length close to the optimum, creates a DR close to the travelling wave. Again, the value of the wave length for producing maximum DR for a pure spatial forcing is far from the optimum when considering the travelling wave. Note that the pure spatial forcing (λ + x = 1320) produces a spatially fluctuating DR. On the other hand, the travelling wave forcing yields a nonoscillating profile of the DR since the spatial variation is averaged out when performing the temporal averaging.
Even lower DR than the values obtained by T + = 176 is obtained by T + = 30 (the green line in Figure 1a). The periods closer to the optimum T + = 100, namely T + = 67 and T + = 132 are close to each other and well below the DR for T + = 100, which reflects that the optimum is likely to be around T + = 100.
The same DR results are plotted with the Reynolds number (Re θ ) as a downstream coordinate is shown in Figure 1b. The figure shows that the drag-reduced boundary layers are growing less downstream than that of the uncontrolled flow; hence, Re θ is lower than for the more strongly drag-reduced boundary layers at the same downstream position (identical x).
In Figure 2, the maximum DR for the pure temporal cases are compared with that of the corresponding channel flow taken from [22]. The DR is consistently slightly lower for the boundary layer as compared with that of the channel flow, which agrees with the earlier mentioned travelling wave results. In addition, the peak (optimal) DR seems to be more pronounced in the boundary layer geometry than for the channel flow, for which the data indicate a plateau where near-optimal DR values are obtained for a wide range of oscillation periods.

Downstream Variation
The special amplitude W + m = 11.3 of the forcing was chosen by [4] to compare with that of the experimental data by [23]. However, the Reynolds number was lower in [4] (Re θ = 505), yielding a higher DR than in the experiments performed at Re θ = 1400. On the other hand, Ref. [9] repeated the DNS of [4] at Re θ = 1400 and obtained identical DR as in the experiments [23]. In addition, Ref. [9] repeated the same simulation but at a third Reynolds number (Re θ = 375) to complement the earlier simulations at Re θ = 505 and Re θ = 1400. The trend of the maximum drag reduction DR max with increasing Reynolds number followed DR max (%) ∝ Re −0.153 θ , with the exponent obtained from the three boundary layer simulations considered. The slope is indicated with the black line in Figure 3, where the maximum DR for the three cases are indicated with the filled circles. In the figure, the development prior to reaching maximum DR was excluded for clarity (also the recovery stage is excluded as previously mentioned). However, the local drag reduction (considering the three boundary layers separately) declined more severely downstream in all three boundary layers. In Figure 3 the three cases are the green, blue, and red solid lines. The blue curve in Figure 3 is the same datum as the cyan line in for the boundary layers at Re θ = 375, Re θ = 505 and Re θ = 1400, respectively. Hence, the conclusion after comparing all three simulations was that the decline of DR downstream of the point of maximum DR is greater than what can be explained by the increase in Reynolds number over the same stretch in the downstream direction (which is the DR max (%) ∝ Re −0.153 θ behaviour discussed above). In addition, a simulation was conducted with the same parameters, but where the amplitude of the oscillation varies downstream such that W + m = 11.3 is constantly downstream of the onset of the forcing. This case is shown as the green broken line in Figure 3. All four cases described are summarised in Table 2, since these data were not presented together before. Included also are the measurements [23] of the same case and the channel flow simulation [22] with the same oscillation parameters. That the DR is slightly lower for the case with constant W + m is natural since the amplitude itself (W m ) must decrease downstream to keep W + m constant. The maximum DR (indicated with the open green circle in Figure 3) occurs slightly upstream of the case with constant W m . The relatively large discrepancy between the maximum DR from numerical [9] and experimental [23] data in Table 2 is slightly misleading, since other data points in the experimental investigation show more agreement with the numerical data (see Figure 2 in [9], where also the DNS data from the four boundary layers are shown in linear scale). Tables 1 and 2 form together with  Table 6 in [11] a complete set of the maximum drag reduction from zero-pressure-gradient boundary layer simulations. Note that OW4 in Table 1 is identical to row number three in Table 2. There are small differences in the numbers recorded here and the original publications, depending on how the maximum is defined for the spatial case (taking the maximum as in the peak value or averaging the troughs and valleys before calculating the maximum), or, in some cases, updated reference values from a simulation of the uncontrolled boundary layer with higher resolution were used. Nevertheless, these slight differences in numbers carry little significance.
The exponent −0.153 obtained from the boundary layer simulations (as described above) is approximately the same as the exponent estimated from channel flow simulations in Figure 4b in [24] and using expression (13). Hence, the decline of the maximum dragreduction in boundary layer flow (indicated by the black line in Figure 3) is similar to that derived from channel-flow simulations at a similar oscillation period. This finding highlights the equivalence between the two geometries when comparing the maximum drag-reduction in the boundary layer with the single drag-reduction margin obtained from channel flow, at least for the Reynolds numbers considered here.

Phase-Wise Drag Reduction
When considering only one of the 36 separate sets of statistics (i.e., considering approximately a single phase in time of the forcing) from the travelling wave case, the wall velocity remains constant in time but varies sinusoidally in space. An example of the wall velocity is given in Figure 4a. Included is also the corresponding skin friction profiles. The C f is fluctuating with half the wavelength, and the crests and troughs correspond roughly to the maximum/minimum and zero-crossings in all velocities, respectively. Thus, the same correlation as found in the spatial forcing by [3] is generated phase-wise by the travelling wave. The less-converged C f profiles are such because only 1/36 of the total statistics is being used. The spatial oscillations seen in the phase-wise statistics cancel each other in the total statistics, which is why the DR profile in Figure 1 (red curve) is nonoscillating.
To investigate more quantitatively the relationship between the wall velocity and C f , all the phases are added together, taking into account the phase shift. The result is shown in Figure 4b. Here a small phase shift is revealed, similar to the results for the temporal forcing in channel flow presented by [25]. By examining the phase shift in Figure 4b carefully, the distance can be quantified to be ∆x = 1.2. Hence, the spatial phase shift is ∆x = λ x /14. From the plots of the time variation of the skin friction and wall velocity in [25], one can calculate that the temporal phase shift in their channel flow is ∆T = T/14, i.e., exactly the same shift in phase between the skin friction and wall velocity occurs.  An interesting observation from the temporal wall oscillation case with the long period of T + = 176 is the wave occurring in the skin friction profile, from the onset of the wall-forcing, as illustrated in Figure 5, where C f is plotted for two phases (φ = 0 • and φ = 90 • ). Hence, there is a propagation wave originating from the onset of the oscillating wall, travelling downstream with a decreasing amplitude, until vanishing roughly where the skin friction stabilises to its drag-reduced levels. The period of the wave is half of the forcing period, which can be seen in the video accompanying this paper. The same period (half of the forcing period) can be seen in Figure 13 in [6]. When studying the corresponding skin friction profiles from the case with T + = 30, it is observed that the wave length is shorter (Figure 6), while the period remains at half of the forcing period (not shown). The inset in the figure is an enlargement of the plot within the transient region of the skin friction reduction. By using data also from the forcing with T + = 100, the wavelength (λ + ) of the propagation wave can be plotted as a function of the forcing period, see Figure 7. However, the values of the wavelength are estimated from plots with nonconverged statistics (see Figures 5 and 6). Nevertheless, the evidence seems to point towards a linear relationship. The wave propagation cannot, however, be seen in the phase-wise statistics from the travelling wave forcing, see Figure 8. This is obviously a consequence of the already wavy variation of the wall forcing, and consequently, the skin friction for each phase. The inset on the bottom-left shows that the maximum and minimum alternate at the same x-position, but the wave propagation is difficult to discern for this forcing case. On the other hand, further downstream, where the skin friction attained its drag-reduced level, a clear wavy profile for each phase is observed, which are the same curves as those illustrated in Figure 4a above. Those waves at the drag-reduced region do not exist for the temporal forcing cases in Figures 5 and 6.  Since each phase constitutes a forcing with a standing wave, we compare the amplitude of the wave from the travelling wave case with the pure spatial case (SW) in Figure 9. The wavelength of the C f differs because the forcing wavelength is different, with λ + x = 1320 for SW and λ + x = 384 for TW. However, the amplitude between valleys and crests is also different, although the forcing amplitude is identical (W + m = 12). The relative amplitude (maximum difference between the values divided by the mean skin friction) is 0.024 and 0.030 for the spatial case (SW) and travelling wave case (TW), respectively.  Figure 9. Skin friction coefficient for the travelling wave forcing (TW) (one phase is shown in red) together with skin friction coefficient for the pure spatial forcing (SW) (blue).

Reynolds Stresses
So far, only the streamwise statistics in the form of skin friction were presented. Among the wall-normal statistical data of interest are the velocity profiles and the Reynolds stresses. The former were shown to be altered by the DR in a well-defined manner, and a relationship between the slope of the logarithmic part of the profile and the amount of DR could be derived [3]. Hence, the velocity profiles are uniquely determined by the DR and not by the form of the wall control causing the DR. On the other hand, for the Reynolds stresses no such relationship exists, and no quantitative description of the alteration of the profiles by the DR has been derived. Therefore, the root-mean-square (rms) profiles of the longitudinal (u + rms -solid curves) and wall-normal velocity fluctuations (v + rms -dashed curves), together with their covariance (Reynolds shear stress u v + -dotted curves), are shown for the cases TW (red), SW (blue) and OW1 (green) in Figure 10. The + scaling again is based on the friction velocity from the reference case; hence, the scaling reveals the absolute values of the quantities. In Figure 10, the black lines are the reference case. The green lines denote OW1 (T + = 176), which is the case with low DR, and the Reynolds stresses should hence be closer to the reference case (black lines) than that of the other cases in the figure, which is true upon inspection (except perhaps for the longitudinal stress that are discussed further down). The blue curves represent SW, which Reynolds stresses should be located far from the reference case since the DR for this case is among the best performing in terms of DR. Again, the figure reveals consistent Reynolds stresses that are significantly lower than the reference case. For TW, which is the case with the strongest DR and visualised with the red curves, the wall-normal and shear Reynolds stresses are indeed far from the reference case, and the maximum values coincide with those from SW (the case with almost identical DR values). However, as remarked above, u + rms is not consistent with the other components, since the maximum is closer to the reference case than in both SW and OW1. This may be due to the fact that the maximum of u + rms is located below the edge of the Stokes layer, while the wall-normal and shear stresses are located above it [4].  Figure 10. Reynolds stresses for travelling wave case (TW) (red), spatial case (SW) (blue), and temporal case (OW1) (green). Black curves are from the reference case. All profiles are taken from x = 450. -u + rms ; --v + rms ; · · · u v + . Friction velocity from the reference case at x = 250 (where oscillations start) was used for + scaling.

Conclusions
The compilation of data presented herein provided evidence for the close resemblance of the drag reduction (DR) in boundary layer and channel flow, although the effects of the variation downstream are lacking in the latter case. The compilation of data from several sources in the literature was combined with newly produced direct numerical simulations (DNS) data, and concrete conclusions are (the third item was already mentioned in [9]):

•
The growth of the drag-reduced boundary layer is weaker the stronger the DR is; • The optimal oscillation period for the pure temporal forcing is more well-defined than that of the channel flow; • The decline of DR downstream is more severe than what can be explained by the increase in Reynolds number; • By keeping W + m constant downstream, both the maximum DR and the values downstream are reduced; • The wavelength of the propagation wave of C f is linearly dependent on the oscillation period; • The amplitude of the C f variation is different for the purely spatial case compared to that of the travelling wave forcing, although the amplitude of the forcing is identical; • The longitudinal velocity fluctuations for the travelling wave forcing are not reduced as much as for the spatial forcing case with a similar drag reduction margin.
Although it may initially seem that the extra complication of having parameters varying downstream is of less importance, the practical applications of a drag reducing strategy are mostly in external flows, where the development of DR due to changing conditions downstream plays a role for the overall energy saving attained by the wall forcing. Hence, before spending effort on an application, the understanding of how the drag reducing effect is influenced by the growth of the boundary layer downstream is important. The present paper is a step in that direction, while further studies are obviously required, with higher Reynolds numbers, more variation of oscillation parameters, and the development of practical techniques that enable flow control without wall motion, such as plasma actuators.