Transition in an infinite swept-wing boundary layer subject to surface roughness and free-stream turbulence

Abstract The instability of an incompressible boundary-layer flow over an infinite swept wing in the presence of disc-type roughness elements and free-stream turbulence (FST) has been investigated by means of direct numerical simulations. Our study corresponds to the experiments by Örlü et al. (Tech. Rep., KTH Royal Institute of Technology, 2021, http://urn.kb.se/resolve?urn=urn:nbn:se:kth:diva-291874). Here, different dimensions of the roughness elements and levels of FST have been considered. The aim of the present work is to investigate the experimentally observed sensitivity of the transition to the FST intensity. In the absence of FST, flow behind the roughness elements with a height above a certain value immediately undergoes transition to turbulence. Impulse–response analyses of the steady flow have been performed to identify the mechanism behind the observed flow instability. For subcritical roughness, the generated wave packet experiences a weak transient growth behind the roughness and then its amplitude decays as it is advected out of the computational domain. In the supercritical case, in which the flow transitions to turbulence, flow as expected exhibits an absolute instability. The presence of FST is found to have a significant impact on the transition behind the roughness, in particular in the case of a subcritical roughness height. For a height corresponding to a roughness Reynolds number $Re_{hh}=461$, in the absence of FST the flow reaches a steady laminar state, while a very low FST intensity of $Tu =0.03\,\%$ causes the appearance of turbulence spots in the wake of the roughness. These randomly generated spots are advected out of the computational domain. For a higher FST level of $Tu=0.3\,\%$, a turbulent wake is clearly visible behind the element, similar to that for the globally unstable case. The presented results confirm the experimental observations and explain the mechanisms behind the observed laminar–turbulent transition and its sensitivity to FST.


Introduction
Aircraft wings are not smooth surfaces but are subjected to various surface imperfections that could cause a premature transition of the laminar boundary-layer flow to a turbulent one and thereby increasing the friction drag. Investigation of effects of the surface imperfections in the case of a swept wing is of particular importance since most current commercial aircraft are built with such a configuration.
The boundary layer over a swept wing exhibits a cross-flow velocity component which is inviscidly unstable due to the presence of an inflection point in that component. The instability modes are vortices having their axes almost parallel to the outer streamlines. If the level of free-stream disturbance is sufficiently low, the dominating vortices are stationary, while in the case of a higher free-stream turbulence (FST) level non-stationary vortices are promoted (for a review on the subject, see Saric, Reed & White (2003)).
Most of the previous works on the discrete roughness elements in three-dimensional boundary layers were conducted for roughness heights smaller than approximately 30 % of the local displacement thickness. Those studies aimed at understanding the receptivity mechanisms to disturbances which promote unstable steady cross-flow vortices (see the review by Kurz & Kloker (2014)). Moreover, utilization of shallow roughness elements as a transition control strategy has been proposed and experimentally confirmed by Saric, Carrillo & Reibert (1998). This strategy has been validated through direct numerical simulation (DNS) by, for example, Wassermann & Kloker (2002) and Hosseini et al. (2013).
The effect of an isolated three-dimensional roughness on transition in two-dimensional boundary layers has been studied since the 1950s. In the work of Gregory & Walker (1955), the authors described the formation of horseshoe vortices in front of the roughness element. The number and characteristics of the vortices where found by Baker (1979) to be mostly dependent on the aspect ratio of the roughness element. From these vortices, further downstream, two streamline-oriented counter-rotating vortices are generated which by lift-up effect (Landahl 1980) create alternated high-and low-speed streaks. The transient growth associated with these streaks has been investigated by various authors such as Choudhari & Fisher (2005), Ergin & White (2006), Denissen & White (2008, 2009 and Cherubini et al. (2013). This growth was found to scale with the square of the roughness Reynolds number, Re hh = U(h)h/ν, where U(h) is the streamwise velocity of the unperturbed flow measured at a wall-normal distance equal to the roughness element height h.
The importance of the parameter Re hh has already been identified in earlier works (e.g. Gregory & Walker 1955;von Doenhoff & Braslow 1961). In particular, von Doenhoff & Braslow (1961) collected existing results and plotted Re hh for which transition was observed (Re hh,cr ) as a function of the aspect ratio η = d/h, where d is the diameter of the roughness element. As is possible to observe in figure 1, when the aspect ratio is increased transition is triggered for lower value of Re hh . It can also be noticed that for a fixed η the spread in the value of Re hh,cr is significant. Loftin (1946) Gregory & Walker (1950) Schwartzberg & Braslow (1952) Klebanoff, Schubauer & Tidstrom (1955) Smith & Clutter (1957) Carmichael ( Doenhoff & Braslow (1961) re-plotted by Kurz & Kloker (2016) together with roughness parameters of experiments by Örlü et al. (2021). Red and green dots correspond to cases where flow in the experiments remained laminar or transitioned to turbulent, respectively, at the lowest FST level. The circles represent cases simulated here (red for the case which transitions and green for those which stay laminar for Tu = 0.0 %). Ergin & White (2006) found that the location of the maximum of the unsteady fluctuations, in the wall-normal and spanwise directions, corresponds to that of the inflection point of the streamwise velocity profile. This suggested that a Kelvin-Helmholtz instability can be the cause of those fluctuations. The measured frequencies were consistent with the previous work by Klebanoff, Cleveland & Tidstrom (1992).
More recently, the onset of laminar-turbulent transition in a flow passing a roughness element has been addressed using global stability analysis (Loiseau et al. 2014;Citro et al. 2015;Kurz & Kloker 2016;Brynjell-Rahkola et al. 2017). Loiseau et al. (2014) performed DNS of the experiments by Fransson et al. (2005) and made an extensive parametric study. They found two different instability modes depending on the aspect ratio of the roughness element: a sinuous one for η < 2 and a varicose one for η ≥ 2. Moreover, they also showed the role of hairpin vortices in triggering of the transition, as shown in previous works (e.g. Acarlar & Smith 1987;Rizzetta & Visbal 2007;Zhou, Wang & Fan 2010;Cherubini et al. 2013). Bucci et al. (2018) reported the appearance of unsteady structures behind a roughness element in the parameter regime that linear stability theory predicts a stable flow. Through numerical simulations and performing optimal forcing analysis, those authors showed that these structures could be triggered by the quasi-resonance of the least stable varicose mode. Moreover, Bucci (2017) demonstrated that the presence of a small amount of FST could trigger the generation of hairpin vortices, similarly to the optimally forced case.
The works on roughness elements with a height of the order of the local displacement thickness immersed in a three-dimensional boundary layer are limited. Kurz & Kloker (2016) studied the transition triggered by the roughness elements on a swept wing for a compressible flow. Those authors compared the vortex system created in threeand two-dimensional boundary layers highlighting effects of the cross-flow velocity component which is responsible for the loss of symmetry. The instability mechanism in the recirculation or in the near-wake area was found to be either absolute or convective, depending on the height of the roughness element. The value of Re hh,cr appeared  Örlü et al. (2021). The numbers inside the cells represent the diameter of the roughness element (mm). Green cases are those with no observable effect, lighter blue those with slight effect, orange those with stronger effect and red those where transition was observed.
not to be significantly influenced by the compressibility or the three-dimensionality of the flow. Brynjell-Rahkola et al. (2017) considered roughness elements placed in a cross-flow-dominated Falkner-Skan-Cooke boundary layer. Those authors found the flow instability shifted from convective to global when the relative height of the roughness was greater than a certain value (h/δ * ≈ 1.3, Re hh ≈ 458, where δ * is the displacement thickness at the roughness location). Moreover, they showed that the numerical solution is very sensitive to the numerical details and in particular to the domain size as the simulated cross-flow vortices were continuously growing. Örlü, Tillmark & Alfredsson (2021) performed an experimental campaign where different roughness dimensions and FST levels were investigated. For each case, based on infrared camera and hot-wire signals, the authors evaluated if transition was taking place behind the roughness. The results are reported in table 1 and compared with data from Gregory &Walker (1955) andvon Doenhoff &Braslow (1961) in figure 1. Increasing the height of the roughness or the the value of the incoming velocity, and thus Re hh , has a destabilizing effect on the wake behind the roughness. The scenario is more complex if we consider effects of the diameter or the level of FST. If we focus on cases with wind tunnel inlet velocity U MTL = 9 m s −1 and h = 0.8 mm, increasing the FST intensity from Tu = 0.03 % to Tu = 0.3 % has a stabilizing effect for d = 32 mm, while for d = 8 mm and d = 16 mm the trend is opposite. If the FST level is kept constant and the diameter is increased we can observe both stabilizing and destabilizing effects. Note that U MTL is different from the free-stream velocity close to the wing due to the installation of contoured side walls.
In the present work, we investigate the configuration and some flow cases from the experiments by Örlü et al. (2021) through DNS. The main objective of our work is to understand the interaction of FST with roughness elements and its effect on the characteristics and nature of the flow instability, as well as explaining some of the observed trends in the experiments. In both the experiments and the cases addressed here, roughness elements with high aspect ratios are considered for which numerical analyses are scarce.

A24-4
Transition in an infinite swept-wing boundary layer

Flow configuration and simulated cases
The flow configuration follows that of the experiment by Örlü et al. (2021) where a wing model with a constant cross-section (modified NACA 67 1 -215) was used. The model was installed at a sweep angle Φ = 35 • and an angle of attack α = −5 • (measured in a plane normal to the leading edge). On the upper surface, at a distance from the leading edge corresponding to 15 % of the chord length (close to the first branch of the neutral stable curve), disc-type roughness elements were placed and their effect on the stability of the boundary layer was studied. Different upstream velocity, roughness dimensions and FST intensity were considered. In figure 2, the airfoil shape and different coordinate systems used here are shown. A Cartesian coordinate system (x, y, z) is used where the x axis is oriented perpendicular to the leading edge of the wing, the z axis parallel to the leading edge and the y axis normal to these two. Some results are presented in a curvilinear reference system (τ, n, z), where τ is the chordwise direction tangent to the wing surface and n normal to the surface of the wing. The corresponding velocity components are (u τ , v n , w). For visualization purposes, a third reference system obtained by rotating the Cartesian one 30 • around the y axis, denoted as (x , y, z ), is also employed.
For the cases studied here, the Reynolds number is Re = U ∞ c/ν = 581 100. Here, U ∞ = 11.1 m s −1 is the total free-stream velocity, ν is the kinematic viscosity and c = 800 mm is the cord length measured normal to the leading edge. This value of the Reynolds number corresponds to the wind tunnel inlet velocity U MTL = 9 m s −1 in table 1.
The main geometrical parameters of the roughness elements studied here are summarized in table 2. With the exception of case D, the sizes of the roughness elements are the same as in the experiment. In particular, case A corresponds to h = 0.4 mm and d = 16 mm, and cases B, C and E have a height h = 0.8 mm and diameters d = 8, 16 and 32 mm, respectively. By choosing these roughness sizes the effects of the height and diameter of the roughness on the flow can be studied independently. These cases have been selected based on the observed flow behaviour in the experiments. Cases B, C and E all have a roughness Reynolds number Re hh = 461 which is close to the critical value reported by Brynjell-Rahkola et al. (2017) and below the critical value reported by Kurz & Kloker (2016). The flow in the experiments shows different behaviour for different aspect ratios d/h, while Re hh is constant. Also a small amount of FST seems to have a strong effect on transition. Therefore, first, simulations without FST were performed to find the threshold  Table 2. Geometrical and flow parameters measured at x = 0.15 for different roughness sizes. Displacement thickness at this location is δ * = 6.96 × 10 −4 . Velocities u τ,max and u τ,min are the maximum and minimum chordwise perturbation velocities measured in the n-z plane at x = 0.2.
for the occurrence of the global instability. In the experiments for case C, depending on the level of FST intensity, either a clear transition or some visible effects on the heat transfer were observed. This makes it a suitable case to explore the effects of FST intensity on instability of the wake behind the roughness. Finally, some cases were studied from a flow stability point of view by means of an impulse-response analysis.

Numerical approach
We assume the flow evolution is described by the incompressible Navier-Stokes equations where u = {u, v, w} T is the velocity vector, p the pressure and f is the body force.
Here, U ∞ is used as reference velocity and c as reference length. In order to integrate (3.1) the spectral element method proposed by Patera (1984) is used. In particular, we employed the code Nek5000 of Fischer, Lottes & Kerkemeier (2008). In this formulation, the equations are solved in the weak form and the numerical domain is partitioned using hexahedral elements. Within each element, Lagrange interpolants on N p + 1 Gauss-Lobatto-Legendre nodes are used to discretize the velocity while N p − 1 Gauss-Legendre nodes are used for the pressure in order to avoid spurious modes. In the present work, we have used N p = 11 for most of the simulations. The equations are advanced in time treating the nonlinear terms explicitly using a third-order extrapolation scheme (EXT3) whereas the viscous terms are discretized with a third-order backward differentiation scheme (BDF3).

Computational domain and boundary conditions
Different computational domains used in the current work are presented in figure 2(b). The solid line denotes the domain used for the laminar base flow simulations and the dashed line denotes the boundary for the simulations including the roughness elements. The spanwise extension of the domain is L z = 0.07. This width is sufficient to ensure that the stability of the recirculation bubble behind the roughness element is not affected by the employed periodic boundary condition in the spanwise direction (von Doenhoff & Braslow 1961). The relation between L z and the spanwise wavelength of the instability modes is discussed in § 3.4. The whole mesh consists of 14 750 spectral elements for the base-flow computations and 216 144 elements for cases including the roughness elements. Its structure at z = 0 is shown in figure 3(a). Details of the mesh around a roughness element are given in figure 3(b). The boundary conditions for DNS are found as follows. First, a Reynolds-averaged Navier-Stokes (RANS) simulation of the experimental set-up (with the wing model placed in the wind tunnel including the test section walls) is performed assuming homogeneous flow along the wing span. In these calculations, transition was prescribed at x ≈ 0.7 on the upper surface and at x ≈ 0.2 on the lower one. These locations correspond to those in the experiments where the boundary layer was tripped. Then, a DNS is performed to obtain the laminar base flow in the domain corresponding to the solid line in figure 2(b). Data from this simulation are then used to set the boundary conditions in the main simulations.
Here, Dirichlet conditions are prescribed for velocities at the inflow boundary. In DNS of the laminar base flow, the following values have been used: At the outflow, a modified version of the natural boundary condition is used: 1 Re where the value of p a is computed from the RANS solution: These boundary conditions are similar to those employed by Shahriari, Kollert & Hanifi (2018) where a similar multiple outflow configuration was used. On the upper and lower boundaries a Dirichlet condition is used for the chordwise and spanwise directions whereas a modified natural condition is prescribed for the cross-stream component v: Further, on the lateral boundaries a periodic boundary condition is set. In all present simulations a sponge region is prescribed at the outflow to avoid numerical instabilities. In this region, the flow is forced to the steady state of the smooth-wing simulations. This is done by adding a forcing on the right-hand side of the momentum equations that reads Here, U f (x) is the target velocity field taken from the RANS simulations, λ f a smooth step function and A f the amplitude of the forcing, chosen to be A f = 100.

Free-stream turbulence
In order to numerically generate the FST, we superimpose a number of Fourier modes to the base flow at the inlet. In order to have a homogeneous and isotropic turbulence, we divide the wavenumber space into 80 concentric shells where each shell represents the amplitude associated with a given wavenumber. Then, 40 points are randomly chosen on each shell giving the components of the wavenumber vectors. Finally, amplitudes of these modes are chosen to match the von Kármán spectrum: (3.7) Here, E is the energy density, k the wavenumber and Λ the integral length scale. This procedure is described in detail in Schlatter (2001) and has been used in simulations of pitching wings (Negi et al. 2018) and wind turbines (Kleusberg 2019). Similar to the experiments, the FST level is evaluated in the free stream above the roughness elements.
Here, we have considered FST intensities Tu = 0.03 % and Tu = 0.3 % which correspond to the natural level of FST in the wind tunnel and the lowest intensity of the grid-generated turbulence studied in the experiments, respectively. The integral length scale was not measured in the experiments. We have used a value of Λ = 2.5 × 10 −3 to make sure that the relevant part of the von Kármán spectrum, including scales with the maximum energy, can be resolved in our simulations. The amplitude of the spectrum at the inflow boundary is chosen such that the desired FST intensity is achieved at the location of the roughness elements. The quality of the developed FST is checked through examination of its decay rate. In figure 4 the FST intensity as a function of chordwise coordinate for the case of Tu = 0.3 % is given. There, a power-law fitting Tu ∝ (x − x 0 ) n with n = −0.83 is also included, validating correct decay rate of the simulated FST. the experiment (Örlü et al. 2021):

Validation of laminar base flow
(3.8) The pressure at x = 0.43 is used as p ref in (3.8), and U ∞ and p ∞ are the values of the incoming velocity and pressure, respectively. As can be seen in figure 5(a), a good agreement between the data is found confirming correct treatment of the boundary conditions. Small differences observed here can be attributed to the presence of the traverse system in the experiments.
In order to check the accuracy of the computed boundary-layer flow, DNS data are compared with the profiles obtained with a boundary-layer code using the pressure distribution from the DNS. In figure 5(b), the chordwise velocity u τ profiles inside the boundary layer are shown. As can be seen, a close agreement is found, indicating the boundary layer in the DNS being well resolved.
Further, grid convergence studies have been performed by varying the polynomial order. Results of these studies are presented in Appendix A.
To further characterize the flow in the absence of roughness elements, the wall-normal maximum value of velocity in the direction perpendicular to the local inviscid streamline, the cross-flow component, is shown in figure 5(c). This quantity reaches its maximum value of 6.9 % for x = 0.041 and then decreases to 5.0 % towards the end of the numerical domain.
3.4. Stability characteristics of the undisturbed flow A nonlocal stability analysis of the mean flow has been performed using the NoLoT PSE code Hein et al. 1994). A summary of the results is given in figure 6. As can be observed there, for the conditions of the experiments, amplification of the steady and unsteady cross-flow modes is low. The steady cross-flow vortices have their maximum growth (N ≈ 0.8) for β ≈ 200 and the unsteady perturbations reach an amplification of N ≈ 2 for β ≈ 350 close to the end of the computational domain at x = 0.6. This indicates that in the absence of the roughness elements these perturbations will not cause transition over the smooth surface. This has actually been the reason behind the choice of flow parameters in the experiments. The disturbance amplifications reported here are in agreement with those measured experimentally by Borodulin et al. (2019) for the same geometry and flow conditions. The spanwise extension of the computational domain, L z = 0.07, corresponds to a spanwise wavenumber of β ≈ 90. The isolated roughness elements used here will also trigger other spanwise wavenumbers. However, none of these would reach high amplitudes due to the linear growth of the cross-flow modes, as suggested by the results presented in figure 6. With the resolution of our simulations, the unstable portion of the spectrum can be computed with a high accuracy.

Results
We start our investigation by comparing the flow behaviour for different sizes of the roughness elements in the absence and presence of FST. Further, the instability mechanisms behind the observed flow behaviour are studied through impulse-response analysis.

Flow behaviour in the absence of FST
First, we examine the critical roughness size for a laminar inflow (Tu = 0). A summary of the obtained results is presented in figure 7 where the flow fields from the analysed cases are shown. In all cases, except case D, a steady laminar solution is reached and the effect of the roughness is limited to the creation of streaky flow patterns which are described in the following sections. The results indicate that for cases with E) the flow remains laminar. In case D, transition is triggered due to the global instability since numerical noise is thought to be extremely low (this is investigated further through impulse-response analysis). Our observations indicate that, for the current flow case, the critical roughness Reynolds number Re hh,cr is between 461 and 712 for all aspect ratios studied. Concentrating on the flow around the roughness elements, we can identify two recirculation zones upstream and downstream of them, similar to those found by Kurz & Kloker (2016). In figure 8(b), the isosurface of the velocity u τ = 0 shows a weak asymmetry due to the presence of the cross-flow component which is about 6 % at the location of the roughness (figure 5c). In figure 8(a), extensions of the recirculation zones obtained with polynomial order N p = 7 and N p = 11 are compared showing a good agreement, confirming grid convergence.

Effects of roughness element height
In this section, we analyse the effects of roughness height by comparing the steady solutions for cases A, C and D. In these cases, the roughness diameter is d = 0.02 (16 mm) and the aspect ratios are d/h = 40, 20 and 14.39, respectively, corresponding to h/δ * = 0.72, 1.44 and 2.0. We visualize the structures created by the roughness elements by subtracting the reference velocity field (smooth geometry) from the steady solutions. As mentioned above, in case D the flow does not reach a steady state and laminar-turbulent transition takes place close to the roughness. In this case, a steady base flow is obtained using the BoostConv algorithm proposed by Citro et al. (2017). For the sake of completeness, a description of the method is given in Appendix B.
The structures of the generated vortices behind the roughness elements are presented in figure 9. The first observation is that these are not symmetric as in the case of the Blasius boundary layer (see Loiseau et al. 2014). This is due to the existence of the cross-flow velocity component. In the studies by Kurz & Kloker (2016) and Brynjell-Rahkola et al. (2017) the generated high-speed streaks behind the roughness elements merged to a single one. But, in the cases studied here, we did not observe the coalescence of the high-speed streaks. This can be due to the fact that the generated vortices in the present work are weaker and present limited growth of their amplitude (figure 10). Moreover, their lateral distance (due to the high aspect ratio) is larger than those in the previous works.
The strengths of the observed streaks are measured at x = 0.2 and reported in table 2. These are defined as the maximum (u τ,max ) and the minimum (u τ,min ) of the chordwise perturbation velocity u τ in the planes normal to the surface. As expected, the strength of the vortices increases with increasing roughness height. However, on increasing the height of the roughness, the velocity deficit grows much faster than the amplitude of the high-speed streaks. As is evident from table 2, the amplitude of the central low-speed streaks, in the range of our study, seems to increase with the square of roughness height. Comparing cases D and C, low-speed streak amplitude increases by a factor of 2 while the ratio of the roughness heights is 1.4. Similar behaviour is observed when comparing cases A and C: the ratio of the low-speed streaks is about 4.6 while the ratio of the roughness element heights is 2. Case A has the lowest roughness height which results in the lowest streak amplitude as reported in table 2. In figure 9(a) (left-hand column), we can see how the streaks develop in the numerical domain. They grow in amplitude with increasing x with the exception of the leftmost high-speed streak whose amplitude starts to decay close to the end of the domain.
In figure 9(b), for case C, we see that at x ≈ 0.2 one major low-speed streak develops behind the roughness element. At its side, two high-speed and two weaker external low-speed streaks develop, too. Initially, the two high-speed streaks next to the central low-speed one have a similar amplitude but the rightmost one becomes dominant while developing downstream. Both the central and the leftmost low-speed streaks become weaker further downstream while the rightmost one increases its amplitude. As can be seen in figure 9(c), the structure of disturbances in case D is very similar to that of case C but with higher streak amplitudes. This is also confirmed in figure 10 where the maximum and minimum values of u τ behind the roughness are reported for cases C and D. In particular, the strength of the high-speed streak does not change significantly for x > 0.2 and in both cases reaches a maximum at x ≈ 0.3.
As mentioned above, case D is the only one in which flow turns to turbulent in the absence of FST. In order to see the differences between the cases studied here, we choose to look at the spanwise shear which was found to be more sensitive to variation of roughness size. Figure 11 shows the spanwise derivative of the tangential velocity (∂u τ /∂z) at x = 0.2 for selected cases. As can be seen there, spanwise shear in case A is much weaker than in other cases. In case C, the maximum value of spanwise shear is located between the central low-speed streak and the right-hand high-speed one. In case D, the negative spanwise shear is the dominating one. Further, in this case a double-vortex structure generating higher shear is observed. This concentrated high spanwise shear gives the major contribution to the growth of the instabilities as indicated by the energy budget reported in § 4.2.2.

Effects of roughness element diameter
In this section we compare cases B, C and E, where the roughness height is kept constant and the diameter is varied between d = 0.01 and 0.04. In all these cases a steady flow was found for Tu = 0.0 %. As is evident from table 2, the values of u τ,max and u τ,min are weakly affected by the variations of the roughness diameter which can also be seen in figure 12. We can see that increasing the diameter results in generation of wider streaks, in particular the central low-speed streak becomes larger. Even here, none of the streaks merge as in contrast to what was observed by Kurz & Kloker (2016) and Brynjell-Rahkola et al. (2017). Further, the central low-speed streak in case E grows all the time while in cases B and C it starts to decay after an initial phase of growth.
In figure 13, values of ∂u τ /∂z for cases B, C and E are plotted. As can be seen there, the maximum and minimum values of the spanwise shear do not differ between those cases. However, increasing the roughness diameter results in slightly lower spanwise shear. This is due to the fact that the streak amplitudes are mainly affected by the roughness height and not by its diameter, and that the streaks get wider when increasing the roughness diameter.

Unsteady solution and effects of FST
As mentioned above, in case D the flow does not reach a steady state and transition to turbulence takes place in the vicinity of the roughness element. Similar flow behaviour is observed for case C in the presence of FST. In this section the mechanisms behind this transition are investigated. Figure 14(a) shows the isosurface of spanwise velocity for case D. The transition process starts with disturbances being generated in the low-speed region behind the roughness which soon become chaotic. This is similar to what was found in simulations by Kurz & Kloker (2016) for the case of a compressible swept-wing boundary layer. For configuration of case C, once a FST of Tu = 0.3 % is imposed at the inflow, the transition is triggered in the wake of the roughness similar to what is observed in case D. The flow structures are plotted in figure 14(b). It can be noticed that the wake of the roughness is a highly receptive flow region as mostly this area is affected by FST. This observation is consistent with the findings of Bucci et al. (2018) for a two-dimensional boundary-layer flow. The transition scenario in case C is significantly modified when a lower value of FST intensity (Tu = 0.03 %) is selected. This value corresponds to the natural level of FST in the MTL wind tunnel in which the experiments of Örlü et al. (2021) were performed. In figure 15, the isosurface of the spanwise velocity is reported for different time instants. As seen there, there is no fixed transition point behind the roughness, but some turbulent spots appear randomly. In the first snapshot, we can observe a large turbulent spot close to the outflow boundary and a smaller one in middle of the domain. In the other snapshots we can follow the evolution of the latter spot in time.

Flow structures
While the developed turbulent spots leave the domain, new ones are generated which in turn propagate downstream. So, although such a low level of FST does not trigger the transition directly after the roughness, it can change the flow behaviour significantly. In figure 16, the instantaneous friction coefficient C f is plotted for case C with Tu = 0.0 % and 0.3 % as well as case D with Tu = 0.0 %. In the first case, we can recognize the footprint of the high-and low-speed streaks demonstrated as the region of the higher and lower C f , respectively. In case C with Tu = 0.3 % and case D with Tu = 0.0 %, first traces of some spanwise vortices are observed. Once transition takes place, a region of elongated chordwise structures with mainly higher values of C f appears and the turbulent region spreads in the spanwise direction as it propagates downstream.
To further investigate the observed unsteady flows, we subtract the steady solution from the instantaneous flow fields and visualize the structures. The isosurfaces of negative and positive chordwise velocity u τ are shown in figure 17. Again, this indicates that initially  velocity at the location of the roughness height, are f h = f * h * /Q * (h * ) = 0.076 and 0.113, respectively. If instead the diameter of the roughness element is used as the reference length, the non-dimensional frequencies are f d = f * d * /Q * (h * ) = 1.5 and 1.6, respectively. Further, case C presents relatively high fluctuations also for 10 < f < 30 which are not present in case D. This may be due to the noise generated by FST. In figure 19, we can see that also the r.m.s. fields present similarities. In both cases the high r.m.s. values are located close to z = 0 where also the highest values of ∂u τ /∂z were found (see figure 11). At this spanwise position, the mean value of the chordwise velocity isū τ ≈ 0.7.

Impulse-response analysis
In order to further study the stability characteristics of the flow behind the roughness elements, we performed an impulse-response analysis similar to the one carried out by Brandt et al. (2003), Peplinski, Schlatter & Henningson (2015) and Brynjell-Rahkola et al. (2017). We introduced a wave-packet disturbance, as in Bech, Henningson & Henkes (1998), upstream of the roughness element and monitored its evolution. This is done in a linear framework where the steady solutions discussed in previous sections are used as the base flow. In particular, we analyse evolution of the disturbance kinetic energy to obtain some insights about the linear behaviour of small perturbations. In figure 20, the evolution of the total kinetic energy K T (integrated in the whole domain) of the wave packet as a function of time is given for cases B, C, D and E. As can be seen there, in all cases K T decreases initially and then it starts growing for t ≈ 0.1. With the  exception of case D, K T decreases again for t > 0.2. The evolution of the wave packet for case D is very different. While the energy of the wave packet in cases B, C and E exponentially decays after it has reached its maximum, in case D it grows exponentially continuously after an initial transient stage. This indicates that flow in case D may be globally unstable.
In order to further analyse the response of the flow to the impulse, we monitor the disturbance kinetic energy of the wave packet K (integrated in the n-z plane) as a function of space and time. Results for cases C and D are presented in figure 21. The same analysis for cases B and E showed a similar behaviour of K. Figure 21(a) shows the evolution of K as a function of chordwise coordinate x for case C. The disturbance energy grows  As shown in figure 21(b), the disturbance energy, which is significantly larger than in case C, grows contentiously in time.
In the case where the flow is globally unstable, the wave packet will travel in both upstream and downstream directions, while if it is convectively unstable, disturbances will propagate downstream and eventually leave the computational domain. In the latter scenario, the velocity of the trailing edge of the wave packet will be positive, whereas in the case of global instability it will be negative (Schmid & Henningson 2001). In order to find the propagation velocity of the wave packet, we plot K as a function of c g = (τ − τ 0 )/t, where τ 0 is the tangential coordinate of a point just behind the roughness element. Figure 21(c) shows the disturbance energy K as a function of c g for different time instants. We can observe that the rear part of the curves for t > 0.15 crosses at a negative value of c g , which indicates that the wave packet travels upstream. Hence, this case appears to be globally unstable. This result is consistent with that of the nonlinear simulations where flow transitioned to turbulent in the absence of external disturbances (Tu = 0.0 %).
In figure 22, the wave packets for cases C and D at t = 0.25 are given showing the isosurfaces of negative and positive chordwise velocity. As seen there, in both cases the wave is located in the low-speed region behind the roughness element. In case C, the wave packet evolves forming elongated low-and high-speed regions inclined to the direction of the vortices created by the roughness element. In case D, the structure is almost perpendicular to the axis of the vortex with a shorter wavelength. One can see a similarity between these structures and those found in our nonlinear simulations illustrated in figure 17. These structures also bear similarities to the global modes found in the works of Kurz & Kloker (2016)  . The one observed in case C corresponds to the low-to moderate-frequency modes (of the z-mode type ) in those works, while the structure in case D is similar to the high-frequency global modes. The similarity to results of global mode analysis in the mentioned works is also confirmed by the results presented in figure 18, where the dominant frequencies in the flow are shown.

Energy budget analysis
Here, we analyse the energy budget of the wave packet discussed above by computing the contributions from the production and dissipation terms of the Reynolds-Orr equations (Schmid & Henningson 2001): where the dissipation term D is given by with u being the perturbation velocity. The production terms are the following: 3) Here, {u τ , u η , u z } are the three velocity components of the stationary solution. In the case where the sum of the production terms is greater than the dissipation, the energy of the perturbation grows in time. The production terms can be either positive or negative causing destabilization or stabilization of the perturbations, respectively. In figure 23(a) we report the ratio between the total production and the dissipation for cases C and D. As seen there, in case D the ratio is always greater than one, apart from at the first instants, and thus we have a positive contribution to the growth of the perturbation. On the contrary, in case C we can see that the dissipation terms are greater than the sum of the production ones for t > 0.22 and so the kinetic energy of the perturbation decreases. Up to that time, also in case C we can observe an initial growth of the perturbation which is in agreement with the study of a two-dimensional boundary layer performed by (Cherubini et al. 2013). These observations are consistent with the behaviour of K T shown in figure 20. For case D, the contribution from each term of (4.1), scaled with the value of the dissipation, is shown in figure 23. It is found that the term Pr 3 gives the largest contribution. This is in agreement with the observation of Loiseau et al. (2014) for the unstable varicose mode generated by a roughness element in a two-dimensional boundary layer. The main difference between their results and those presented in figure 23 is the significant contribution from the term Pr 9 . This is due to the fact that here the spanwise velocity component of the base flow has values comparable to those of the chordwise one.

Summary and conclusions
The flow over a swept wing in the presence of disc-type roughness elements has been investigated by means of DNS. The main goal of this study is to investigate the nature of the instability caused by the roughness elements with high aspect of ratio and their interaction with FST. Through our numerical simulations, we try to explain the behaviour of wake flow behind the roughness elements (the heat transfer investigated using infrared images) and its sensitivity to FST intensity observed in the experiments of Örlü et al. (2021). To do so, FST with the same intensity and integral length scale as in the experiments is imposed at the inflow.
Nonlinear simulations performed with zero FST showed that for roughness elements with different widths and a height up to 0.8 mm (h/δ * = 1.44, Re hh = 461), flow behind the roughness was steady and no transition was observed. On the contrary, simulations for h/δ * = 2.0 (Re hh = 712) showed a turbulent flow just behind the roughness element. In the experiments, with a FST level of Tu = 0.3 % transition was observed at the location of the roughness for h/δ * = 1.44 and d/h = 20. Even at the lowest FST level (Tu = 0.03 %), a strong effect on transition was observed for this case. These observations were confirmed by our simulations when the same levels of FST were used. The common feature of the flow behind the roughness elements, when non-transitional, is the existence of a series of high-and low-speed streaks formed due to the lift-up effect. Due to the presence of the cross-flow component of the base flow, the spanwise distribution of these streaks is not symmetric as in the case of a two-dimensional boundary layer. The rightmost high-speed streak, as well as the low-speed streaks next to it, grows in amplitude along the streamwise direction. The central low-speed streak on the contrary decays. Different from the observation of Kurz & Kloker (2016) and Brynjell-Rahkola et al. (2017), the vortices generated behind the roughness elements did not merge here. This can be due to the wider roughness elements used in our studies. Moreover, we have observed that increasing the diameter does not cause a significant change in the stability of the flow.
The transition caused by large roughness elements is thought to be a result of the global instability of the flow generated by these elements (Loiseau et al. 2014;Kurz & Kloker 2016;Brynjell-Rahkola et al. 2017). Here, we have investigated the instability characteristics of the flow behind the roughness elements by means of an impulse-response analysis. For all cases with h/δ * = 1.44, the wave packet generated by the impulse grows weakly downstream of the roughness until x ≈ 0.25 and then decays and exits the computational domain. The scenario is very different for the case with h/δ * = 2.0, where the wave packet grows by several orders of magnitude and expands in the space behind the roughness. Further, the trailing edge of the generated wave packet has a negative propagation velocity meaning the flow is absolutely unstable in this case. The origin of the flow instability is found to be located in the low-speed flow region behind the roughness and creates hairpin-vortex-like structures which then rapidly break down. An energy budget analysis of case D showed, as expected, that the total contribution of the production terms is greater than the dissipation resulting in growth of the total kinetic energy. The terms associated with the spanwise gradient of the base flow are the ones giving the major contributions. On the contrary, in case C after the initial transient the dissipation is higher than production and so the perturbation is damped.
Although the impulse-response analysis for case C revealed a weak growth of the wave packet, a continuous forcing with a FST level of Tu = 0.3 % was found to be sufficient to trigger the flow transition behind the roughness element. This transition takes place in the low-speed region as in the unstable case, but the flow pattern appears to be more chaotic. A comparison between the time signals and their r.m.s. values retrieved in the transitional region showed that the fluctuations for case C with Tu = 0.3 % and case D have a similar spectral content and spatial support. This suggests that in both cases those fluctuations may have been governed by similar instability modes. Outside the wake the flow does not appear to be significantly affected which indicates that the low-speed region is highly receptive to the forcing caused by the FST. We also have found that for Tu = 0.03 % flow exhibits travelling turbulent spots which randomly appear due to the convectively unstable disturbances. These observations were found to be consistent with experimental findings (Örlü et al. 2021) where for h/δ * = 1.44 and d/h = 20 with Tu = 0.3 % a clear transition was observed, while for Tu = 0.03 % a strong effect was found but not a clear transition. Our analysis confirmed a large impact of FST level on the flow behaviour behind the subcritical roughness elements. A similar conclusion was drawn in Bucci et al. (2018) and Bucci (2017) for the case of a two-dimensional boundary layer. briefly describe the basic idea of the algorithm. For further details, readers are referred to Citro et al. (2017) who also applied the method to stabilize the flow behind a hemispherical roughness element in a two-dimensional boundary layer. Let us consider the linear system If we apply an iterative method to compute the solution, we can write where r n = b − Ax is the residual at step n and B the matrix that describes the iterative procedure. Then, (B2) can be rewritten as The convergence of the iterative procedure depends on the spectrum of the matrix (I − BA). Often, the convergence is dictated only by a small part of the spectrum, either the least damped or the most amplified eigenvalues, depending on the system being stable or unstable, respectively. The algorithm aims to modify this part of the spectrum. Now, we consider the evolution equation for the residual: which can be obtained from (B2) applying −A and adding b on both sides. The central idea is to replace the residual r n with a new one ξ n (r n ). We need to choose ξ n such that ξ n → 0 only when r n → 0. If we substitute r n , equation (B4) reads r n+1 = r n − ABξ n (r n ).
Now, the task is to to choose ξ n such that r n+1 is minimized. To do so, a least-squares method is adopted to approximate the solution of ABξ n = r n . The value of ξ n is approximated using the subspace of AB built through its repeated action. It should be noticed that the described procedure requires only the knowledge of the residuals r n to compute ξ n and for that reason can be included as a black box in the original procedure.