3.1. Direct numerical simulations

The numerical code used for the DNS is NEK5000, an open-source massively parallel spectral element code. The spatial discretization is based on a local tensor product of Gauss–Lobatto–Legendre quadrature points within each hexahedral element. Pressure stabilization is obtained by a $P_{N}{-}P_{N-2}$ formulation (in the present simulations $N\geqslant 9$ ). A third-order backward-differentiation formula in time is used for the time derivatives. Concerning the other terms of the equations, the viscous ones are treated implicitly while the nonlinear convective ones are forward extrapolated in time with third-order accuracy so as to treat them explicitly in the time-advancing scheme. Details on the code can be found following the link http://nek5000.mcs.anl.gov/index.php/Main_Page.

The simulations of the steady BL flow including the free-stream vortices were carried out with a computational domain starting at $X_{v}$ , which is a variable in the range $0\leqslant \mathit{Re}_{X_{v}}=X_{v}U_{\infty }/{\it\nu}\leqslant 2.74\times 10^{4}$ , where $\mathit{Re}_{X_{v}}$ is the Reynolds number at $X_{v}$ and ending at $X_{f}$ . The height of the computational domain is $H=86~\text{mm}$ , giving $H/{\it\delta}\geqslant 274$ at the inlet section where ${\it\delta}(x)=\sqrt{{\it\nu}x/U_{\infty }}$ is the characteristic length scale of the BL. The simulations are limited to computing only one pair of free-stream vortices, i.e. with a spanwise computational domain width of $1{\rm\Lambda}$ , but with the imposition of periodic boundary conditions in the spanwise direction. Two different streamwise extents of the computational domains were used: the first domain (D1) extends for a normalized length equal to ${\rm\Delta}\mathit{Re}_{X}=\mathit{Re}_{X_{f}}-\mathit{Re}_{X_{v}}\simeq 4.24\times 10^{5}$ and employs approximately $6\times 10^{7}$ degrees of freedom (dofs); the second one (D2), used for the configurations of slowly evolving streaks, has ${\rm\Delta}\mathit{Re}_{X}=\mathit{Re}_{X_{f}}-\mathit{Re}_{X_{v}}\simeq 1.00\times 10^{6}$ with approximately $10^{8}$ dofs.

The imposed velocity distribution at the inflow at $X_{v}$ consists of a superposition of the velocity field induced by a pair of vortices in the free stream to a Blasius BL profile. This implies that the axis of the vortices at $X_{v}$ is parallel with the plate. No-slip boundary conditions are used at the solid wall and null-stress outflow boundary conditions are imposed on the outflow and on the top boundary. The spatial resolution has been chosen on the basis of the simulations reported in Camarri et al. (Reference Camarri, Fransson, Talamelli, Talamelli, Oberlack and Peinke2013) and of dedicated convergence tests which are not shown here for the sake of brevity. Moreover, the results of the DNS have been further validated against those obtained by a three-dimensional BL solver and employing a completely different numerics (mixed finite-difference/finite-element discretization on unstructured grid).

Simulations of transition, triggered by forcing high-initial-amplitude TS waves, have also been performed in order to investigate the control effect by the free-stream vortices on transition delay. In these simulations the domain was extended to $10{\rm\Lambda}$ in the spanwise direction, thus including 10 pairs of counter-rotating vortices to avoid any influence of the boundary conditions on the transition location. In the streamwise direction the domain extends from $\mathit{Re}_{X_{v}}=5.28\times 10^{4}$ to $\mathit{Re}_{X_{f}}=6.11\times 10^{5}$ and the normalized forcing frequency of the TS wave is $F=10^{6}(2{\rm\pi}f{\it\nu}/U^{2})=140$ . The TS waves are excited in the DNS by a time-periodic volume force oriented in the wall-normal direction and acting in a narrow streamwise strip centred at the point $X_{TS}=150~\text{mm}$ . The forcing frequency matches the experimental C01 case and the wall-normal disturbance amplitude further downstream compares well with the corresponding eigenfunction of the Orr–Sommerfeld equation. The adopted boundary conditions are the same as for the simulation of the base flow. Concerning the spatial discretization, convergence tests have been carried out leading to a grid that is globally coarser than the one used for a single vortex pair and employs approximately $2\times 10^{8}$ dofs. This grid regeneration was carried out because of the different streamwise extension of the domain for the transition simulations and, at the same time, in order to keep the number of dofs for the transition simulation affordable from a computational point of view.

3.2. Biglobal stability analysis

The linear stability of the streaky BL induced by the free-stream vortices has been characterized by performing a local biglobal stability analysis, similar to the analysis by Piot, Casalis & Rist (Reference Piot, Casalis and Rist2008) on the flow past a periodic row of roughness elements. Since the flow is fully three-dimensional, the discretization of the stability problem has to be performed in the cross-sectional $yz$ -plane at each downstream location. The local stability problem is governed by the linearized Navier–Stokes equations for an incompressible flow and a generic solution of the streamwise ( $u$ ), wall-normal ( $v$ ) and spanwise ( $w$ ) velocity disturbances as well as the pressure perturbation ( $p$ ) is searched for in the following modal form: $\{u,v,w,p\}=\{\tilde{u} ,\tilde{v},\tilde{w},\tilde{p}\}(y,z)\text{e}^{\text{i}({\it\alpha}x-{\it\omega}t)}$ . Here, $\tilde{~}$ denotes the disturbance amplitude distribution in the $yz$ -plane, ${\it\alpha}$ is the streamwise wavenumber and ${\it\omega}$ is the angular frequency. The disturbance is assumed to vanish on and far away from the wall, and to be periodic in the spanwise direction (subharmonic modes have also been checked using multiples of the periodicity length of the base flow). When the linearized Navier–Stokes equations are discretized in space, a generalized eigenvalue problem is obtained, which is quadratic in ${\it\alpha}$ for the spatial analysis ( ${\it\omega}$ fixed) and linear in ${\it\omega}$ for the temporal analysis ( ${\it\alpha}$ fixed). Spatial discretization is here obtained by a finite-element formulation on unstructured grids, employing Taylor–Hood elements. The open-source code FreeFem $++$ (http://www.freefem.org) has been used to this purpose. The stability has been characterized by a spatial stability analysis. The initial guess for the solution of the resulting nonlinear eigenvalue problem has been estimated by performing a temporal stability analysis and applying the Gaster transformation to the results. Both the linear and nonlinear eigenvalue problems have been solved by a Krylov–Schur method with a shift-invert technique, using the parallel implementation available in the SLEPc library (see http://www.grycap.upv.es/slepc/).

4.1. Exploration of the free parameters

In the present flow configuration there are a number of free parameters, as detailed in § 2 and illustrated in figure 1. The geometrical and vortical parameters are ( $X_{v}$ , $h_{v}$ , $l_{v}$ , ${\rm\Lambda}$ ) and ( $r_{0}$ , ${\rm\Gamma}_{0}$ , $l_{p}$ ) respectively, with $l_{p}$ only being part of the Batchelor vortex model. As already discussed, we intend to demonstrate the use of free-stream vortices for BL flow control, and hence an optimization of all the above parameters is beyond the scope of the present study. However, since modelled vortices are introduced it is important to grant the physical realizability of the flow by varying the free parameters within physically meaningful ranges. These have been determined on the basis of previous experience with MVGs, as described below.

Physically relevant values of $r_{0}$ , ${\rm\Gamma}_{0}$ and $l_{v}$ have been found by analysing the vorticity generated past the MVGs in case C01 (see Shahinfar et al. Reference Shahinfar, Fransson, Sattarzadeh and Talamelli2013), which has been extracted from available DNS data (Camarri et al. Reference Camarri, Fransson, Talamelli, Talamelli, Oberlack and Peinke2013). The streamwise vorticity distribution inside one of the two counter-rotating vortices past an MVG pair is shown in figure 2(a) at a distance of 13 mm downstream of the MVG. The circulation has been calculated on circuits coinciding with the isocontours of axial vorticity (in figure 2 a) and is plotted in figure 2(b) versus the averaged distance $r$ of each isocontour from the centre of the vortex. Plausible values for ( $r_{0}$ , $l_{v}$ , ${\rm\Gamma}_{0}$ ) have been derived to be ( $0.8~\text{mm}$ , $3.12~\text{mm}$ , $9\times 10^{-4}~\text{m}^{2}~\text{s}^{-1}$ ) respectively from the synthetic data in figure 2(a,b). The parameters $l_{v}=3.12~\text{mm}$ and ${\rm\Lambda}=13~\text{mm}$ have been kept constant for all cases. We note that the only parameters that cannot be extracted on the basis of case C01 are $X_{v}$ and $h_{v}$ , since they are related to the placement of a hypothetical free-stream device generating the vortices. However, a constraint on $h_{v}$ is that we are interested in generating the vortices outside the BL and $X_{v}$ should be chosen such that an appreciable streak amplitude is obtained in the region of convective instability for the uncontrolled BL.

Figure 2. (a) Streamwise vorticity in the vortex 13 mm behind one MVG blade (shown in the background) for configuration C01 (isocontours equally distributed from 0 to $1.0\times 10^{-3}~\text{m}^{2}~\text{s}^{-1}$ ). (b) Vortex circulation, computed on the isocontours of vorticity, plotted against the area-weighted mean radius $r$ .

On the basis of the above considerations an exploratory set of simulations, varying some of the free parameters in the ranges $(0.65\leqslant r_{0}\leqslant 1.3)~\text{mm}$ , $(1.30\leqslant H\leqslant 7.80)~\text{mm}$ , $(0\leqslant X_{v}\leqslant 235.00)~\text{mm}$ and $(6\times 10^{-4}\leqslant {\rm\Gamma}_{0}\leqslant 27\times 10^{-4})~\text{m}^{2}~\text{s}^{-1}$ , were carried out in domain D1. A representative subset of the obtained results, using the Rankine model for the inlet conditions, is reported in figure 3, where the conventional streak amplitude measure,

(4.1) $$\begin{eqnarray}A_{ST}^{\mathit{max}}=\max _{y}\Big[\max _{z}\{U(X,y,z)\}-\min _{z}\{U(X,y,z)\}\Big]\Big/(2U_{\infty }),\end{eqnarray}$$

is plotted versus $x$ . Figure 3(a) illustrates the effect of $h_{v}$ on the generated streaks ( $X_{v}=26.0~\text{mm}$ ). The three values $h_{v}=1.95$ , 2.60, 3.90 mm correspond to $h_{v}/{\it\delta}(X_{v})=8.8$ , 11.7, 17.6 respectively. It should be recalled that the Blasius BL edge corresponds to approximately $5{\it\delta}$ . As $h_{v}$ increases, the global $A_{ST}^{\mathit{max}}$ maximum increases slightly, but with its corresponding downstream position appearing further downstream. This is an expected behaviour, since the further away from the BL edge the vortices are, the lower the dissipation rate becomes, and at the same time they induce a smaller modulation of the BL. The effect on the streak amplitude evolution by varying $X_{v}$ is rather weak, as illustrated in figure 3(b). It is worth noting that a variation in $X_{v}$ brings approximately a variation of $h_{v}/{\it\delta}(X_{v})$ for a constant $h_{v}$ , which may explain the variations among the cases in figure 3(b). In the same figure, two different ${\rm\Gamma}$ -values are also considered, namely $6\times 10^{-4}$ (case E1) and $9\times 10^{-4}~\text{m}^{2}~\text{s}^{-1}$ (case E2), plotted with blue and black lines respectively. As expected, the evolution of the streak amplitude is quite similar but with a lower amplitude for a lower ${\rm\Gamma}$ -value. From here on, we will focus on these two cases (E1 and E2) with the parameters E1: $X_{v}=52~\text{mm}$ , $r_{0}=1.3~\text{mm}$ , $h_{v}=3.9~\text{mm}$ (or $h_{v}/{\it\delta}(X_{v})=12.5$ ) and ${\rm\Gamma}_{0}=9\times 10^{-4}~\text{m}^{2}~\text{s}^{-1}$ ; E2: as E1 but with ${\rm\Gamma}_{0}=6\times 10^{-4}~\text{m}^{2}~\text{s}^{-1}$ . Case E2 has been designed to obtain a global $A_{ST}^{\mathit{max}}$ maximum approximately equal to that in C01, which has proven to successfully attenuate TS waves and accomplish transition delay.

Figure 3. Streamwise evolution of the streak amplitude on varying (a)  $h_{v}$ ( $X_{v}=26~\text{mm}$ ) and (b)  $X_{v}$ ( $h_{v}=3.90~\text{mm}$ ). The Rankine vortex model is used here.

Finally, as concerns the sensitivity to the actual vortex model, case E1 was also computed using the Batchelor model (case E1b). In this case the helical pitch $l_{p}$ was set to 2.6 mm, in agreement with the experimental results in Velte et al. (Reference Velte, Hansen and Okulov2008), where a similar configuration to ours was investigated. The results for E1b (discussed later in figure 5 a) show that the difference in streak amplitude with respect to the two different vortex models (Rankine and Batchelor) corresponds to a maximum of 2 percentage units in terms of the global $A_{ST}^{\mathit{max}}$ maximum. Thus, we can conclude that the streamwise amplitude evolution is rather insensitive to the selected vortex model.

4.2. Streamwise evolution of streaks

Three downstream locations of the streamwise velocity (colourmap) together with the streamwise vortices, identified using isocontours of the ${\it\lambda}_{2}$ vortex identification criterion, are shown in figure 4 for the E2 case. The three plots clearly show that the streamwise vortices successively modulate the BL in the streamwise direction. These vortices drive high-momentum fluid from the free stream towards the wall in between the two vortices of a pair, while low-momentum fluid is elevated in the region between two contiguous vortex pairs. Since $l_{v}<{\rm\Lambda}$ , the resulting high-velocity streaks are initially narrower and stronger than the low-velocity streaks. Moving downstream along the plate, the streamwise vortices progressively approach and eventually enter the BL.

Figure 4. Case E2: cross-sectional contours of the mean streamwise velocity field $U/U_{\infty }$ : (a)  $x=310~\text{mm}$ , (b)  $x=850~\text{mm}$ and (c)  $x=1675~\text{mm}$ . The dashed black lines are contours of constant velocity (0.1:0.1:0.9). The solid black lines represent the ${\it\lambda}_{2}$ vortex identification criterion.

In figure 5(a) the experimental $A_{ST}^{\mathit{max}}$ for case C01 is plotted and compared with the induced BL streak amplitude by the free-stream vortices (E1 and E2 cases). This shows that when MVGs are used, the BL streaks are formed markedly earlier since the vortices are intense and already inside the BL, but the streaks also dissipate more rapidly. When vortices equivalent to those past the MVGs are positioned outside the BL (case E1), they generate higher-amplitude streaks because the lift-up mechanism is active for a longer distance. This is quantified in figure 6(a), where a significantly more rapid decay of vortex circulation for the BL vortices is shown compared with the free-stream vortex decay of the E1 and E2 cases. Since the induced velocity modulations are initially weaker due to the longer distance from the wall, the generation process takes a longer streamwise distance to reach the same amplitude. Thus, to obtain a global $A_{ST}^{\mathit{max}}$ maximum equal to that in case C01, less intense vortices should be used, as in case E2. An additional comparison between the cases E2 and C01 is shown in figure 5(b), where the streamwise velocity and the streamwise vortices of E2 (left) and C01 (right) are reported at the point of global $A_{ST}^{\mathit{max}}$ maximum of both cases. The figure clearly highlights the differences between the two cases, especially in the shape of the velocity streaks and the different positions of the streamwise vortex within the BL.

Figure 5. (a) Streamwise streak amplitude evolution for the cases C01 (▫, from experiments), E1 (blue line), E1b (dashed line) and E2 (black line). (b) Cross-sectional contours of normalized mean streamwise velocity $U/U_{\infty }$ , for E2 (left) and for C01 (DNS, right) at the position of global streak amplitude maximum. Dashed contour lines correspond to (0.1:0.1:0.9). The solid black lines represent the ${\it\lambda}_{2}$ vortex identification criterion.

Figure 6. (a) Streamwise evolution of the circulation of the streamwise vortices for E1 (blue line), E2 (black line) and C01 (DNS, red line). (b) The solid blue line corresponds to the neutral curve for E1 and the dashed line to the neutral curve for the Blasius BL in the $F{-}\mathit{Re}$ -plane. Case E2 is stable in the range of $\mathit{Re}_{x}$ considered in the figure.

4.3. Spatial stability analysis

Here, we present the results of the biglobal spatial stability analysis of the cases E1 and E2. The results are condensed in the $F{-}\mathit{Re}$ stability diagram shown in figure 6(b), where the solid blue line enclosing the striped area corresponds to the E1 case and the dashed line to the neutral stability curve of the reference Blasius BL. For the E1 case the unstable region is enlarged, with respect to the Blasius BL, in terms of a broader frequency range, but with a downstream shift of the critical Reynolds number ( $Re_{x}^{cr}$ ) down to $Re_{x}^{cr}\approx 4\times 10^{5}$ , i.e. delayed by a factor of four. Inside the unstable region (striped area) the growth factors for both cases are comparable (not shown in the figure), indicating that the E1 streaky base flow case has the potential to accomplish transition delay with respect to the Blasius BL. The onset of instability of the E1 case starts around the streamwise position where the global $A_{ST}^{\mathit{max}}$ maximum is located (see figure 5). With the current knowledge on streamwise streak instability we suspect that the instability in the E1 case, with a global $A_{ST}^{\mathit{max}}$ maximum of approximately 29 % of $U_{\infty }$ , is driven by an instability of the streaks themselves, whose amplitude is too high.

For the E2 case, on the other hand, the modulated BL is stable for all the sections of the computational domain, i.e. the critical Reynolds number is greater than $Re_{x}^{cr}=8\times 10^{5}$ . With a global $A_{ST}^{\mathit{max}}$ maximum of less than 22 % of $U_{\infty }$ , it is close to the successful experimental case C01 and well below the critical $A_{ST}^{\mathit{max}}$ of 26 % of $U_{\infty }$ for the onset of sinusoidal secondary instability (see Andersson et al. Reference Andersson, Brandt, Bottaro and Henningson2001). This explains why in E2, where the streak amplitude is lower, the convective instability of the BL is completely stabilized throughout the length of the computational domain.

The results presented here demonstrate that there is a large margin of improvement for the stability of the controlled BL, since it is possible to modify the streak amplitude and its shape by acting on several parameters. For instance, the amplitude of the generated streaks can be controlled both by acting on ${\rm\Gamma}$ , as done in switching from E1 to E2, and by varying the initial distance $h_{v}$ of the vortices from the wall, so as to reduce the intensity of the induced velocities at the BL edge by the lift-up mechanism. Despite this the free parameters are not optimized, that is with the present configurations E1 and E2, the delay in the onset of convective instability inside the BL is so evident that a gain in terms of transition delay is expected by the proposed control.

4.4. BL transition delay

Previously reported experiments (Shahinfar et al. Reference Shahinfar, Fransson, Sattarzadeh and Talamelli2013) and now DNS simulations have been carried out for the configuration C01, where transition was triggered by forcing a TS wave at $F\simeq 140$ . This frequency was chosen in the experiment after having locked the location of the MVG array and successfully established a stable streaky BL. The desire was to locate branch I of the neutral stability curve of the reference Blasius BL without the MVGs at the location of the global $A_{ST}^{\mathit{max}}$ maximum, which has shown to be effective in attenuating TS waves (see Bagheri & Hanifi Reference Bagheri and Hanifi2007). In the DNS documented here the TS waves are forced as detailed in § 3 at $F=140.6$ . The forcing has been introduced in a DNS simulation of the reference Blasius BL, and the forcing amplitude has been varied and tuned so as to observe transition inside the computational domain. The result is shown in figure 7(a), where the temporal root mean square fluctuations of the normalized streamwise velocity ( $u_{rms}/U_{\infty }$ ) are plotted in the $xz$ -plane 1 mm above the wall. The figure shows that the initially forced amplitude at $X_{TS}=150~\text{mm}$ decays until it reaches branch I ( $x\approx 250~\text{mm}$ ), from where the $u_{rms}$ starts to grow, which it does beyond branch II ( $x\approx 600~\text{mm}$ ) since the flow is no longer linear with this high forcing amplitude. After some downstream location the $u_{rms}$ level decreases slightly, whereafter the flow undergoes a secondary instability, triggering transition after $x\simeq 850~\text{mm}$ . In a second DNS simulation the same forcing was applied as in the Blasius BL, but with the presence of the free-stream vortices (case E1). The result is presented in figure 7(b), showing that the disturbance is rapidly attenuated after an initial transient growth. The forced TS waves reach such a low value that it is not possible to detect any successive growth, at the scale of the plot, when they again enter the convectively unstable region, which is almost at the end of the computational domain (recall that $Re_{X_{f}}=6.11\times 10^{5}$ ).

Figure 7. Dimensionless temporal root mean square streamwise velocity fluctuations $u_{rms}/U_{\infty }$ in the $xz$ -plane at $y=1.0~\text{mm}$ (TS waves at $F\simeq 140$ ): (a) without free-stream vortices; (b) with free-stream vortices (case E1).

The test carried out in the present section shows that, even though it is the E1 case that is considered, which is characterized by the stability curve illustrated in figure 6(b), it is possible to obtain an evident transition delay with respect to the reference Blasius BL. A DNS with TS wave forcing of the E2 case has not been carried out here, but an even better result is expected since the onset of convective instability in the BL is even further delayed compared with the E1 case (see figure 6 b).