1 Introduction

Viscous drag in turbulent wall-bounded flows is one of the major contributors to the overall drag of flow over slender bodies in general and passenger planes in cruise flight in particular. In its simplest form, the incompressible flow around a two-dimensional infinite wing section consists only of friction and pressure drag, both due to the viscosity of the fluid. Also in realistic flow conditions, the friction drag of the—most often turbulent—boundary layer can take a considerable share of the drag budget for civil transportation aircraft, especially when also the fuselage is taken into account (Schrauf 2005). Lowering turbulent friction drag is therefore essential to meet future \(\text {CO}_2\) reduction goals. Besides preventing fully turbulent flow and benefiting from the considerably lower laminar drag (Spalart and McLean 2011), there is substantial past and ongoing research in the field of turbulent drag reduction. Unlike active techniques, which require energy introduction to the system, passive techniques such as riblets (Walsh and Weinstein 1978; Bechert et al. 1985; Walsh et al. 1989; García-Mayoral and Jiménez 2011b, a) and compliant surfaces (Choi et al. 1997; Kim and Choi 2014; Luhar et al. 2015; Zhang et al. 2017), which are optimized for single operating conditions, yield reduced skin friction without any added energy. From an energetic point of view, this is the perfect solution. Unfortunately, the design conditions are not always satisfied. Active techniques, however, can be adjusted to a range of flow conditions such that at least certain approaches can achieve higher net power saving. These results, however, hold mostly in canonical flow setups like turbulent channel flows under laboratory conditions, i.e., at extremely low technology readiness levels. Additionally, the investigation and identification of the mechanisms of active drag reduction techniques might also benefit from the development of improved or new passive approaches, e.g., passive wavy wall undulations (Ghebali et al. 2017) inspired by the active technique of streamwise oscillations of spanwise wall velocity (Viotti et al. 2009).

In the following, we briefly review active flow control techniques that are most relevant to the current study. Particularly, we focus on methods that use either in-plane wall motion, such as forcing parallel to the wall or out of plane wall motion. In this study, we employ actuation that belongs to the latter category. The objective is to discuss its similarities and peculiarities over existing techniques, and present its drag reduction and net power saving potentials, which reach 26 and 10% compared to the unactuated flow.

Inspired by turbulence suppression by temporary pressure gradient variations (Moin et al. 1990) Jung et al. (1992) performed the first simulations of spanwise wall oscillations which resulted in significantly lowered friction drag. The method was subsequently investigated in detail in the following years for Poiseuille flow (Choi et al. 1998; Quadrio and Ricco 2004) and turbulent boundary layer flow (Ricco and Wu 2004; Yudhistira and Skote 2011; Lardeau and Leschziner 2013). Detailed analyses indicated an interaction of the oscillating spanwise shear with the near-wall velocity streaks (Touber and Leschziner 2012; Agostini et al. 2014). Furthermore, it was found that the maximum drag reduction in turbulent boundary layer flow is moderately lower than in turbulent channel flow (Lardeau and Leschziner 2013). Motivated by this simple but effective approach, other forms of spatio-temporal forcing have been developed, which is excellently discussed by Quadrio (2011). A modified variant of the purely temporal oscillations of spanwise velocity (Touber and Leschziner 2012) are spanwise traveling waves of spanwise forcing (Du and Karniadakis 2000; Du et al. 2002) or spanwise traveling waves of a flexible surface (Zhao et al. 2004). Although the techniques are different in their actuation principle, the effect of introducing oscillating spanwise shear close to the wall is alike.

Another actuation variant is spanwise traveling transversal surface waves (Itoh et al. 2006). Instead of directly introducing spanwise velocity, the surface is wavily deflected in the wall-normal direction to generate a secondary flow field of periodic wall-normal and spanwise fluctuations. Positive drag reduction using this technique was achieved experimentally (Itoh et al. 2006; Tamano and Itoh 2012; Li et al. 2018) and numerically for channel flow (Tomiyama and Fukagata 2013), boundary layer flow (Klumpp et al. 2011; Koh et al. 2015a, b; Ishar et al. 2019), and airfoil flow (Albers et al. 2019). In the experiments for turbulent boundary layer flow, the wavelengths were in the range \(2570 \le \lambda ^+ \le 4440\), the actuation periods \(T^+ \ge 110\) (Tamano and Itoh 2012; Li et al. 2018), and the amplitudes \(A^+ \le 30\) (Tamano and Itoh 2012). The \(+\) sign denotes inner coordinates, i.e., the wave parameters are scaled by the friction velocity \(u_\tau \) and the kinematic viscosity \(\nu \). Generally, an increased drag reduction for decreasing period and increasing amplitude was found. Due to the increasing computational cost for larger domains the numerical investigations generally use smaller wavelengths \(35 \le \lambda ^+ \le 1000\) (Tomiyama and Fukagata 2013; Ishar et al. 2019) but a greater period range \(10 \le T^+ \le 227.2\) (Klumpp et al. 2011; Tomiyama and Fukagata 2013). The same drag reduction trend as in the experiments was observed. Tomiyama and Fukagata (2013) observed a possible shielding effect of quasi-streamwise vortices from the wall by the wave-like deformations and showed that a combination of the thickness of the Stokes layer, i.e., the actuation period, and the actuation velocity amplitudes scales reasonably well with drag reduction.

However, the question remains what happens at higher amplitudes and wavelengths and lower periods, especially considering the vast gap between the mostly relatively short wavelength setups in numerical simulations and the large wavelengths in all experimental setups limited by mechanical actuator constraints. We will investigate if the trend of higher drag reduction for longer wavelengths (Du et al. 2002) can be confirmed. Furthermore, an optimum forcing period \(T^+\) in inner scaling was not determined for this technique and it remains an open question if one exists and if so if it is in the range of other techniques, e.g., \(T^+ \approx 70\) for spanwise oscillating wall in turbulent boundary layer flow (Lardeau and Leschziner 2013). In this study, we address these questions. We investigate the higher reduction trends with longer wavelengths and examine the flow sensitivities over the space spanned by the three actuation parameters, i.e., wavelength, wave period, and wave amplitude, using high-resolution large-eddy simulations (LES) of turbulent boundary layer flow. In total, 80 configurations are computed. The objective is to achieve drag reduction and net energy saving in the range of other actuation techniques and to compare the flow response to that from pure spanwise oscillations.

The paper has the following structure. First, the numerical method is concisely described in Sect. 2. Then the flow setup and all flow and actuation parameters are specified in Sect. 3. The results are discussed in Sect. 4. Finally, the essential results are summarized in Sect. 5.

2 Numerical Method

The actuated turbulent boundary layer flow is computed by solving the unsteady compressible Navier-Stokes equations by a large-eddy simulation (LES) formulation. To capture the temporal variation of the geometry, the equations are written in the Arbitrary Lagrangian–Eulerian (ALE) formulation (Hirt et al. 1974) such that the actuated wall can be represented by an appropriate mesh deformation. Additional volume fluxes are determined to satisfy the Geometry Conservation Law (GCL).

The discrete solution is based on a finite-volume approximation on a structured body-fitted mesh. A second-order accurate formulation of the inviscid fluxes using the advection upstream splitting method (AUSM) is applied. The cell-surface values of the flow quantities are reconstructed from the surrounding cell-center values using a Monotone Upstream Scheme for Conservation Laws (MUSCL) type strategy. The viscous fluxes are discretized by a modified cell-vertex scheme at second-order accuracy. The time integration is performed by a second-order accurate five-stage Runge-Kutta scheme, rendering the overall discretization second-order accurate.

The subgrid scales in the LES are implicitly modeled following the monotonically integrated large-eddy simulation approach (Boris et al. 1992), i.e., the numerical dissipation of the AUSM scheme models for the viscous dissipation of the high wavenumber turbulence spectrum (Meinke et al. 2002). Thus, the small-scale structures are not explicitly resolved in the whole flow domain and the grid is used as a spatial filter resolving the large energy-containing structures in the inertial subrange.

The numerical method has thoroughly been validated by computing a wide variety of internal and external flow problems (Rütten et al. 2005; Alkishriwi et al. 2006; Renze et al. 2008; Statnikov et al. 2017). Analyses of drag reduction have been performed for riblet structured surfaces (Klumpp et al. 2010a) and for traveling transversal surface waves in canonical turbulent boundary layer flow (Klumpp et al. 2011; Koh et al. 2015a, b; Meysonnat et al. 2016) and in turbulent airfoil flow (Albers et al. 2019). The quality of the results confirms the validity of the approach for the current flow problem.

3 Computational Setup

The zero-pressure gradient (ZPG) turbulent boundary layer flow over a wall actuated by a sinusoidal wave motion is defined in a Cartesian domain with the x-axis in the main flow direction, the y-axis in the wall-normal direction, and the z-axis in the spanwise direction. The velocity vector in the Cartesian frame of reference \(\mathbf {x} = (x,y,z)\) is denoted by \(\mathbf {u} = (u,v,w)\), the pressure is given by p, and the density by \(\rho \). The flow variables are non-dimensionalized using the flow quantities at rest, the speed of sound \(a_0\), and the momentum thickness of the boundary layer at \(x_0 = 0\) such that \(\theta (x_0=0) = 1\). The momentum thickness based Reynolds number is \(Re_\theta = u_{\infty } \theta / \nu = 1000\) at \(x_0\), where \(u_\infty \) is the freestream velocity and \(\nu \) is the kinematic viscosity. This corresponds to a friction velocity based Reynolds number \(Re_\tau = 360\) and a momentum thickness in inner units \(\theta ^+ = 44\). The Mach number is \(M = 0.1\), i.e., the flow is nearly incompressible. Note that unlike standard ZPG turbulent boundary layer flow, the actuated flow is statistically three-dimensional due to the wave propagating in the z-direction.

An overview of the setup is given in Fig. 1.

Fig. 1
figure 1

Overview of the physical domain of the actuated turbulent boundary layer flow. The quantities \(L_x, L_y,\) and \(L_z\) are the dimensions of the domain in the Cartesian directions, \(\lambda \) is the wavelength of the spanwise traveling wave, and \(x_0\) marks the onset of the actuation. The surface area \(A_\mathrm {surf}\) for the integration of the wall-shear stress \(\tau _w\) is shaded red

The dimensions of the physical domain are \(L_x = 190\,\theta \), \(L_y = 105\,\theta \) in the streamwise and wall-normal direction. For the spanwise direction, five domain widths, \(L_z \in [21.65\,\theta ,\)\(25.98\,\theta ,\)\(34.64\,\theta ,\)\(38.97\,\theta , 64.95\,\theta ]\) are used. The mesh resolution is \(\varDelta x^+ = 12.0\) in the streamwise direction, \(\left. \varDelta y^+\right| _{\text {wall}} = 1.0\) in the wall-normal direction with gradual coarsening off the wall up to \(\varDelta y^+ = 16.0\) at the boundary layer edge, and \(\varDelta z^+ = 4.0\) in the spanwise direction. This yields a DNS-like resolution near the wall. Away from the wall, the resolution requirements are lower such that overall an LES resolution is achieved.

At the inflow of the domain, the reformulated synthetic turbulence generation (RSTG) method by Roidl et al. (2013) is used to prescribe a fully turbulent inflow distribution with an adaptation length of less than five boundary-layer thicknesses \(\delta _{99}\). A fully turbulent boundary layer is achieved at \(x_0\), which marks the onset of the actuation. Characteristic outflow conditions are applied at the downstream and upper boundaries, whereas periodic conditions are used in the spanwise direction. On the wall, no-slip conditions are imposed and the wall motion is described by

$$\begin{aligned} y^+|_\mathrm {wall}(z^+,t^+) = g(x) A^+\cos \left( \frac{2\pi }{\lambda ^+} z^+ - \frac{2\pi }{T^+}t^+\right) , \end{aligned}$$
(1)

where \(A^+ = A u_\tau / \nu \) is the amplitude, \(\lambda ^+ = \lambda u_\tau / \nu \) is the wavelength, and \(T^+ = T u_\tau ^2 / \nu \) is the period. Variables in inner scaling based on the friction velocity of the non-actuated reference case \(N_1\), i.e., \(u_{\tau ,\mathrm {ref}}\), are denoted by the superscript \((\cdot )^+\). If the actual friction velocity of an actuated case \(u_{\tau ,\mathrm {act}}\) is used, the variables are denoted by the superscript \((\cdot )^*\).

In total, 80 variations of \(A^+ \in [0,78]\), \(T^+ \in [20,120]\), and \(\lambda ^+ \in [200,3000]\) are simulated. A detailed list of all parameter combinations can be found in Table 2 in the “Appendix”. Note that the narrowest domain has a spanwise extent of \(L_z^+ = 1000\) such that for all wavelengths \(\lambda ^+ < 1000\) multiple wavelengths are considered. A sketch of all wavelengths and the respective maximum amplitude at each wavelength is illustrated in Fig. 2.

Fig. 2
figure 2

Overview of the sinusoidal wall function with all used wavelengths and the corresponding maximum amplitude. The spanwise extent is varied to fit an integer number of wavelengths into the domain

To enable a smooth spatial transition from the stationary flat plate to the deflected wall and vice versa, the piecewise defined function

$$\begin{aligned} g(x)= & {} \left\{ \begin{array}{ll} 0 &{} \text {if} \quad x< -5 \\ \frac{1}{2}\left[ 1-\cos \left( \frac{\pi (x+5)}{10}\right) \right] &{} \text {if} \quad -5 \le x< 5\\ 1 &{} \text {if} \quad 5 \le x< 130\\ \frac{1}{2}\left[ 1+\cos \left( \frac{\pi (x-130)}{10}\right) \right] &{} \text {if} \quad 130 \le x < 140\\ 0 &{} \text {otherwise} \end{array} \right. \end{aligned}$$
(2)

is used in Eq. 1.

The drag is integrated over the wall surface within the streamwise interval \(x \in [50.0, 100.0]\) and over the entire spanwise extent. This area is colored in Fig. 1. Hence, the drag is only computed in the region where the flow is fully influenced by the traveling wave actuation. The actuated boundary layer is not impacted by the flow upstream and downstream of the actuated surface.

The computing strategy is such that first, the non-actuated reference case is simulated for \(t u_\infty / \theta \approx 650\) convective times until a quasi-steady state is observed in the integrated drag. All actuated cases are then initialized by the flow field of the reference case and the transition from a flat plate to an actuated wall flow is performed via a temporal decay controlled by \(1-\cos (t)\). Once a new quasi-steady state is observed all simulations are averaged over \(t u_\infty / \theta \approx 1250\) times.

4 Results

In the following, the results of the parameter study will be investigated in detail. First, a grid convergence study is performed in Sect. 4.1. Then, the wall-shear stress reductions as a function of the wave parameters are thoroughly discussed in Sect. 4.2. The findings are compared with data from the literature for the same and similar drag reduction techniques. This analysis is followed in Sect. 4.3 by a discussion of the variation of the total drag, i.e., the wall-shear stress multiplied by the wetted surface, since the wetted surface changes for different parameter setups. Support vector regression is used to predict drag reduction and to determine the drag reduction sensitivities for varying actuation settings in Sect. 4.4. The statistics of the non-actuated and actuated flow field are compared, and links to drag reduction mechanisms are drawn in Sect. 4.5. The spanwise shear distribution in the near-wall region with special focus on the periodic Stokes shear and its relevance for drag reduction is analyzed in Sect. 4.6. The spanwise force and flow fields are discussed in Sect. 4.7 and the net energy balance results are presented in Sect. 4.8.

4.1 Grid Convergence

To ensure a sufficient grid resolution for the large-eddy simulation of the actuated turbulent boundary layer flow a grid convergence study is conducted. The data of the three meshes that are compared are listed in Table 1.

Table 1 Summary of the grid spacings \(\varDelta x^+\), \(\varDelta y^+_\mathrm {wall}\), and \(\varDelta z^+\), the number of cells inside the boundary layer \(N_{BL}\), the number of cells in the coordinate directions, the total number of cells, and the computed drag reduction \(\varDelta c_d\) for the coarse, standard, and fine grid

Besides the mesh data, the drag reduction \(\varDelta c_d\), which is defined in Sect. 4.3, is given. The actuation parameters in inner coordinates are \(\lambda ^+ = 1000\), \(T^+ = 40\), and \(A^+ = 40\). They define case \(N_{24}\) in Table 2 in the “Appendix” for which a moderate drag reduction is achieved. It is evident from the results in Table 1 that the drag reduction values on the standard and the fine grid are nearly identical. On the coarse grid, however, a clear deviation is determined.

A comparison of the symmetric stresses and the shear-stress component of the Reynolds stress tensor at the streamwise location of \(x/\theta = 50\), i.e., on the non-actuated surface at \(Re_\theta = 1033\), is shown in Fig. 3.

Fig. 3
figure 3

Comparison of the wall-normal distributions of the symmetric stresses and the shear-stress component of the Reynolds stress tensor on the coarse, standard, and fine grids for the non-actuated reference case and the actuated case \(N_{24}\) with DNS data by Schlatter and Örlü (2010)

The distributions of the standard and the fine mesh nearly collapse for the non-actuated reference case and for the actuated case \(N_{24}\) for all four components. Note that similar results were determined for other actuation parameter configurations. Only on the coarse mesh larger deviations are obtained. Furthermore, the data of the non-actuated reference case for the standard and the fine mesh shows good agreement with DNS data (Schlatter and Örlü 2010) of a turbulent boundary layer at a similar Reynolds number, i.e., \(Re_\theta = 1006\).

In conclusion, the analysis of the data shows that the resolution of the standard grid can be considered sufficient to accurately predict actuated turbulent boundary layer flow.

4.2 Wall-Shear Stress Reductions

The skin-friction reduction in the streamwise direction \(\varDelta c_f\) is defined in percent by

$$\begin{aligned} \varDelta c_f = \frac{\tau _{w,x,\mathrm {ref}} - \tau _{w,x,\mathrm {act}}}{\tau _{w,x,\mathrm {ref}}} \cdot 100\>, \end{aligned}$$
(3)

where the wall-shear stress is computed by

$$\begin{aligned} \tau _{w,x} = \left[ \left( \overline{\overline{\tau }}\cdot \mathbf {t}_{1,x}\right) \cdot \mathbf {n}\right] _{\mathrm {wall}}\,. \end{aligned}$$
(4)

The quantity \(\mathbf {t}_{1,x}\) denotes the unit tangential vector on the surface in the streamwise direction, i.e., \(\mathbf {t}_{1,x}\) has a zero z-component, \(\mathbf {n}\) is the unit normal vector of the surface, and \(\overline{\overline{\tau }}\) is the stress tensor. Note that the total wall-shear stress vector consists also of a second component which corresponds to the tangential vector in the spanwise direction, generating a net force in the spanwise direction which is considered in detail in Sects. 4.3 and 4.7. The wall-shear stress is averaged over the shaded surface \(A_\mathrm {surf}\) in Fig. 1 and the values for \(\varDelta c_f\) of the 80 cases are listed in Table 2 in the “Appendix”. The dependence of \(\varDelta c_f\) on the various parameters, i.e., the wavelength, period, amplitude, and amplitude velocity, is summarized in Fig. 4. The highlighted and numbered distributions are mainly from cases of the upper envelope of the wall-shear stress reduction, i.e., the maximum \(\varDelta c_f\) values for the wavelength, the period, and the amplitude are emphasized. The discussion and illustration in Fig. 4 summarize the pronounced varying dependence of the wall-shear stress on the different actuation parameters.

Fig. 4
figure 4

Dependence of the skin-friction reduction \(\varDelta c_f\) on a the wavelength \(\lambda ^+\), b the period \(T^+\), c the amplitude \(A^+\), and d the actuation velocity \(V^+ = 2\pi A^+/T^+\). For clarity, only the cases that define the upper envelope of the skin-friction reduction \(\varDelta c_f\) are shown in (a-c)

In Fig. 4a a quasi-linear increase of the skin-friction reduction \(\varDelta c_f\) as a function of the wavelength \( \lambda ^+\) is observed. Note, however, that this quasi-linear distribution is achieved by changing simultaneously the amplitude \(A^+\) and the period \(T^+\). Especially the latter has to undergo quite a non-linear variation to obtain such an approximately linear \(\varDelta c_f\) growth.

The dependence of \(\varDelta c_f\) on the wave period \(T^+\) at various \(A^+\) and \(\lambda ^+\) is presented in Fig. 4b. Due to the coupling between the forcing strength and the actuation period, which is in contrast to other actuation methods like spanwise traveling waves of spanwise forcing (Du et al. 2002) and traveling waves of flexible wall (Zhao et al. 2004), the optimum period \(T^+\) is determined by an internal, i.e., fluid mechanical, and an external, i.e., actuator, related condition. The ideal \(T^+\) is defined by the streak formation time scale (Touber and Leschziner 2012), i.e., for oscillatory spanwise forcing the period must be small enough to disrupt the reorganization of the streaks, and a sufficient strength of the forcing, which increases with decreasing period, is required. The dependence of the skin-friction reduction on the period in Fig. 4b shows that the optimum period among all 80 cases is on the order of \(T^+ = \mathcal {O}\left( 50 \right) \), which is slightly lower than the optimum period \(T^+ \approx 70\) of a spanwise oscillating wall in turbulent boundary layer flow (Lardeau and Leschziner 2013).

Note that likewise tendencies can be found for spanwise traveling oscillatory forcing (Du et al. 2002) such as increased drag reduction with higher wavelengths (cf. Fig. 4a). The longest wavelength considered in this study (\(\lambda ^+ = 3000\)) is comparable to that used in the experimental setups by Tamano and Itoh (2012) and Li et al. (2018). Although their lowest investigated period is \(T^+ \approx 110\) and thus considerably higher than the optimum found in this study, their results corroborate the tendency of higher wall-shear stress reduction at lower periods in the regime \(110 \le T^+ \le 302.5\).

The distribution of the skin-friction decrease as a function of the amplitude in Fig. 4c shows that the maximum skin-friction reduction is directly coupled to the amplitude. This is to some extent expected since the velocity and thus, the strength of the actuation \(V^+ = 2\pi A^+ / T^+\) is directly related to the amplitude. That is, at a given period the strength of the actuation is determined by the amplitude. This is confirmed by the experimental findings of Li et al. (2018). They obtain in a lower amplitude range a monotonic increase of the skin-friction reduction for increasing amplitude. Figure 4d presents the skin-friction reductions as a function of the velocity amplitude of the actuation \(V^+\), where the scaling shows a quasi-linear behavior for larger wavelengths \(\lambda ^+ > 1000\).

A similar scaling for \(\varDelta c_f\) was proposed by Tomiyama and Fukagata (2013) by combining the amplitude of the actuation velocity \(V^+\) (cf. Fig. 4d) and the thickness of the Stokes layer \(\sqrt{T^+/(2 \pi )}\) such that \(\varDelta c_f = f(A^+\sqrt{2\pi /T^+})\) which is plotted in Fig. 5. For shorter wavelengths (cf. Fig. 5a), i.e., \(\lambda ^+ \le 1000\), a linear scaling is only observed for small values \(A^+\sqrt{2\pi / T^+} < 10\). Note that the current overall reductions are lower than those in Tomiyama and Fukagata (2013) which is likely due to the higher Reynolds number in this study (\(Re_\tau = 360\) versus \(Re_\tau = 180\) in Tomiyama and Fukagata 2013) and due to a generally lower skin-friction reduction efficiency in turbulent boundary layers compared to turbulent channel flow (Ricco 2004). For higher scaling factor values, the distribution is more scattered. For larger wavelengths \(\lambda ^+ > 1000\) (cf. Fig. 5b), the skin-friction reduction scales almost linearly over the entire range. Above a certain value \(A^+\sqrt{2\pi / T^+} \gtrsim 20\), however, the linear behavior of the skin-friction reduction deteriorates. We believe the reason for this degradation is the large momentum injection into the boundary layer via too high a velocity amplitude. This increases the spanwise velocity component which leads to an amplified turbulent exchange. This deterioration is generated by flow separation occurring for the flow in spanwise direction. The separation is observed for smaller values of the parameter \(A^+\sqrt{2\pi / T^+}\) and for large amplitude to wavelength ratios \(A^+ / \lambda ^+\). This is further discussed in Sect. 4.7 in which the generation of spanwise force and flow is addressed.

Fig. 5
figure 5

Relative skin-friction reduction of all cases as a function of the scaling parameter \(A^+\sqrt{2 \pi / T^+}\) for a\(\lambda ^+ \le 1000\) and b\(\lambda ^+ > 1000\)

4.3 Drag Reduction

Introducing a wave motion of the surface means that the area of the moving wall increases with the amplitude and the wavelength of the wave. The data in Table 2 in the “Appendix” shows that this change of the wetted surface \(\varDelta A_\mathrm {surf}\) can be quite substantial, especially at small wavelengths and high amplitudes. At high wavelengths, this variation becomes rather small. Since the friction drag is defined by the product of the wall-shear stress and the surface interacting with the fluid, the variation of the wetted surface has to be taken into account leading to a non-linear relation between wall-shear stress and friction drag reduction.

The averaged drag reduction is defined as

$$\begin{aligned} \varDelta c_d = \frac{c_{d,x,\mathrm {ref}} - c_{d,x,\mathrm {act}}}{c_{d,x,\mathrm {ref}}} \cdot 100 \end{aligned}$$
(5)

where \(c_{d,x}\) is the drag coefficient in the streamwise direction. Since spanwise traveling transversal surface waves also generate a net spanwise force it is reasonable to also define a spanwise drag coefficient \(c_{d,z}\) which is discussed in detail in Sect. 4.7. The two drag coefficients are given by

$$\begin{aligned} c_{d,x} = \frac{2 F_x}{\rho _\infty u_\infty ^2A_\mathrm {ref}} \quad \text {and}\quad c_{d,z} = \frac{2 F_z}{\rho _\infty u_\infty ^2A_\mathrm {ref}}\,, \end{aligned}$$
(6)

where \(F_x\) and \(F_z\) are the streamwise and spanwise components of the force vector \(\mathbf {F}\), which is computed by integrating pressure and friction forces over the shaded surface \(A_\mathrm {surf}\) in Fig. 1

$$\begin{aligned} \mathbf {F}&= \mathbf {F}_p + \mathbf {F}_f = \int _{A_\mathrm {surf}} -p\mathbf {n}dA + \int _{A_\mathrm {surf}} \overline{\overline{\tau }}\cdot \mathbf {n} dA \end{aligned}$$
(7)

and \(A_\mathrm {ref} = 1\) is the reference surface.

To verify the computation of the drag, the fluxes of the integral momentum equation are integrated over the faces of a volume extending from \(x = 50\) to \(x = 100\), i.e., the streamwise interval over which the drag is computed, and in the wall-normal direction from the wall to a point outside the boundary layer at \(y = 52\). For the highest drag reduction case \(N_{80}\), a difference of \(0.023\%\) was determined compared to the direct integration over the wall surface, i.e., Eq. 7. This confirms the validity of the numerical data discussed in this study.

The data in Table 2 in the “Appendix” evidences the differences between \(\varDelta c_f\) and \(\varDelta c_d\), especially at small wavelengths. The highest drag reduction is \(\varDelta c_d = 26\%\) for a wavelength of \(\lambda ^+ = 3000\), a period of \(T^+ = 50\), and an amplitude of \(A^+ = 78\). Due to the large wavelength, the increase of the wetted surface is only \(\varDelta A_\mathrm {surf} = 0.7\%\). The highest drag increase \(\varDelta c_d = -27\%\) with the highest corresponding skin-friction coefficient increase \(\varDelta c_f = -7\%\) is observed for \(\lambda ^+ = 200\), \(T^+ = 20\), and \(A^+ = 30\). As stated before, this configuration with a small wavelength suffers considerably from a drastic increase of the wetted surface \(\varDelta A_\mathrm {surf} = 19.4\%\).

The temporal distributions of the instantaneous drag of the ”best” and ”worst”, i.e., highest and lowest drag reduction \(N_{80}\)\(\varDelta c_d = 26\%\) and \(N_2\)\(\varDelta c_d = -27\%\), which is a massive drag increase, are compared as an example in Fig. 6,

Fig. 6
figure 6

Temporal evolution of the instantaneous drag reduction \(\varDelta c_d\) for the non-actuated reference case, the actuated case with the highest drag reduction (\(N_{80}\), Table 2 in the “Appendix”), and the actuated case with the highest drag increase (\(N_2\), Table 2 in the “Appendix”)

with the instantaneous drag of the reference case. The temporal fluctuations of the drag appear stronger for the ”worst” case. Note that due to the larger wavelength the drag of the ”best” case is numerically integrated over a three times larger spanwise extent. Due to the spanwise averaging, this leads to the temporally smoother distribution of the \(N_{80}\) (\(\varDelta c_d = 26\%\)) case compared to the \(N_2\) (\(\varDelta c_d = -27\%\)) case.

4.4 Drag Reduction Modeling and Sensitivity Analysis

In the following, the drag reduction \(\varDelta c_d\) is modeled as a function of the actuation parameters \(\lambda ^+\), \(T^+\), and \(A^+\). For this task, LES simulations provide only a sparse data set comprising 80 points. The modeling is further complicated by the fact that these points are far from regularly distributed in the solution space spanned by \(\lambda ^+\), \(T^+\), and \(A^+\). A dense coverage of the actuation space using expensive LES simulations is hardly feasible.

The modeling is performed using a powerful regression solver from machine learning: support vector regression (SVR) (Cortes and Vapnik 1995). The algorithm is chosen for its prediction accuracy and its smooth response distribution. SVR is a supervised learning algorithm that predicts the response \(\hat{y}\), based on the feature vector \(\varvec{x}\), from training data points \((\varvec{x}_m,y_m)\). It is formulated as

$$\begin{aligned} \hat{y}(\varvec{x}) = \mu + \sum \limits _{m=1}^M\omega _m K (\varvec{x},\varvec{x}_m) = \mu + \varvec{\omega }^\mathrm {T} \varvec{K} (\varvec{x} ), \quad m=1,\ldots ,M \end{aligned}$$
(8)

where \(\hat{y}\) is the modeled output, \(\varvec{\omega }\) is the weight vector, and \(\varvec{K} (\varvec{x} )\) designates the Gaussian kernel functions used in this study. Here, the features are the actuation parameters \(\lambda ^+\), \(T^+\), and \(A^+\), and the responses are the averaged drag reduction \(\varDelta c_d\). By construction, \(\hat{y}(\varvec{x}) \) asymptotes to \(\mu \) as the kernel decays to zero outside the training range. In other words, SVR, as most machine learning algorithms, loses accuracy beyond the training parameter range.

Note that due to the highly non-linear response behavior and the scarcity of data points at very low wavelengths, cases with wavelengths \(\lambda ^+ < 500\) are ignored during the modeling process. This range is also of little interest as the best drag reduction is found at larger wavelengths. This data exclusion effectively yields 71 data points instead of 80. Overfitting is prevented with a 5-fold cross-validation. The SVR model has a coefficient of determination of \(R^2 = 0.93\) indicating an excellent prediction accuracy.

The SVR model from the LES data is employed to visualize a continuous actuation response in the investigated parameter range of \(A^+ \in [0,78]\), \(T^+ \in [20,120]\), and \(\lambda ^+ \in [500,3000]\). Figure 7 shows the isosurfaces of three drag reduction levels: \(\varDelta c_d=15\), 20, and 25%.

Fig. 7
figure 7

Actuation response surface. The gray surfaces represent drag reduction at three levels: \(\varDelta c_d = 15\), 20, and \(25\%\). The thick black line illustrates the ridgeline. Its solid portion is interpolated using all LES data. The dotted ridgeline between points A and B extrapolates a better actuation response at given \(\lambda ^+\) for amplitudes \(A^+>78\), i.e. beyond the investigated parameter range. Point C corresponds to the actuation parameters of the LES simulation \(N_{80}\) with the largest drag reduction. The cloud of small dots indicates the parameter datapoints of all considered simulations

Within this parameter range, the best performance of 26.5% is achieved at \(\lambda ^+=3000\), \(T^+=38\) and \(A^+=78\), which is slightly higher than the best simulated LES configuration \(N_{80}\) with \(\varDelta c_d=26.3\)%. This location indicates that better drag reduction could be achieved by increasing amplitude and wavelength.

An extrapolation of better performance outside the investigated parameter range is obtained with a ridgeline. In every \(\lambda ^+=\text {const}\) plane, the drag reduction \(\varDelta c_d\) features a single maximum \((A_r^+,T_r^+)\) with respect to the actuation amplitude \(A^+\) and period \(T^+\). The curve of \((A^+,T^+, \lambda ^+)\) connecting all these \(\lambda ^+\)-dependent \(\varDelta c_d\) maxima is the ridgeline, which is illustrated as a thick black curve in Fig. 7. Variables on this ridgeline are denoted by the subscript ‘r’.

In the range \(\lambda ^+ \in [560,1865]\), this maximum is inside the modeled \(T^+\) and \(A^+\) data range. It is illustrated by the solid black curve. However, the ridgeline leaves this modeled data range through the exit point A at the top surface \(A^+ =78\) near \(\lambda ^+ \approx 1865\). The rigdeline is extrapolated outside the data range and depicted as dotted curve between points A and B. The extrapolation method is detailed in Fernex et al. (2019). Along the ridgeline, \(\varDelta c_d\) monotonously increases from 7% at \(\lambda ^+=560\) to the maximum of 27.9% at \(\lambda ^+=3000\) (point B). Note that this point is outside the current LES parameter range.

The ridgeline defines a ‘skeleton’ of the parametric behavior as illustrated in Fig. 8. Figure 8a shows its projection in the \(\lambda ^+-T^+\) and \(\lambda ^+-A^+\) planes. Like in Fig. 7, the dotted sections correspond to the extrapolated ridgeline between points A to B. The amplitude \(A_r^+\) and period \(T_r^+\) along the ridgeline monotonously increase with the wavelength \(\lambda ^+\). The period asymptotes rapidly towards 44. The amplitude continually increases with the wavelength although at a decreasing rate.

An intriguing physical insight about the drag-reduction mechanism is revealed in Fig. 8b complementing Fig. 5b. The relative drag reduction \(\varDelta c_{d}\) along the ridgeline is shown as a function of the scaling parameter proposed by Tomiyama and Fukagata (2013) based on a Stokes layer of a transverse wall oscillation. \(\varDelta c_{d}\) clearly exhibits a linear behavior along the ridgeline in the scaling parameter range between 15 and 30. Away from the ridgeline, the scaling shows scatter on the order of that observed in Fig. 5b.

Fig. 8
figure 8

a Projection of the ridgeline in the \(\lambda ^+\)\(T^+\) and \(\lambda ^+\)\(A^+\) planes. b Drag reduction along the ridgeline as a function of the scaling proposed by Tomiyama and Fukagata (2013). The solid circles marked A, B, and C correspond to the likewise marked points in Fig. 7

4.5 Turbulent Flow Statistics

In the following, the mean statistics of a drag reduced flow will be investigated in detail. For this analysis, the case with the highest drag reduction, i.e., the case \(N_{80}\), will be considered. If data from other cases is used, it is explicitly indicated. All presented wall-normal distributions are considered at the streamwise position \(x = 90\, \theta \), which is located in the actuated region. The actuated fully developed turbulent flow possesses a momentum based Reynolds number \(\mathrm{Re}_\theta = 1077\). The flow field of the actuated cases is phase averaged in the spanwise direction. Therefore, a triple decomposition (Hussain and Reynolds 1970) of the flow variables is used

$$\begin{aligned} \phi = \underbrace{\overline{\phi } + \tilde{\phi }}_{\langle \phi \rangle } + \phi '', \end{aligned}$$
(9)

where \(\overline{\phi } = f(x,y)\) is the temporal and spanwise average representing phase independent quantities, \(\tilde{\phi } = f(x,y,z,t)\) are periodic fluctuations generated by the actuation, i.e., the secondary flow field, \(\langle \phi \rangle = \overline{\phi } + \tilde{\phi }\) are phase averaged quantities, and \(\phi '' = f(x,y,z,t)\) are stochastic fluctuations. The periodic fluctuations \(\tilde{\phi }\) and therefore also the phase average \(\langle \phi \rangle \) are constant when considering a constant phase angle but time-dependent in the z-coordinate direction. The detailed description of the phase averaging for spanwise traveling transversal surface waves can be found in Tomiyama and Fukagata (2013). Spanwise averages are obtained along lines of constant distance from the wall, i.e., along the curved mesh lines. This calculation of the spanwise average suffers from some uncertainty for short wavelengths with high amplitudes, where the traveling wave massively intrudes into the boundary layer. For spanwise averages of long wavelengths as in the \(N_{80}\) case, where the local perturbation of the viscous sublayer and the buffer layer is less drastic, this problem does not occur.

A first overall impression of the impact of the wave actuation on the turbulent coherent structures is given in Fig. 9 by comparing contours of the \(\lambda _2^+\)-criterion Jeong and Hussain (1995) for the random velocity fluctuations \(u''_i\). It is evident that the total number of vortical structures in the near-wall region is significantly reduced for the actuated flow. Extended regions of little to hardly any structures occur in Fig. 9b.

Fig. 9
figure 9

Contours of the \(\lambda _2^+\)-criterion (Jeong and Hussain 1995) colored by the random velocity fluctuations \(u''^+\) in the near-wall region \(y^+ < 20\) of a the non-actuated reference case, b the actuated case with the highest drag reduction \(N_{80}\), i.e., \(\lambda ^+ = 3000\), \(T^+ = 50\), and \(A^+ = 78\), normalized by the friction velocity of the non-actuated reference case \(u_{\tau ,\mathrm {ref}}\), and c the actuated case \(N_{80}\) normalized by the actual reduced friction velocity \(u_{\tau ,\mathrm {act}}\)

A closer look evidences that unlike the structures of the non-actuated reference flow in Fig. 9a, the structures of the actuated flow are inclined to the left and right depending on the phase angle of the traveling wave. It will be discussed in Sect. 4.6 that this wave determined orientation of the flow structures is an important feature related to drag reduction (Touber and Leschziner 2012). When the lowered friction velocity of the actuated case \(u_{\tau ,\mathrm {act}}\) is used for the non-dimensionalization, the contours of \(\lambda _2^*\) depicted in Fig. 9c show a less significant reduction of the vortical structures. Furthermore, a less homogeneous distribution is apparent and the preferential inclination angle of the structures depending on the wave phase persists.

To highlight the influence of the actuation on the instantaneous flow field, Fig. 10 shows contours of the instantaneous random fluctuations of the velocity component in the streamwise direction \(u''\) in a region \(0< y^+ < 20\) above the wall.

Fig. 10
figure 10

Contours of the random streamwise velocity fluctuations for \(u''^+ = -3\) (blue) and \(u''^+ = 3\) (red) in the near-wall region \(0< y^+ < 20\) of a the non-actuated reference case, b the actuated case with the highest drag reduction \(N_{80}\), i.e., \(\lambda ^+ = 3000\), \(T^+ = 50\), and \(A^+ = 78\), normalized by the friction velocity of the non-actuated reference case \(u_{\tau ,\mathrm {ref}}\), and c the actuated case \(N_{80}\) normalized by the actual reduced friction velocity \(u_{\tau ,\mathrm {act}}\)

For the non-actuated flow, regions of localized high-speed and low-speed fluid, i.e., streaks, are illustrated whose length and width are on the order of \(\mathcal {O}(10^3)\) and \(\mathcal {O}(10^2)\) in inner units. The actuated flow field, non-dimensionalized by the friction velocity of the reference case \(u_{\tau ,\mathrm {ref}}\) (cf. Fig.  10b), shows much less pronounced regions of high- and low-speed fluid. That is, the distinctive structure of thin meandering streaks is considerably alleviated compared to the non-actuated reference case. The non-dimensionalization by the lowered friction velocity of the actuated case \(u_{\tau ,\mathrm {act}}\) shown in Fig. 10c results in a distribution with still a considerably weakened streak strength. Furthermore, the previously observed preferential alignment of the streaky structures depending on the phase angle becomes even more apparent.

The wall-normal distributions of the phase averaged streamwise velocity \(\langle u \rangle \) above the wave crest and in the wave trough and the mean velocity \(\overline{u}\) are shown in Fig. 11.

Fig. 11
figure 11

Wall-normal distributions of the phase averaged streamwise velocity \(\langle u \rangle \) above the crest and in the trough and the spanwise averaged mean velocity \(\overline{u}\) for the non-actuated reference case and the actuated case \(N_{80}\); in a the \(N_{80}\) distributions are non-dimensionalized by the friction velocity of the non-actuated reference case and in b the \(N_{80}\) distributions are non-dimensionalized by the friction velocity of the actuated case \(N_{80}\)

The scaling with the friction velocity of the non-actuated reference case \(u_{\tau ,\mathrm {ref}}\) in Fig. 11a illustrates the decrease of the velocity in the near-wall region. The wall-normal gradient at the wall is lowered, which results in drag reduction. Scaling the velocities with the friction velocity of the actuated wall \(u_{\tau ,\mathrm {act}}\) in Fig. 11b leads to an offset of the velocity profiles in the logarithmic region with respect to the non-actuated reference case by \(\varDelta B^* \approx 3.8\). Based on the idea of the analysis of the impact of roughness on fully turbulent flow (Nikuradse 1933; Clauser 1956), Gatti and Quadrio (2016) suggested the offset \(\varDelta B^*\) to predict drag reduction at higher Reynolds numbers

$$\begin{aligned} \varDelta B^* = \sqrt{\frac{2}{c_{f,0}}} \left[ \left( 1-\varDelta c_f \right) ^{-1/2} - 1\right] - \frac{1}{\kappa }\ln \left( \frac{\mathrm{Re}_\tau }{\mathrm{Re}_{\tau ,0}}\right) \,. \end{aligned}$$
(10)

Note, however, that this equation cannot be further simplified since for actuated turbulent boundary layer flow the term \(\frac{\mathrm{Re}_\tau }{\mathrm{Re}_{\tau ,0}}\) is neither constant, as for constant pressure gradient turbulent channel flow, nor can it be substituted by the drag reduction rate, as for constant flow rate turbulent channel flow. Thus, \(\varDelta c_f\) cannot be directly determined by equation (10). Nevertheless, using the local values of \(c_{f,0}\), \(\varDelta c_f\), \(\mathrm{Re}_\tau \), and \(\mathrm{Re}_{\tau ,0}\) at \(x = 90\theta \) the calculated offset from Eq. (10) is \(\varDelta B^* = 4.07\), which reasonably agrees with the result \(\varDelta B^* = 3.8\) shown in Fig. 11b. The velocity profiles in Fig. 11 show that for the current actuation neither the non-actuated nor the actuated friction velocity scaling—regardless from crest, trough, or spanwise averaged wall shear scaling—result in a collapsed distribution over the entire boundary layer.

Next, the components of the Reynolds stress tensor are depicted in Fig. 12.

Fig. 12
figure 12

Wall-normal distributions of the symmetric components \(\overline{u''_iu''_i}\) and shear stress \(\overline{u''v''}\) of the Reynolds stress tensor for the non-actuated reference case and the actuated case \(N_{80}\) scaled by the friction velocity of a the non-actuated reference \(u_{\tau ,\mathrm {ref}}\) and b the actuated case \(u_{\tau ,\mathrm {act}}\). Spanwise averaged values are shown as lines and the shaded regions illustrate phase variations of the depicted quantity

If scaled by the friction velocity of the non-actuated reference case, the symmetric Reynolds stresses \(\overline{u''_iu''_i}\) and the Reynolds shear stress \(\overline{u''v''}\) shown in Fig. 12a are significantly lowered with only minor phase variations. For the case \(N_{80}\), the reductions at \(y^+ = 14.2\), which defines the location of the peak of the streamwise fluctuations and the location of the maximum streamwise velocity streak intensity, are approx. \(39\%\) for the streamwise component and \(62\%\) for the shear-stress component. This suggests that the turbulent motion close to the wall is massively damped. Touber and Leschziner (2012) have reported similarly large reductions in this region for spanwise wall oscillations without normal deflection. They emphasize the importance of the reduced near-wall Reynolds shear stress and drag, as characterized by the Fukagata, Iwamoto, Kasagi (FIK) identity (Fukagata et al. 2002), i.e., for the shear-stress contribution \(c_{f,\mathrm {RSS}} \sim \int _0^{\delta _{99}} (1-y)(-\overline{u''v''})dy\). If the friction velocity of the actuated case is used for non-dimensionalization (cf. Fig. 12b), the distributions are much closer to those of the non-actuated reference case. However, especially in the near-wall region a considerable difference between the distributions of the streamwise and shear-stress components remains. This corroborates the observations in the illustration of the streak contours in Fig. 10. That is, no self-similarity of the drag reduced boundary layer with the reference boundary layer can be determined.

Note that the stronger suppression of the streamwise fluctuations compared to the other components leads to a shift of the actuated distribution from one-dimensional to isotropic turbulence. This was also found for spanwise oscillations in Touber and Leschziner (2012).

The distributions of the joint probability density function (PDF) of the streamwise and the wall-normal stochastic velocity fluctuations \(u''\) and \(v''\) are presented in Fig. 13.

Fig. 13
figure 13

Joint PDF of \(u''\) and \(v''\) at \(y^+ = 12\) for the non-actuated reference case and the actuated case \(N_{80}\)

High values in the upper left quadrant (negative \(u''\) and positive \(v''\)) denote ejections of fluid from the near-wall region towards the outer flow, whereas high values in the lower right quadrant (positive \(u''\) and negative \(v''\)) indicate sweeps of high-speed fluid from the outer flow towards the near-wall region. As can be seen in Fig. 13, an overall attenuation of the fluctuations is observed with a strong damping of the sweeps and ejections in the second and the fourth quadrant. Again, this is in agreement with the results from spanwise oscillating walls without normal deflection (Agostini et al. 2014; Touber and Leschziner 2012).

Spanwise premultiplied energy spectra of the velocity fluctuations \(\kappa E_{u_i''u_i''}\), where \(\kappa = 2\pi / l_z\) is the wavenumber, are presented in Fig. 14.

Fig. 14
figure 14

Spanwise premultiplied energy spectra of the velocity fluctuations \(u_i''u_i''\) normalized by the total resolved kinetic energy of the related case for (left) the non-actuated reference case and (right) the actuated case with the highest drag reduction \(N_{80}\); a, b streamwise, c, d wall-normal, and e, f spanwise velocity component

Each spectrum is normalized by the total resolved kinetic energy of the related case, i.e., the non-actuated reference case and \(N_{80}\). No general decrease in the energy peak of the actuated case is observed, only a shift in the energy distribution as a function of the structural wavelength and wall-normal coordinate can be seen. A comparison between the two differently normalized spectra for the streamwise component (cf. Fig. 14a, b) shows an energy decrease especially for the small scales and in the near-wall region. In other words, for the actuated case \(N_{80}\) the energy is accumulated further off the wall in the larger scale turbulent structures. The peak of the non-actuated reference case at \(\lambda ^+ \approx 100\), which is associated with the typical spacing of the near-wall streaks of \(l_z^+ = \mathcal {O}(100)\), becomes less pronounced for the actuated case \(N_{80}\) and is shifted off the wall. This observation corroborates the visual impression from Fig. 10 of a strong reduction of the near-wall streaks for the actuated case. Similar tendencies are observed for the wall-normal (cf. Fig. 14c, d) and the spanwise velocity component (cf. Fig. 14e, f). Additionally, a stronger concentration of the energy in the length scale range of the near-wall streaks is observed for the spanwise velocity component. There is a sharper peak for the actuated case in comparison to a broader energy distribution in the non-actuated case.

In Fig. 15 the phase averaged and spanwise averaged distributions of the vorticity fluctuations \(\omega _i'' \omega _i''\) are presented as a function of the wall-normal distance. The comparison of the profiles of each component shows that the major attenuation is observed in the wall-normal and the spanwise components, whereas the streamwise component shows only minor changes. Generally, for all cases with \(\lambda ^+ > 1000\) a good correlation between the decrease of the skin-friction and the decrease of the peak of the wall-normal (\(R = 0.96\)) and spanwise (\(R = 0.98\)) vorticity fluctuations was found.

Fig. 15
figure 15

Wall-normal distributions of the phase and spanwise averaged vorticity fluctuations \(\overline{\omega ''_i\omega ''_i}\) for the non-actuated reference case and the actuated case \(N_{80}\); a streamwise, b wall-normal, and c spanwise component

Again, similar vorticity trends were reported for spanwise oscillating walls (Touber and Leschziner 2012). The drag reduction was discussed to be caused by the weakening of velocity streaks near \(y^+ \approx 10\) (Agostini et al. 2014). That is, at least for actuation with large wavelengths \(\lambda ^+ > 1000\), a direct interference with quasi-streamwise vortices (Tomiyama and Fukagata 2013) has a minor effect on drag reductions.

The comparison of the vorticity fluctuation contours for four cases, i.e., the non-actuated reference case, the case with the highest drag increase \(N_2\), a case with moderate drag reduction \(N_{24}\), and the case with the highest drag reduction \(N_{80}\), in Fig. 16 shows details about the phase dependence of the overall structure of the vorticity field. It is obvious that the drag increase is associated with strongly enlarged and highly phase dependent vorticity contours. Due to the high amplitude and short wavelength sinusoidal wall motion, the boundary layer flow is massively perturbed. This is completely different for the two drag reduction cases, where the overall boundary layer structure is maintained but with reduced values of the wall-normal and spanwise vorticity component. Phase variations occur especially for the wall-normal component with the highest decrease above the wave crest and the lowest decrease in the wave trough. Note, however, that despite the clear phase variations, the overall lowered vorticity levels are maintained throughout the entire actuation period, as shown in Fig. 15.

Fig. 16
figure 16

Contours of the vorticity fluctuations \(\overline{\omega ''_i\omega ''_i}\) in the y-z-plane at \(x = 90\,\theta \) for (first row ac) the non-actuated reference case, (second row df) the case with the highest drag increase \(N_{2}\), (third row gi) a case with moderate drag reduction \(N_{24}\), and (fourth row jl) the case with the highest drag reduction \(N_{80}\); (left) streamwise, (center) wall-normal, and (right) spanwise vorticity component. Note that case \(N_{80}\) has a wavelength of \(\lambda ^+ = 3000\), the images for this case are thus compressed for lack of space

4.6 Spanwise Shear

To investigate the secondary flow strength and its effect on the near-wall turbulent structures, the Stokes strain \(\partial \tilde{w}^+ / \partial y^+\), i.e., the rate of change in the wall-normal direction of the periodic fluctuations of the spanwise velocity component, is considered in Fig. 17 for cases with wavelengths \(\lambda ^+ = 200\), \(\lambda ^+ = 1000\), \(\lambda ^+ = 1800\), and \(\lambda ^+ = 3000\). Based on the data summarized in Table 2 in the “Appendix”, for each \(\lambda ^+ = \text {const}\) set the cases with the highest, medium, and lowest skin-friction reduction are shown. The Stokes strain is used to characterize the Stokes layer that develops above an oscillating wall without any wall-normal deflection. However, similar to a configuration with pure spanwise oscillating walls (Jung et al. 1992; Touber and Leschziner 2012) a Stokes-like layer is also generated by a transversal wave motion of the surface. Through the introduction of a periodic wall-normal velocity \(\tilde{v}\), a periodic spanwise velocity component \(\tilde{w}\) is induced via mass conservation resulting in a wall-normal shear distribution defining a Stokes layer. Cases with a high drag reduction generally show a high level of symmetric, i.e., positive and negative, spanwise shear very close to the wall, whereas less symmetric shear distributions yield lower drag reduction. When the Stokes drag significantly increases near the wall, i.e., a singular-like distribution occurs, the drag reduction reduces drastically. This observation supports the assumption that the drag reduction mechanism is strongly linked to spanwise oscillations which are generated either by wave oscillations (Touber and Leschziner 2012), traveling waves of spanwise forcing (Du et al. 2002), spanwise velocity at the wall (Zhao et al. 2004), or spanwise transversal surface waves (Klumpp et al. 2010b; Koh et al. 2015b).

Fig. 17
figure 17

Phase averaged spanwise shear above the wave crest (—) and in the wave trough (\(\ldots\)) for wavelengths a\(\lambda ^+ = 200\), b\(\lambda ^+ = 1000\), c\(\lambda ^+ = 1800\), and d\(\lambda ^+ = 3000\). The cases listed in Table 2 in the “Appendix” are representative for high (\(N_8\), \(N_{24}\), \(N_{70}\), \(N_{80}\)), medium (\(N_7\), \(N_{29}\), \(N_{75}\), \(N_{82}\)), and low (\(N_{20}\), \(N_{68}\), \(N_{83}\)) skin-friction reduction or skin-friction increase (\(N_{2}\)) at each wavelength

The assumption of the importance of the oscillating spanwise shear for skin-friction reduction also yields an explanation for the increasing skin-friction reduction with growing wavelength, which agrees with an observation of Du et al. (2002). For short wavelengths, e.g., \(\lambda ^+ = 200\), the integration of the spanwise shear distribution over the spanwise width of a near-wall streak, i.e., \(l_z^+ = \mathcal {O}(100)\), results in only a minor force in the spanwise direction acting on the streaks. For wavelengths \(\lambda ^+ > 1000\), however, the spanwise shear varies only slightly over the width of a streak, therefore a considerable spanwise force is determined by the integration over the spanwise streak width. Note that this behavior does not occur for spanwise wall oscillations, since the periodic spanwise shear does not depend on the spanwise coordinate.

A comparison of spanwise and streamwise shear is shown in Fig. 18.

Fig. 18
figure 18

Wall-normal distributions of a the phase averaged streamwise and spanwise shear for the actuated case \(N_{80}\) and b the ratio of the phase averaged spanwise and streamwise shear for cases with high (\(N_{80}\)), medium (\(N_{82}\)), and low (\(N_{83}\)) drag reduction

Touber and Leschziner (2012) discuss a certain optimal scenario for oscillatory forcing in turbulent channel flow, where the ratio of spanwise to streamwise shear reaches values of up to \(\frac{\partial \tilde{w}^+ / \partial y^+}{\partial u^+ / \partial y^+} \approx 3\). Figure 18b shows that a similar value of this ratio is obtained for the case with the highest skin-friction reduction \(N_{80}\), whereas a lower ratio is obtained for the cases with medium (\(N_{82}\)) and low (\(N_{83}\)) skin-friction reduction. For all cases, the change of the skin-friction is well correlated with a spanwise-to-streamwise shear ratio of \(\frac{\partial \tilde{w}^+ / \partial y^+}{\partial u^+ / \partial y^+} = 3.1\).

Overall, the results for the spanwise traveling transversal waves in Fig. 17 underline the similarities to other drag reduction techniques based on periodic spanwise forcing. The results of the cases with lower wavelength in combination with high amplitude and high frequency deviate from this observation due to the increased wetted surface.

4.7 Spanwise Flow

Unlike symmetric actuation techniques such as a spanwise oscillating wall spanwise traveling transversal surface waves generate a net spanwise force and flow. It can be observed in Fig. 19 that this force, when compared to the streamwise drag of the non-actuated reference case, can be on the same order of magnitude and reach values up to \(c_{d,z}/c_{d,x}\cdot 100 = 60\%\) for case \(N_2\). It is evident from Fig. 19b, c that the friction and the pressure drag have opposite signs such that they partially cancel each other, as is shown in Fig. 19a. The cancellation is also observed for streamwise traveling transversal waves traveling at speeds higher than the mean flow velocity (Shen et al. 2003).

Fig. 19
figure 19

Ratio of the spanwise to the streamwise drag coefficient \(c_{d,z}/c_{d,x}\cdot 100\) for a the total drag, b the friction drag, c the pressure drag, and d the pressure drag compensated by the number of wavelengths in the domain \(\lambda ^+ / L_z^+\) and the actuation amplitude \(A^+\)

The distributions of the pressure drag show a much steeper gradient for smaller wavelengths at increasing values of the actuation velocity amplitude \(V^+\) in comparison to large wavelength setups. That is, for similar values of the velocity amplitude \(V^+\) a much lower pressure drag is obtained for a long wavelength setup than for a short wavelength setup. However, this is also due to a larger projected surface of the wave in the xy-plane. When this varying projected surface is compensated by dividing by \(A^+\) and the number of wavelengths in the computational domain \(\lambda ^+ / L_z^+\), the result depicted in Fig. 19d delivers a straight view on the net pressure force over one wavelength in the z-direction. Still, the longest wavelengths with \(\lambda ^+ = 3000\) show the smallest value of the compensated pressure drag in comparison to much larger values for smaller wavelengths for the same actuation velocity amplitude \(V^+\). The difference in the pressure drag is due to a spanwise flow which might cause separation on the leeward side of the wave. For instance, this is observed for case \(N_2\). The separation which is apparent in Fig. 16e, f leads to a drastic increase of the pressure drag. This explains the earlier breakdown of the linear scaling of the skin-friction reduction with the parameter \(A^+ \sqrt{2\pi /T^+}\) in Fig. 5a, b for smaller compared to larger wavelengths.

To investigate the net spanwise flow, the wall-normal distributions of the mean spanwise velocity for a few relevant cases are depicted in Fig. 20.

Fig. 20
figure 20

Wall-normal distribution of the mean velocity in the spanwise direction for several actuated cases

For the case with the strongest overall spanwise force \(N_2\) and the case with the strongest pressure force \(N_{20}\), a considerable net spanwise flow is generated. It is on the same order of magnitude in the near-wall region as the streamwise flow (cf. Fig. 11a). In contrast, the case with the highest drag reduction \(N_{80}\) and the case with the highest net energy saving \(N_{84}\) show a weaker spanwise flow. That is, naturally a positive net energy saving is only achieved if the flow in the spanwise direction, which needs to be transported by additional power input, is small.

4.8 Energy Saving Analysis

The previous discussion has shown that considerable drag reduction rates have been obtained. However, drag reduction is not the only metric of interest. From a prospective application point of view, the question of net energy saving and its relation to drag reduction must be addressed. The ideal net energy saving is defined as

$$\begin{aligned} \varDelta P_{\mathrm {net}} = \frac{P_{d,\mathrm {ref}} - (P_{d,\mathrm {act}} + P_{\mathrm {control}})}{P_{d,\mathrm {ref}}} \cdot 100\,, \end{aligned}$$
(11)

where \(P_{d,\mathrm {ref}/\mathrm {act}} = u_\infty (F_{f,x} + F_{p,x})\) is the power necessary to overcome the friction \(F_{f,x}\) and pressure forces \(F_{p,x}\) of the non-actuated \(P_{d,\mathrm {ref}}\) and actuated surface \(P_{d,\mathrm {act}}\) in the streamwise direction. The power spent on control \(P_{\mathrm {control}}\), i.e., on deflecting the surface for the traveling wave, is computed by an integration over the actuated surface \(A_\mathrm {surf}\)

$$\begin{aligned} P_\mathrm {control} = \int _{A_\mathrm {surf}} \left[ -p\mathbf {v} + \overline{\overline{\tau }}\cdot \mathbf {v} \right] \cdot \mathbf {n}dA\,. \end{aligned}$$
(12)

hence \(P_{\mathrm {control}}\) is a sum of the pressure and friction forces multiplied by the speed of the wall motion. Note that the same definition was also used in Shapiro et al. (1969) and Floryan and Zandi (2019). The values of \(P_{\mathrm {control}}\) and \(\varDelta P_\mathrm {net}\) are depicted in Fig. 21 for all cases.

Fig. 21
figure 21

Dependence of a the power spent \(P_{\mathrm {control}}\) and b the net power saving \(\varDelta P_{\mathrm {net}}\) on the cube of the actuation velocity amplitude \((V^+)^3 = (2\pi A^+/T^+)^3\); a zoom of the red rectangle in (b) is shown in (c). To indicate which cases possess net power saving the notation for three selected cases is given in (c)

The data for \(\varDelta P_\mathrm {net}\) are also listed in Table 2 in the “Appendix”. Fig. 21a shows the expected approximately linear dependence between the power spent \(P_\mathrm {control}\) and the actuation velocity cubed \((V^+)^3 = \left( 2\pi A^+/T^+\right) ^3\).

It is evident from Fig. 21b that net power saving is only obtained for a few cases with a maximum of \(\varDelta P_{\mathrm {net}} = 10\%\) for case \(N_{84}\). Most cases clearly show no net power saving but net power loss. For instance, for \(N_{20}\)\(\varDelta P_{\mathrm {net}}\) is \(\varDelta P_{\mathrm {net}} = -289\%\), i.e., almost the fourfold \(P_{d,\mathrm {ref}}\) has to be invested.

A closer look at the net-power-saving cases in Fig. 21c shows that high drag reduction rates are no indicator for high net power saving. That is, there is no linear relation between drag reduction and net power saving. Instead, a high value of the scaling parameter \(A^+ \sqrt{2\pi / T^+}\) obtained by a low amplitude speed \(2 \pi A^+ / T^+\) leads to positive net power saving. Thus, as expected, there is a trade-off between the minimum power input to effectively influence the turbulent boundary layer and a maximum power input above which the energy costs grow tremendously. In the current parameter range, the optimum energy saving solution, i.e., \(\varDelta P_\mathrm {net} = 10\%\), is achieved for the \(N_{84}\) case with \(\lambda ^+ = 3000\), \(T^+ = 90\), and \(A^+ = 66\), which possesses just a medium drag reduction of \(\varDelta c_d = 16\%\). Note that the parameters that result in high net energy saving are in the upper range of the interval. Furthermore, the data in Table 2 in the “Appendix” indicates that the sensitivity of \(\varDelta P_\mathrm {net}\) is less pronounced for larger wavelength and above a wave period of 60.

5 Conclusions

To analyze drag reducing effects and the net energy saving potential of spanwise traveling transversal surface waves high-resolution large-eddy simulations were conducted. The parameter space defined by the wave amplitude, wave period, and wavelength was investigated based on 80 wave parameter setups for purely spanwise traveling waves. The variation of skin-friction reduction, i.e., mean wall-shear stress alteration, drag reduction, i.e., surface integrated wall-shear stress, and net energy saving was analyzed. In brief, a maximum drag reduction and net energy saving of \(26\) and \(10\%\) was found.

The highest skin-friction reduction was achieved for a period of \(T^+ \approx 50\), which is lower than the one reported for spanwise oscillating wall and within the range of the streak formation time scale. Larger wavelengths and amplitudes yielded higher skin-friction reduction. For wavelengths larger than 1000 plus units, a scaling with the Stokes layer height and the velocity amplitude was found to predict skin-friction reduction reasonably well. Additionally, the difference between skin-friction reduction and drag reduction, i.e., the increase of the wetted surface was taken into account, was found to be substantial for short wavelengths in combination with high amplitudes. A drag-reduction model was derived from the sparse dataset using optimized support vector regression. From the model, a tendency to an asymptotic behavior of amplitude and period could be identified, supporting the assumption of an optimum period in the range \(40 \le T^+ \le 50\) for large wavelengths. Moreover, a ridgeline behavior of optimum drag reduction in the high wavelength regime was extracted from the model.

The statistical results of the turbulent flow field confirmed this result for high wavelength configurations, where similar effects of the actuation on the near-wall region compared to spanwise oscillating walls were observed. That is, considerable reductions of the near-wall velocity streak strength were found for the cases with high drag reduction. For the highest drag reduction case, the smaller wall-shear stress was coupled to a substantial decrease of the Reynolds shear stress in the near-wall region. Generally, for large wavelength cases \(\lambda ^+ > 1000\) the decrease of the wall-normal and spanwise vorticity fluctuations strongly correlated with skin-friction reduction and drag reduction. A comparison among several configurations revealed that for unfavorable combinations of short wavelength and high amplitude, a considerable increase of the turbulent exchange resulting in skin-friction and drag increase was observed, whereas large wavelengths circumvented this effect and led to drag reduction. The periodic secondary flow field generated by the wavy surface motion approximated that of Stokes flow. Similar oscillating spanwise shear distributions were observed for many drag reducing cases, although no perfectly symmetrical oscillatory excitation of the near-wall structures is achieved. The generation of a net spanwise force and flow was shown where strong spanwise forces and considerable net spanwise flow are obtained for cases with drag increase. Furthermore, the strength of the spanwise pressure force depends on the wavelength which is explained by a lower susceptibility to separation for large wavelengths due to a smaller amplitude–to–wavelength ratio.

No linear relationship between drag reduction and net energy saving was determined. That is, due to the non-linear response of the near-wall flow to the actuation the highest drag reduction does not result in the highest net energy saving. The maximum net energy saving \(\varDelta P_\mathrm {net} = 10\%\) was achieved for a drag reduction of \(\varDelta c_d = 16\%\), which is clearly lower than the maximum drag reduction of \(\varDelta c_d = 26\%\). A high value of the product of the actuation amplitude speed and the thickness of the Stokes layer at low amplitude speed results in positive net energy saving. The susceptibility of \(\varDelta P_\mathrm {net}\) is less pronounced for larger wavelength, which is a promising observation for prospective applications.