Global stability analysis and direct numerical simulation of boundary layers with an isolated roughness element

Global stability analysis and direct numerical simulation (DNS) are performed to study boundary layer flows with an isolated roughness element. Wall-attached cuboids with aspect ratios $\eta=1$ and $\eta=0.5$ are investigated for fixed ratio of roughness height to displacement boundary layer thickness $h/\delta^*=2.86$. Global stability analysis is able to capture the frequency of the primary vortical structures. For $\eta=1$, only varicose instability is seen. For the thinner roughness element ($\eta=0.5$), the varicose instability dominates the sinuous instability, and the sinuous instability becomes more pronounced as $Re_h$ increases, due to increased spanwise shear in the near-wake region. The unstable modes mainly extract energy from the central streak, although the lateral streaks also contribute. The DNS results show that different instability features lead to different behavior and development of vortical structures in the nonlinear transition process. For $\eta=1$, the varicose mode is associated with the shedding of hairpin vortices. As $Re_h$ increases, the breakdown of hairpin vortices occurs closer to the roughness and sinuous breakdown behavior promoting transition to turbulence is seen in the farther wake. A fully-developed turbulent flow is established in both the inner and outer layers farther downstream when $Re_h$ is sufficiently high. For $\eta=0.5$, the sinuous wiggling of hairpin vortices is prominent at higher $Re_h$, leading to stronger interactions in the near wake, as a result of combined varicose and sinuous instabilities. A sinuous mode captured by dynamic mode decomposition (DMD) analysis, and associated with the `wiggling' of streaks persists far downstream.


Introduction
Surface roughness has important effects on boundary layers. The laminar-turbulent transition in boundary layer flows can be greatly modified by the presence of localized or distributed roughness. Understanding the roughness-induced transition process is therefore important for the control of flows in an engineering context.
The roughness Reynolds number is one of the important parameters in roughnessinduced transition, and empirical correlations have been developed to predict the transition onset. One commonly used definition is Re h = U e h/ν, where U e is the boundary layer edge velocity, h is the roughness height and ν is the kinetic viscosity of the fluid. This definition, however, does not account for the impact of the relative location of the roughness element in a boundary layer. Another definition Re hh = u h h/ν based on the Blasius velocity solution at the roughness tip location u h has been suggested to best characterize roughness-induced transition (Klanfer & Owen 1953). Von Doenhoff & Braslow (1961) have suggested a transition diagram that correlates the roughness aspect ratio η = d/h (where d is the roughness width and h is the roughness height) with Re hh , using experimental dataset from different types of distributed surface roughness. Their diagram highlights the crucial roles of η and Re hh on the transition onset.
Isolated, three-dimensional (3-D) roughness elements may be considered as the primary models to be generalized and extended for more complex roughness geometries. The effects of isolated roughness on transition have been investigated experimentally by Gregory & Walker (1956). The main flow pattern is observed to be horseshoe vortices that wrap around the roughness element, and whose legs trail downstream and give birth to the streamwise vortices farther downstream. Baker (1979) has experimentally studied the vortex system around an isolated cylindrical roughness element, and shown the dependency of the horseshoe system dynamics on Re D = U e D/ν and D/δ * , where D is the cylinder diameter and δ * is the displacement boundary layer thickness. The streamwise vortices induced by the 3-D roughness elements create longitudinal streaks downstream that are lifted upwards (Landahl 1980;Reshotko 2001). These streamwise longitudinal streaks are related to the disturbance transient growth, which can be unstable and cause transition downstream of the roughness (Fransson et al. 2004(Fransson et al. , 2005. The concept of optimal perturbation was introduced by Böberg & Brösa (1988) and Butler & Farrell (1992) to define these 'most dangerous' initial perturbations that generate the maximum energy growth. Luchini (2000) provided a numerical method to determine the optimal perturbation and explain that the linear growth of initially small disturbances can excite nonlinear interactions and cause transition.
Both symmetric (termed 'varicose') and anti-symmetric (termed 'sinuous') streak instabilities have been detected and are of importance in turbulent boundary layers. The varicose type is associated with horseshoe vortices that originate from a normal inflectional instability in the streamwise velocity profile (Robinson 1991;Asai et al. 2002;Skote et al. 2002). The sinuous streak instability is correlated with a base state with a spanwise inflection and contributes to the regeneration of near-wall turbulence (Jiménez & Moin 1991;Skote et al. 2002). Asai et al. (2002) observed that wider streaks more easily undergo varicose breakdown while narrower streaks are more likely to undergo the sinuous breakdown.
The strength and stability properties of the streamwise streaks also play important roles in roughness-induced transition, and are dependent on roughness characteristics, such as its shape, height and aspect ratio. White et al. (2005) conducted experiments to investigate the effects of surface roughness characteristics on transient growth and suggested a strong dependence of the resulting transient growth on roughness diameters. Piot et al. (2008) used bi-global stability approach and DNS to investigate the local stability of streamwise streaks past a single smooth roughness element. A stabilizing effect of a 'pre-streaky' structure associated with the counter-rotating streamwise vortices is identified in the near wake. Cherubini et al. (2013) performed a global optimization analysis to search for the optimal perturbation inducing the largest transient growth for a boundary layer past smooth 3-D roughness elements. They investigated bumps with different heights at different target times, and their results indicate that when the bump is sufficiently large, it can destabilize the wake flow on a short time scale.
With large-scale linear algebra computations now being possible, global linear stability theory (Theofilis 2011) has been performed on roughness-induced transition. Global stability is especially useful for non-parallel flows such as roughness wakes, thus is a promising tool to predict and analyze roughness-induced transition. Loiseau et al. (2014) used global stability theory to investigate the flow past a cylindrical roughness element. They suggested that the frequencies associated with the dominant fluid dynamics are well predicted by global stability analyses, and that the unstable nature of the central low-speed streak is of crucial importance in the transition process. Citro et al. (2015) presented the direct and adjoint global eigenmodes for boundary layer flows past a hemispherical roughness element, and found that the critical Reynolds number is constant when the ratio of the roughness height and the displacement boundary layer thickness h/δ * is less than 1.5. Kurz & Kloker (2016) used DNS and global stability analysis to investigate the effects of discrete surface roughness with various roughness height and background disturbance on a swept-wing boundary layer. Their results suggest that large elements are able to trip turbulence by either a convective or a global instability in the near-wake region. Puckert & Rist (2018) conducted experiments to detect global instability in the 3-D isolated cylindrical roughness cases of Loiseau et al. (2014). They report that the critical Reynolds number is higher than the transition Reynolds number and suggest that global instability may not be the decisive mechanism for transition.
While linear instability is usually detected in the near wake of roughness elements using local or global stability analysis, secondary instability could appear farther downstream if the streak amplitude is sufficiently large. Andersson et al. (2001) found that secondary instabilities appear at large amplitudes of the primary streaks and suggested that the sinuous modes of instability are dominant and most often reported for streak breakdown, while the varicose instability only occurs for larger streak amplitude and thus is barely observed in natural transition. Denissen & White (2013) studied the stability of steady roughness-induced transient growth to unsteady fluctuations, and found that the transition could be caused via secondary instability in the mid-wake region when the roughness size is sufficiently large. Vadlamani et al. (2018) found that secondary sinuous instability developed on the low-speed streak is evident in the transition process induced by distributed roughness. They also suggested that this sinuous like breakdown is reminiscent with the sinuous instability observed in the context of transition under the effects of free-stream turbulence Hack & Zaki 2014).
Understanding global varicose and sinuous instabilities in roughness-induced transition and their dependency on roughness configuration is important. De  conducted bi-global and three-dimensional parabolized stability analyses to investigate the transition induced by a sharp-edged isolated roughness element in a supersonic boundary layer. Their results suggest that the varicose mode is associated with the entire 3-D shear layer while the sinuous mode is a consequence of the lateral streaks. Loiseau et al. (2014) suggest that the sinuous global mode is similar to the von Kármán vortex street in the 2-D cylinder flow and the varicose mode is associated with the hairpin vortices. They also investigate the dependence of instability types on aspect ratios and suggest that varicose instability is observed for wider roughness elements and sinuous instability is observed for thinner roughness elements. Puckert & Rist (2018) reported experimental observation of sinuous oscillations in roughness-induced transition using PIV. They found that for thin roughness elements (i.e., η = 1), the sinuous mode competes with the varicose mode and becomes dominant in the supercritical regime. Bucci et al. (2021) fixed the aspect ratio of the cylinder to unity and studied the effect of freestream turbulence on roughness-induced transition. They noted that the roughness Reynolds number and aspect ratio might not be the only important parameters for flow characteristics, the shear ratio also plays a crucial role in the onset and symmetry of the primary global instability.
The joint effects of the parameters mentioned above make the instability characteristics and transition process highly sensitive to the flow configuration. Although the ratio h/δ * seems to play a crucial role in the determination of the dominant instability, the dependency and sensitivity of the onset of sinuous instability on η and Re h need further investigations to complement the current understanding. While it has been observed that the sinuous instability is related to a sinuous wiggling similar to von Kármán vortices, the influence of sinuous instability and its interplay with the varicose instability on the nonlinear patterns and dynamics needs further analysis. To address these points, we perform global stability and adjoint sensitivity analyses to study global instability of boundary layer flows over a cuboid with two aspect ratios. The ratio of the cuboid height to the displacement boundary layer thickness is 2.86, which is larger than most past work. The thin cuboids with aspect ratios η = 1 and 0.5 are considered to examine the dominant instability and the sensitivity of the onset of sinuous instability in terms of η and Re h for roughness elements with high h/δ * . We also perform DNS to examine the dependence of Re h and η on the laminar-turbulence transition process, and use DMD analysis to study the development of vortical structures and associated non-linear dynamics corresponding to different global instability characteristics.
The numerical methodology is introduced in §2 and validations of the global stability and adjoint sensitivity solver are shown in §3. The flow configuration, base flow computation, grid convergence and domain length sensitivity are demonstrated in §4. The results and discussions are presented in §5. Finally, the paper is summarized in §6.

Numerical methodology
The governing equations and numerical method are briefly summarized. An overview of modal linear stability, adjoint sensitivity and details regarding the iterative eigenvalue solver are provided.

Direct numerical simulation
The incompressible Navier-Stokes (N-S) equations are solved using the finite volume algorithm developed by Mahesh et al. (2004): where u i and x i are the i-th component of the velocity and position vectors respectively, p denotes pressure divided by density, ν is the kinematic viscosity of the fluid and K i is a constant pressure gradient (divided by density). Note that the density is absorbed in the pressure and K i . The algorithm is robust and emphasizes discrete kinetic energy conservation in the inviscid limit which enables it to simulate high-Re flows without adding numerical dissipation. A predictor-corrector methodology is used where the velocities are first predicted using the momentum equation, and then corrected using the pressure gradient obtained from the Poisson equation yielded by the continuity equation. The Poisson equation is solved using a multigrid pre-conditioned conjugate gradient method (CGM) using the Trilinos libraries (Sandia National Labs). The DNS solver has been validated for a variety of problems on wall-bounded flows, including: realistically rough superhydrophobic surfaces (Alamé & Mahesh 2019), random rough surfaces (Ma et al. 2021) and response of a plate in turbulent channel flow (Anantharamu & Mahesh 2021).

Linear stability analysis
Linear stability analysis enables the investigation of the linearized dynamics of infinitesimal perturbations evolving on a base state. In the present work, the incompressible Navier-Stokes equations are linearized about a base state, u i and p. The flow can be decomposed into a base state subject to a small O( ) perturbationũ i andp. The linearized Navier-Stokes (LNS) equations are obtained by subtracting the base state from equation (2.1) and can be written as follows: The same numerical schemes as the N-S equations are used to solve the LNS equations. The LNS equations can be rewritten as a system of linear equations, where A is the LNS operator andũ i is the velocity perturbation. The solutions to the linear system of equations (2.3) are: Re(ω) is defined as the growth rate and Im(ω) is the temporal frequency. The linear system of equations can then be transformed into a linear eigenvalue problem: where ω j = diag(Ω) j is the j-th eigenvalue andû j i = U i [j, :] is the j-th eigenvector. For the global stability analysis, the computational cost to solve the eigenvalue problem using direct methods is very expensive. Instead, a matrix-free method, the implicitly restarted Arnoldi method (IRAM) is usually used. We make use of the IRAM implemented in the PARPACK library to solve for the leading eigenvalues and eigenmodes.

Adjoint sensitivity analysis
Adjoint sensitivity analysis solves for the dominant eigenvalues and eigenmodes of the adjoint LNS Equations, which yields the dominant sensitivity modes corresponding to the direct modes. According to the definition of the continuous adjoint to the LNS equations by Hill (1995), the adjoint equations are: Similar to the direct problem, the adjoint systems of linear equations can be simplified to an eigenvalue problem: (2.7) Hill (1995) suggested that the adjoint perturbation velocity field would highlight the optimal locations where the largest response to unsteady point forcing occurs. In the present work, our aim is to use the global adjoint sensitivity analysis in conjunction with the global stability analysis to determine the most sensitive flow regions to point forcing and the inception of instability.

Validation
The global stability solver is developed in the present work on 3-D structured platforms. First, the global stability of a 3-D lid-driven cavity is validated against Regan & Mahesh (2017). Then, the global stability and adjoint sensitivity analyses are performed for laminar channel flow, where the results are compared to the parallel flow stability of Poiseuille flow conducted by Juniper et al. (2014). The global adjoint sensitivity results are also examined.

3-D lid-driven cavity
3-D lid-driven cavity is studied as the first validation case for the global stability solver. The Reynolds number based on the cavity height and the lid velocity is 1000.
The leading eigenvalues solved from the global stability solver show good agreement with Regan & Mahesh (2017) in table 1. The isocontours of the real part of the leading eigenmodes in figure 1 show good qualitative agreement with the results in Gomez et al. (2014) and Regan & Mahesh (2017). Regan & Mahesh (2017) Present Table 1. The leading eigenvalues from global stability analysis for a stable 3-D lid-driven cavity at Re = 1000 compared to Regan & Mahesh (2017).

Laminar channel flow
Global stability and adjoint sensitivity analyses of a laminar channel are performed. The global stability results are compared to results from parallel flow stability analysis. The Reynolds number Re τ = 44.7 is based on the friction velocity u τ and the channel half height δ ch . The domain length is 4πδ ch in the streamwise direction and (4π/3)δ ch in the spanwise direction. Since the streamwise and spanwise wavenumbers α and β are not specified in global stability analysis, any combination of those can be present in the global stability results. Thus we can extractû j i for a selective combination of α j and β j using streamwise and spanwise Fast Fourier transforms. The selected combination of α j and β j can be used as the input into the parallel flow stability. The leading eigenvalues from global stability and adjoint sensitivity analyses show agreement with those from the parallel flow linear stability analysis in table 2. The non-zero components of the first leading direct and adjoint eigenmodes are shown in figure 2. For the contour plots, qualitative agreement is shown compared to the results from Regan & Mahesh (2017. A quantitative comparison between the global stability and parallel flow linear stability results is shown in figure 2(c). Good agreement is obtained for |ŵ|. The results of the second leading direct and adjoint eignmodes are shown in figures 3 and 4 respectively. Both qualitative and quantitative agreement are obtained.

Problem formulation
In this section, the simulation set-up is shown, the base flow computation is described, and a study of grid convergence and domain length sensitivity is performed.

Flow configuration
The flow configuration, the computational domain and the roughness geometries are depicted in figure 5. At the inflow, a laminar Blasius boundary layer profile is prescribed. The cuboid with the height h and width d is centered at the origin of the Cartesian coordinate system. The ratio of the roughness height to the displacement thickness of (a)   the boundary layer h/δ * is fixed at 2.86. Two aspect ratios η = d/h = 1 and 0.5 are investigated. The roughness height is h = 1, the reference length in the simulations. The streamwise extent of the computational box L x is 45h for global stability analyses, and is extended in the DNS to examine the transition process farther downstream. The spanwise extent is L z = 10h to ensure that the roughness element behaves as isolated, and the wall-normal extent is L y = 15h. The distance from the inlet of the computational domain to the center of the roughness element is denoted by l = 15h. The Blasius laminar boundary layer solution is specified at the inflow boundary, and convective boundary conditions are used at the outflow boundary. Periodic boundary conditions are used in

Base flow computation
Linear stability analysis requires a stationary base flow. The time-invariant state can be either the equilibrium, or the time-averaged (mean) flow. For flows at moderate Reynolds numbers, the equilibrium state can be obtained using the selective frequency damping (SFD) method (Åkervik et al. 2006) or the BoostConv algorithm (Citro et al. 2017). For turbulent flows at higher Reynolds numbers, the equilibrium state is difficult to obtain; instead, the time-averaged mean flow can be used as the base state for stability analysis to seek meaningful physical interpretation (Turton et al. 2015;Tammisola & Juniper 2016). Barkley (2006) shows that linear stability analysis on cylinder wake flow using the mean flow as the base state is able to track the Strouhal number of vortex shedding, but yields a marginally stable state with a small growth rate. In the present work, we use SFD to compute the base flow, compare this base flow to the time-averaged mean flow, and compare their global stability results in §5.
SFD introduced byÅkervik et al. (2006) is a useful technique to artificially settle the flow towards a steady equilibrium. The main idea is to apply a temporal low-pass filter (a) x/h = 0 to damp the oscillations due to the unsteady part of the solutions, and is achieved by introducing a linear forcing term on the right-hand side of the Navier-Stokes equations. An encapsulated formulation of the SFD method developed by Jordi et al. (2014) is used in the present work. The problem is considered to have converged when ||q−q|| inf 10 −8 according to Jordi et al. (2014), where q is the filtered state. When using SFD, the control coefficient χ and the filter width ∆ play important roles in the convergence process. The control coefficient χ should be positive and larger than the growth rate of the desired mode, while the filter cut-off frequency ω c = 1/∆ must be lower than all of the flow instabilities to ensure the unstable disturbances are well within the damped region. For example, χ = 0.5 and ∆ = 2 are used for the unstable case (Re h , η)=(600, 1), and the convergence history is shown in figure 6.

Grid convergence and domain length sensitivity
Global stability results show strong sensitivity to grid sizes and domain lengths, highlighted by Loiseau et al. (2014) for roughness wake flow and Peplinski et al. (2015) for a jet in crossflow. A study on grid convergence and domain length sensitivity is thus performed in this section. Three different grids are used in the grid convergence study which are referred to as 'coarse', 'medium' and 'fine'. Simulation details are listed in table 3. For all cases presented in table 3, uniform grids are used in both streamwise and spanwise directions while non-uniform grids are used in the wall-normal direction. Compared to the coarse grid, the medium grid is refined in the wall-normal direction. In the finest grid, the grid spacing in all three directions is reduced. Table 3 presents ∆y spacing at the wall (denoted by ∆y wall ) and ∆y spacing at the roughness height location (denoted by ∆y top ). Note that the roughness element is resolved by 43, 86 and 172 grid points in the wall-normal direction for the coarse, medium and fine cases respectively.
The streamwise velocity profiles of the base flow are examined at three different stations Table 3. Simulation parameters for grid convergence and domain length sensitivity study, and comparison of the direct leading eigenvalue for Case (Re h , η)=(600, 1). in figure 7. The results show significant deviation of the solution for the coarse grid, while the differences between the medium and fine grids are small, indicating grid convergence. The leading eigenvalues obtained from the global stability analysis also show convergence in table 3, suggesting that the medium grid is adequate for global stability analyses on the present case. The influence of streamwise domain length on the global stability results is examined in the simulation with L x = 75h (denoted by Case L x 75). Simulation details are listed in table 3. Note that the grid sizes in Case L x 75 are comparable to the medium grid, already proven to sufficiently resolve the flow. The leading eigenvalue shows good agreement with that of Case Medium in table 3, suggesting that the streamwise domain length L x = 45h is adequate for the present case. The leading eigenmodes in Case Medium and Case L x 75 are also depicted in figure 8. The results are identical between both cases. The global mode decays appreciably before reaching the outflow boundary, which guarantees convergence in the global stability results. Based on these conclusions, the medium grid and the domain length L x = 45h are used for the cases presented in §5.1 and §5.2. Figure 9. Contour plots at the spanwise mid-plane of (a) time-averaged streamwise velocity field for Case (Re h , η)=(600, 1), streamwise velocity field of the base flow obtained from SFD for (b) Case (Re h , η)=(600, 1), (c) Case (Re h , η)=(475, 1), (d) Case (Re h , η)=(600, 0.5) and (e) Case (Re h , η)=(800, 0.5).

Base flow
The differences between the base flow (obtained from the SFD method) and the time-averaged mean flow (obtained from DNS) are first investigated. Figures 9(a) and 9(b) show that the streamwise velocity field upstream of the cube is identical between the base flow and the mean flow, but the flow field downstream of the cube shows some differences. Although the reversed flow regions behind the cube are similar, with increasing downstream distance, the base flow demonstrates a stronger wall-normal gradient corresponding to the shear layer generated from the top edge of the cube. In contrast, this strong wall-normal shear is not prominent in the mean flow. Figure 10 compares the base and time-averaged mean flows for Case (Re h , η)=(600, 1), using streamlines and contours of streamwise velocity in different cross-flow planes. Qualitative agreement is seen between the base flow and the mean flow immediately downstream of the roughness element (x 4h). At x = 0, as shown in figure 10(a), two pairs of streamwise vortices are observed on the lateral sides of the cube in both base and mean flows. The pair very close to the cube is referred to as the symmetry plane vortices (SP) (Iyer & Mahesh 2013) or the rear pair vortices (Ye et al. 2016;Bucci et al. 2021). They push low-momentum flow upwards, move closer to the symmetry plane, give rise to the low-speed region behind the roughness, and are dissipated farther downstream. The other counter-rotating vortex pair is formed away from the symmetry plane, which is referred to as the off-symmetry plane vortices (OSP) by Iyer & Mahesh (2013). They are the continuation of the vortex tubes from the horseshoe vortex system upstream. At x = 4h ( figure 10(b)), hairpin vortices (H) and secondary wall-attached vortices (SW) are observed in both the base and mean flows. As the streamwise location increases farther downstream (figures 10(c) and 10(d)), the central low-speed region is weakened and secondary vortical structures are intensified in the mean flow, due to the enhanced nonlinear interactions. In the base flow however, only one pair of wall-attached vortices is observed below the primary hairpin vortices since the unsteady oscillations are damped out. The influence of the non-linear saturation in the mean flow results in differences in the global stability results between the base and mean flows, as further discussed in §5.3.2.
The dependence of the base flow features on different η and Re h is examined in figures 9(b)-9(e). The spanwise vortices observed upstream of the roughness element correspond to the horseshoe vortex system induced by the stagnation effect of the roughness. Baker (1979) suggested that the stability and topology of the horseshoe vortex system is mostly dependent on Re h and h/δ * . For η = 1, the location of the horseshoe vortex moves slightly farther from the front face of the roughness as Re h increases, shown in figures 9(b) and 9(c), consistent with the observations by Daniel et al. (2017). Also, the shear layer induced by the roughness lifts up and shows a stronger wall-normal gradient as Re h increases. For η = 0.5, shown in figures 9(d) and 9(e), the regions corresponding to the upstream spanwise vortices and the dowmstream reversed flow are smaller due to thinner roughness geometry. The Re h dependence for η = 0.5 is similar to what is observed for η = 1.
The high-and low-speed streaks are examined in figure 11, using isosurfaces of the streamwise velocity deviation u d = u − u bl . For Case (Re h , η)=(600, 1), the central lowspeed streak and two lateral low-speed streaks are illustrated in figure 11(a). The central low-speed streak, which occurs symmetrically with respect to the mid-plane, originates from the flow separation downstream of the roughness element. The lateral low-speed streaks are associated with the counter-rotating vortices. High-speed streaks close to the wall appear farther downstream. Figure 11(b) shows that for Case (Re h , η)=(600, 0.5), the thinner roughness geometry leads to thinner and less sustainable central and lateral low-speed streaks, and the high-speed streaks are absent farther downstream. For Case (Re h , η)=(800, 0.5), figure 11(c) shows that the strength of the central and lateral lowspeed streaks gets amplified as Re h increases. In contrast to the other two cases, the high-speed streaks are prominent in the near-wake regions, indicating increased spanwise shear that would contribute to the sinuous instability examined in §5.2. Combining the above results and the smaller h/δ * results of Loiseau et al. (2014), it can be concluded that: first, larger h/δ * , larger η and higher Re h lead to a stronger wall-normal shear and a more sustainable central low-speed streak; Second, increasing Re h for thin roughness could result in an increased spanwise shear in the near wake region.

Global stability analysis
Global stability analysis has been performed for cases with η = 1 and η = 0.5 at different Re h , and the leading eigenvalues are shown in figures 12(a) and 12(b) respectively. For η = 1, one leading eigenvalue is obtained at each Re h , as shown in figure 12(a). The case at Re h = 450 is absolutely stable, consistent with the steady flow field observed from the DNS results. As Re h increases, both the growth rate and the temporal frequency are increased. The critical Re h can be identified when the growth rate of an eigenvalue becomes positive. The flow at Re h = 475 is marginally stable which suggests that the critical Re h is close to 475 for this configuration.
For η = 1, the eigenmodes of the leading eigenvalues are all varicose for the various Re h investigated. The real part of the leading eigenmodes is shown for Re h = 475 and Re h = 600 in figures 13(a) and 13(b). Both the leading stable and unstable global modes exhibit a varicose symmetry with respect to the spanwise mid-plane. As shown by the 3-D view of the eigenmode, the shape and location of the modes are consistent with those of the central low-speed streak observed in figure 11(a). The varicose mode demonstrates the unstable nature of the central low-speed region induced by the roughness element. Compared to the stable mode at Re h = 475, the unstable mode at Re h = 600 is more notably lifted, corresponding to the more raised shear layer for higher Re h observed in figure 9. For η = 0.5, a different unstable behavior is shown in figure 12(b). One leading stable eigenvalue is seen at Re h = 450 and its associated mode is varicose. Two leading eigenvalues are obtained at higher Re h . The eigenvalue with larger growth rate and lower frequency is a varicose mode, and the other eigenvalue with smaller growth rate and higher frequency is a sinuous mode. For the thinner roughness geometry, the sinuous instability becomes more prominent as Re h increases. The associated varicose and sinuous eigenmodes of the leading eigenvalues for Case (Re h , η) = (800, 0.5) are visualized in figures 13(c) and 13(d). While the varicose mode is associated with the central low-speed streak observed in figure 11(c), the sinuous mode shows a larger streamwise extent along the central region. These results indicate that both varicose and sinuous oscillations exist in the wake flow, and the effect of sinuous instability could be more persistent on the The contour levels are shown within the range from −1.0e −7 (blue) to 1.0e −7 (red). The localized shear is depicted by the solid lines of us = ((∂u/∂y) 2 + (∂u/∂z) 2 ) 1/2 from 0 to 2. The orange dashed lines show the location of the element.
transition process. It thus can be concluded that for thin roughness with large h/δ * , while the varicose instability is dominant, the sinuous instability can also be present. The onset of sinuous instability results from the interplay of small η and increased Re h , corresponding to the enhanced spanwise shear observed in the near wake of the base flow with increasing Re h .

Production of disturbance kinetic energy
The production of disturbance kinetic energy provides insight into how and where the global modes extract their energy from the base flow. As illustrated by De  and Loiseau et al. (2014), the main contributions to the production of disturbance kinetic energy are the two terms The streamwise variation and spatial distribution of these two dominant terms are examined for Cases (Re h , η)=(600, 1) and (Re h , η)=(800, 0.5). The spatial variations of P y and P z in cross-flow planes are depicted in figure 14. In The lateral low-speed streaks also make a contribution to the dominant production terms when h/δ * is large. The mode extracts energy from the lateral streaks, as shown at x = 10h in figure 14(b). The top views of P y and P z for Case (Re h , η)=(600, 1) demonstrated in figure 15(a) display the contributions of the two lateral streaks more clearly. The large shear ratio h/δ * leads to stronger central and lateral streaks in the present case. Although the varicose mode extracts most of energy from the central lowspeed streak, the contribution of the lateral streaks can not be neglected for cases with large shear ratios.
In contrast, the contours of P y and P z for the leading varicose and sinuous modes of Case (Re h , η)=(800, 0.5) are shown in figures 14(c) and 14(d) respectively. The distribution of P y demonstrates that while the varicose mode extracts the energy from the top edge of the central streak, the sinuous mode extracts its energy from the lateral parts of the central streak. These results are consistent with the observation by Loiseau et al. (2014) for small h/δ * cylindrical roughness. For the thinner geometry (η = 0.5), there is less fluid passing above the roughness element, and a stronger spanwise shear is seen corresponding to the longer wall-normal extent for the lateral parts of the central streak, shown in figure 14(d). This suggests that the sinuous instability occurs due to the fact that it could extract more energy from the spanwise shear. The contour plots of P y and P z at y = 0.75h are shown in figures 15(b) and 15(c). The P y and P z distributions of the sinuous mode show a longer streamwise extent than those of the varicose mode, implying the influence of sinuous instability on the wake flow could last farther downstream. Both the varicose and sinuous modes are able to extract some energy from the lateral streaks. The contribution of the lateral streaks is associated with the strength of the lateral streaks which is more likely dependent on the shear ratio.

Adjoint sensitivity analysis
The adjoint perturbation velocity field highlights the most receptive regions to momentum forcing. The leading adjoint eigenvalues are computed for Cases (Re h , η)=(600, 1) and (Re h , η)=(800, 0.5). The results show good agreement with their associated direct eigenmode counterpart in table 4. The streamwise velocity component of the leading adjoint modes is depicted in figure 16. The adjoint modes are located immediately upstream of the roughness element as well as on the top edge of the separation region directly above and downstream of the roughness element. The receptive region for Case  (Re h , η)=(800, 0.5) is smaller than that for Case (Re h , η)=(600, 1) due to the thinner geometry. The adjoint mode is symmetric with respect to the spanwise mid-plane, corresponding to the direct varicose mode, and is anti-symmetric corresponding to the direct sinuous mode. Due to the large differences between the spatial distribution of direct and adjoint modes, neither direct nor adjoint solution alone can describe the whole picture. The product for each jth pair of direct and adjoint global modes computed as W j (x, y, z) = ||û j ||||û †,j || max(||û j ||||û †,j ||) , determines the region where the eigenvalues of the LNS operator are most sensitive to localized feedback (Giannetti & Luchini 2007), -also called the "wavemaker" regions. Locations where W ≈ 1 are sensitive to localized feedback, corresponding to the instability core. The value of W can be interpreted as quantification of a possible change in the eigenvalues as a result of applied forcing in the given region of the flow (Ilak et al. 2012). Figure 17 depicts the wavemaker regions for the leading modes of Cases (Re h , η)=(600, 1) and (Re h , η)=(800, 0.5). As shown in figure 17(a), the wavemaker for the varicose mode in Case (Re h , η)=(600, 1) is prominent on the top edge of the reversed flow region and over an extended region along the central low-speed streak. The maximum value of the wavemaker at each z − y plane is plotted along the streamwise direction in figure 18. It is shown that the wavemaker has its maximum value within the separation region corresponding to the instability core, and drops to the order of 10 −1 as it passes through the reversed flow region. Similar features are also seen for the varicose mode of Case (Re h , η)=(800, 0.5), as depicted in figure 17(b). However, for the thinner roughness geometry, the spatial growth shows a much smaller streamwise extent and a sharper decline is observed in figure 18. This could be related with the weaker and shorter central streak for the thinner roughness element. For the sinuous mode in Case (Re h , η)=(800, 0.5), figure 17(c) shows that the instability core is also located in the downstream reversed flow region, but mainly associated with the lateral sides of the reversed flow region. Unlike the varicose mode that shows one primary sensitivity region, the sinuous mode shows two lateral sensitivity regions. The wavemaker maximum of the sinuous mode shows an immediate drop behind the reversed flow region but a gradual decrease compared to that of the varicose mode. It thus can be concluded that a thinner roughness geometry results in a weaker spatial extension of the wavemaker, and the wavemaker strength of the varicose mode diminishes more quickly than that of the sinuous mode.

Nonlinear breakdown
To understand the effects of η and Re h on transition, and the role of different instability characteristics in the non-linear evolution downstream of the roughness element, DNS are performed for cases with η = 1 and η = 0.5 at different Re h . Longer streamwise domain lengths are used to examine the transition process.

Joint effect of η and Re hh on transition
Dependence of the transition process on η and Re hh is examined for the present cases by reproducing the von Doenhoff-Braslow transition diagram (Von Doenhoff & Braslow 1961) in figure 19. This transition diagram shows the correlation between η and Re hh and provides an approach to predict the transition characteristics. The cases located in Region (i) (below the lower curve) are expected to have a steady wake flow, the cases fitting into Region (ii) (between the lower and upper curves) indicate that the wake flow is unsteady and transition occurs, and the cases situated in Region (iii) (above the upper curve) means that transition occurs immediately downstream of the roughness. The correlation between η and Re hh revealed by figure 19 shows that for roughness elements with 0.6 η 2, η plays a more important role on the onset of unsteadiness than the onset of immediate transition downstream of roughness, while the opposite effect is seen for roughness with 0.3 η 0.6.
As shown in figure 19, for η = 1, the stable cases at Re h = 450 and 475 are located in Region (i), consistent with a stable wake flow observed in both the global stability and DNS results. The unstable cases at Re h = 600 and 800 are within Region (ii) and the transition process is examined in figures 20(a) and 20(b) respectively. For Case (Re h , η)=(600, 1), both the near and farther wakes are symmetric with respect to the spanwise mid-plane, corresponding to the varicose mode obtained in the global stability analysis.
Compared to Case (Re h , η)=(600, 1), symmetric fluid motions with smaller length scales are seen in the near-wake region for Case (Re h , η)=(800, 1), indicating that nonlinear breakdown occurs more closely downstream of the roughness as Re h increases. The perturbations farther downstream are more prominent than those at Re h = 600. Intense shear between streaks results in spanwise oscillations in the farther wake. Although only the varicose global instability is detected for this configuration, a sinuous like breakdown could happen when Re h is sufficiently high. This sinuous breakdown might be related or subsequently lead to secondary sinuous instabilities observed by Denissen & White (2013) and Vadlamani et al. (2018), which would destabilize the shear layer and promote transition to turbulence. Note that whether or not the unstable cases undergo transition to turbulence is not revealed in this transition diagram since the effects on other configuration parameters, such as spanwise spacing and Re δ , need to be considered. Case (Re h , η)=(1100, 1) is located in Region (iii), suggesting that immediate transition downstream of the roughness is expected. This is verified in figure 20(c), as Re h increases to 1100, the non-linear breakdown occurs immediately downstream of the roughness. Also, streamwise streaks can be identified farther downstream, implying that transition to fully-developed turbulence might have occurred.
For η = 0.5, the stable case at Re h = 450 and the unstable case at Re h = 800 are situated in Regions (i) and (ii) respectively, as shown in figure 19, which is consistent with the global stability results. The unstable case at Re h = 600 is slightly off from Region (ii), which could be due to the fact that the present roughness has sharp edges and would have a lower Re hh for unsteadiness to occur than other smoother roughness elements.
The effect of small η and Re h dependence on the transition process is examined for η = 0.5. Figure 20(d) shows that the wake flow at Re h = 600 displays a thinner symmetric central streak compared to that of η = 1. There are no sinuous oscillations observed, corresponding to the global stability results that the sinuous mode is marginally stable at Re h = 600. As Re h increases to 800 (figure 20(e)), the anti-symmetric oscillations in the spanwise direction become evident in both the near and farther wakes, associated with the more prominent sinuous instability and indicating a persistent effect of sinuous oscillations farther downstream. In summary, cases with two thin roughness geometries in the present work are well fitted into the classification by this transition diagram, and the interplay between η and Re hh leads to different flow behavior in the transition process, in accordance with the global instability characteristics.

Non-linear evolution of vortical structures
Cases (Re h , η)=(600, 1) and (Re h , η)=(800, 0.5) are examined to better understand how non-linear saturations are triggered and how vortical structures develop following different types of global instability. Figure 21 shows the vortical structures using isocontours of Q (Chong et al. 1990). For Case (Re h , η) = (600, 1), as shown in figure 21(a), the vortical motions induced by side edges of the cube interact with the shear layer generated over the cube, giving birth to the hairpin vortices. Both the primary hairpin vortices and secondary wall-attached vortices are observed downstream of the roughness element. These vortical structures are amplified and fragment into small structures beyond x = 18h, which is a manifestation of transition. The vortex "head" and "legs" are advected, stretched downstream, and diminish after x = 44h. Even though the unsteadiness is noticeable and the unstable nature of the longitudinal streaks is revealed at Re h = 600, transition to turbulence may not happen since Re δ for the boundary layer is low in the present case.
In contrast, Case (Re h , η)=(800, 0.5) shows different vortical motions in the wake flow in figure 21(b). As the horseshoe vortices wrap around the roughness element and interact with the shear layer, anti-symmetric distribution of vortical structures is seen in the immediate vicinity, downstream of the roughness. This indicates that sinuous oscillations occur just downstream of the roughness element. The primary hairpin vortices modulated by the sinuous oscillations of the central streak also exhibit a sinuous wiggling. While the hairpin vortices are advected and stretched downstream, and break down at about x = 18h, the sinuous wiggling of the streaks continues farther downstream. It is clear that for Case (Re h , η)=(800, 0.5), both the varicose and sinuous instabilities have influences on the behavior and development of vortical structures. The varicose instability is associated with the hairpin vortices, thus has a limited streamwise extent as the hairpin vortices break down. The sinuous instability originates from the immediate vicinity of the roughness, is correlated with the wiggling of the central streak, and has a more persistent effect than the varicose instability on the wake flow.
The time history of streamwise velocity probed at three stations is examined for Case (Re h , η)=(600, 1) in figure 22(a). Periodic oscillations with a circular frequency ω = 1.088 are seen at different streamwise stations, corresponding to the periodic shedding of hairpin vortices. These self-sustained oscillations independent of external noise, are also a sign of global instability (Puckert & Rist 2018), and their frequency is close to the temporal frequency of the leading global varicose mode. For Case (Re h , η)=(800, 0.5), stronger fluctuations with multiple frequencies are observed in the immediate vicinity of the roughness element, and smaller amplitude of velocity variations is seen at the farther stations.
To understand the dynamics behind this different behavior for Cases (Re h , η)=(600, 1) and (Re h , η)=(800, 0.5), dynamic mode decomposition (DMD) was performed and compared to the energy spectra and global stability results. DMD is a data-driven modal decomposition technique that identifies a set of modes from multiple snapshots of the observable vectors. Each of the DMD modes has an assigned eigenvalue that describes its temporal growth/decay rate and oscillation frequency. DMD is a useful tool to isolate the regions associated with a particular frequency and provide information on system dynamics. For the present work, we use a novel DMD algorithm developed by Anantharamu & Mahesh (2019) that is suitable for analysis of large datasets. The basic idea behind DMD is that the set of snapshot vectors of flow variables {ψ i } N −1 i=1 can be written as a linear combination of DMD modes {φ i } N −1 i=1 as where λ j are the eigenvalues of the projected linear mapping and c j are the jth entries of the first vector ψ 1 . The detailed derivation of the algorithm can be obtained from Anantharamu & Mahesh (2019). To ensure the accuracy of the results, N=200 snapshots of the flow field were taken with ∆tU e /h=0.15 between them for Case (Re h , η) = (600, 1), and N=700 snapshots were taken with ∆tU e /h=0.1 between them for Case (Re h , η) = (800, 0.5). Also, power spectral density (PSD) are examined at the mid-plane for different streamwise stations downstream of the roughness. For Case (Re h , η)=(600, 1), the PSD shows a primary peak at the Strouhal number St = 0.175 in figure 23(a), corresponding to the shedding frequency of the main hairpin vortices and the secondary wall-attached vortices observed in figure 21(a). The interaction between different vortical structures results in the higher harmonics at St = 0.35 and St = 0.525. Note that similar peaks are also identified in the DMD spectra (figure 23(c)). Table 5 demonstrates a comparison of the eigenvalues and Strouhal numbers obtained from global stability and DMD analyses. The Strouhal numbers obtained from global stability analysis and DMD analysis show good agreement. It is worth noting that using the mean (time-averaged) flow as the base state for global stability analysis can still capture the shedding frequency of hairpin vortices, but the mode is marginally stable with a small growth rate. This discrepancy in the growth rates between base flow and  mean flow is similar to the observations by Barkley (2006) for the linear stability analysis of the cylinder wake flow. The state of marginal stability is due to the strong nonlinear saturation of the mean flow observed in figure 10. Sipp & Lebedev (2007) conducted a global weakly nonlinear analysis for cylinder flow and provided theoretical explanation for the marginal stability of mean flows: the zeroth harmonic is much stronger than the second harmonic. This could explain the fact that the mean flow is marginally stable in the present work. The associated global unstable mode of the mean flow and the DMD mode are examined in figure 24. They both demonstrates varicose features and show good qualitative agreement. Compared to Case (Re h , η)=(600, 1), a combination of multiple frequencies is distributed in the energy spectra and the DMD spectra for Case (Re h , η)=(800, 0.5), indicating more complicated flow behavior. Figure 23 The results indicate that the interactions between the hairpin vortices and the general sinuous oscillations are significant in the near wake, and diminish as the vortical structures develop farther downstream. As the streamwise station increases farther downstream, a peak at a low frequency St = 0.120 in figure 23(b) gets amplified. This peak is also captured in the DMD spectra ( figure 23(d)). The corresponding DMD mode in figure  25(a) shows a sinuous symmetry. This sinuous mode is associated with the wiggling of (a) the streaks observed farther downstream in figure 21(b). It is thus clear that different roughness geometries associated with different instability characteristics lead to different wake flow behavior in the transition process.

Mean flow characteristics
The transitional flow behavior is examined using the time-averaged flow. Figure 26 shows the streamwise variation of the time-averaged skin friction at three different stations across the span for cases with η = 1 and η = 0.5 at different Re h . For Case (Re h , η)=(600, 1), shown in figure 26(a), the C f value at z = 0 shows a prominent increase behind the roughness location, corresponding to the progress from the reversed flow region to the downstream region. At z = 0.5h and z = h, the C f profiles show peaks around x = 18h, which is associated with the evolution of the lateral wall-attached lowspeed streaks observed in figures 10 and 11. As the streamwise distance increases farther downstream, the mean skin friction at three stations does not collapse to the same level, indicating the flow experiences unsteadiness but the saturation is not sufficiently strong, and transition to turbulence may not happen eventually.
As Re h increases to 800, the peaks of C f profiles at three stations move closer to the roughness and drop to a similar level farther downstream, shown in figure 26(b). This suggests that as Re h increases, the onset of unsteadiness occurs more closely to the roughness and the wake flow becomes more homogeneous. Figure 26(c) shows that for Re h = 1100, the increase of C f profiles occur in the immediate vicinity behind the roughness element compared to those at a lower Re h . As the streamwise location increases beyond x = 20h, the C f profiles at three stations collapse and remain constant, suggesting that the wake flow becomes homogeneous and transition to turbulence may occur downstream. Note that whether or not transition to turbulence occurs could also depend on Re δ of the boundary layer. For the high shear ratio h/δ * considered in the present work, Re δ corresponding to a certain Re h is relatively low. The transition to turbulence is thus unlikely to happen at moderate Re h .
In contrast to η = 1, the C f profiles for Case (Re h , η)=(800, 0.5) are examined in figure 26(d). Due to a thinner geometry, the C f profile at z = 0 demonstrates a lower level, and its sharp rise occurs more closely to the roughness element compared to Case (Re h , η)=(800, 1). The C f profiles at z = 0.5h and z = h remain at a lower level since the width of the wake flow and the spacing of two lateral streaks are smaller for a thinner roughness geometry. Note that as the streamwise distance increases, the C f profile at z = h increases and the C f values at z = 0.5h and z = h are slightly higher than that at the mid-plane in the late stages of transition. This corresponds to the wiggling streaks observed in figure 21(b), indicating that the effect of the wiggling streaks on the mean flow persists and contributes to the transition process.
The boundary layer evolution from laminar to turbulent states is examined in figure  27(a) using the mean velocity profiles in wall units at different streamwise locations downstream of the roughness element for Case (Re h , η)=(1100, 1). The time-averaged streamwise velocity at the mid-plane is normalized by the local friction velocity u τ , where u τ is computed from the C f profile at z = 0 for each x location. The wall-normal coordinate in wall units is y + = yu τ /ν. The results show that all profiles collapse well in the viscous sublayer and follow the correlation U + = y + . From x = 5h to x = U + y + y + log(y + )/0.40 + 5.25 40h, significant increase is observed above the viscous layer, which is due to the lift-up behavior of the shear layer. The mean velocity profile above the viscous sublayer reaches its maximum magnitude at x = 40h and decreases to approach the log-law profile as the x location increases farther downstream. Agreement with the logarithmic law is seen beyond x = 100h, indicating that the inner layer is fully-developed. As the x location increases even farther, the profiles at x = 110h and x = 130h show agreement in both the inner and outer layers, suggesting that fully-developed turbulent flow is established in both the inner and outer layers. The velocity fluctuations and Reynolds shear stresses at x = 130h are depicted in figure 27(b) using wall scaling. The velocity fluctuations u rms , v rms and w rms are normalized by u τ,ave , and the Reynolds shear stress u v is normalized by u 2 τ,ave , where u τ,ave is the spanwise-averaged friction velocity computed from the spanwise-averaged C f at x = 130h. The results show good agreement with the results of a turbulent zero-pressure gradient boundary layer from Schlatter et al. (2010).

Conclusions
Global stability analysis and direct numerical simulation are performed to study roughness-induced transition. Isolated cuboids with aspect ratios η = 1 and η = 0.5 immersed in laminar boundary layers at different Re h are investigated. The ratio h/δ * between the roughness height and the local displacement boundary layer thickness is 2.86, which is higher than most past studies.
Differences between the base flow computed using the SFD method and the mean flow obtained from DNS are examined. The base flow shows a stronger wall-normal shear farther downstream than the mean flow, while non-linear interactions are more evident in the mean flow. Either using the base flow or the mean flow as the base state for global stability analysis is able to capture the shedding frequency of the primary vortical structures. However, the mean flow evolves to a marginally stable state due to the strong non-linear saturation, in contrast to an unstable base flow. We thus use the base flow as the base state for global stability analysis in the present work.
The effects of η and Re h on the base flow are investigated. As Re h increases, the downstream shear layer lifts up and shows a stronger wall-normal gradient. For a thinner roughness geometry, the central and lateral low-speed streaks are thinner and less sustainable compared to the thicker roughness at the same Re h . It can be summarized that higher Re h , larger η and higher h/δ * lead to a stronger wall-normal shear and a more sustainable central streak. Also, as Re h increases for the thinner roughness, high-speed streaks below the central streak become prominent in the near-wake region, indicating an increased spanwise shear that contributes to sinuous instability.
Global stability analysis shows that when the shear ratio is sufficiently high (h/δ * = 2.86), the varicose instability is dominant for the roughness element with small aspect ratios (η 1). For η = 1, both the stable and unstable modes exhibit varicose symmetry. For η = 0.5, the varicose instability is dominant at different Re h , and the sinuous instability becomes more pronounced as Re h increases. These results highlight how the onset of sinuous instability is highly dependent on the joint effects of h/δ * , η and Re h . At smaller h/δ * , smaller η and higher Re h , the sinuous instability is more likely to occur.
The production of disturbance kinetic energy shows that the varicose mode extracts energy from the wall-normal and spanwise shear of the central streak. In contrast, the sinuous mode extracts its energy from the lateral parts of the central streak. A longer wall-normal extent of the central streak for the thinner roughness geometry leads to a stronger spanwise shear, and the sinuous instability is able to extract more energy from the spanwise gradient of the base flow and becomes prominent. The two lateral streaks also make a contribution to energy extraction for large shear ratios.
Global adjoint sensitivity analysis is performed to examine the receptivity and the inception of global instability. The most sensitive region to the point forcing is located immediately upstream of the roughness element and at the top edge of the separation region downstream. The wavemaker results show that the instability core is located in the reversed flow region for both varicose and sinuous modes. While the varicose mode displays one primary wavemaker region along the central streak, the sinuous mode shows two lateral wavemaker regions. The spatial growth of wavemaker is stronger for thicker roughness element. For thinner roughness, the strength of the spatial growth for varicose mode decreases faster than that for the sinuous mode.
The impact of η and Re h on the transition process associated with different instability characteristics is investigated by performing DNS. The results are compared to the transition diagram by Von Doenhoff & Braslow (1961), and transition features are seen to agree with their classification. For η = 1, the peak corresponding to the shedding of the primary hairpin vortices is obtained in both the energy and DMD spectra, in accordance with the eigenfrequency of the global varicose mode. As Re h increases, transition occurs closer to the roughness element, sinuous like breakdown is seen farther downstream, destabilizing the shear layer and promoting transition to turbulence. When Re h is sufficiently high, a fully-developed turbulent flow is established in both the inner and outer layers farther downstream. For η = 0.5, the sinuous wiggling of hairpin vortices becomes prominent in the near wake as Re h increases. Multiple peaks including the peaks corresponding to the varicose and sinuous instabilities are seen in the energy and DMD spectra. Stronger non-linear interactions between the hairpin vortices and the sinuous oscillations of the central streak are seen in the near wake. After the hairpin vortices break down, a sinuous mode associated with the wiggling of streaks persists farther downstream.