1. Introduction
The subcritical transition of leading-edge boundary layer near the attachment line of swept wings plays an important role in aerodynamics, which means the boundary layer may undergo transition to turbulence below the critical Reynolds number predicted by linear stability theory (LST). This phonomenon is especially critical because turbulent flow that starts at the leading edge of a swept wing can propagate downstream, affecting extensive regions of the wing’s chord and compromising its overall aerodynamic performance.
As the actual flow is three-dimensional (3-D) in nature, to simplify the problem, it is common to employ the swept Hiemenz boundary layer past a flat plate (Rosenhead Reference Rosenhead1963; Schlichting & Gersten Reference Schlichting and Gersten2017) as an approximation model for the actual three-dimensional boundary layer (Poll Reference Poll1979; Hall et al. Reference Hall, Malik and Poll1984; Theofilis Reference Theofilis1998; Theofilis et al. Reference Theofilis, Fedorov, Obrist and Ch. Dallmann2003) around the leading edge. Based on this model, the LST performed by Hall et al. (Reference Hall, Malik and Poll1984) gives a linear critical Reynolds number of
$Re_{crit}\approx 583.1$
, which is in good agreement with the previous experimental finding (Poll Reference Poll1979) as well as the numerical simulation of Spalart (Reference Spalart1988). However, in many experimental tests (Gaster Reference Gaster1967; Poll Reference Poll1979; Arnal et al. Reference Arnal, Juillen, Reneaux and Gasparian1997), transitions are often observed at a significantly lower value,
$Re_{tr} \approx 250$
, if the boundary layer is subject to sufficiently large external disturbulences. To understand the discrepancy between linear stability results and experimental findings, finite amplitude perturbations and nonlinear processes have to be taken into account. Obrist et al. (Reference Obrist, Henniger and Kleiser2012) and John et al. (Reference John, Obrist and Kleiser2014, Reference John, Obrist and Kleiser2016) carried out a direct numerical simulation (DNS) study on a incompressible swept Hiemenz boundary layer with a pair of stationary counter-rotating streamwise vortex-like structures of finite amplitude. A bypass transition scenario has been identified, which can explain the occurrence of subcritical transition in experiments. The initial pair of stationary counter-rotating vortex-like structures lead to the transient growth of streaks according to the lift-up effect, and then the damped primary vortices and streaks interacts with unsteady perturbations, causing secondary instabilities and leading to the final transition to turbulence.
However, the aforementioned conclusions are based on a simplified model in incompressible flow only. When compressible effects (such as Mach number, shock waves, wall temperature, etc.) are taken into account, the problem becomes significantly more complex. Based on previous studies (Theofilis et al. Reference Theofilis, Fedorov and Collis2006; Li & Choudhari Reference Li and Choudhari2008; Mack et al. Reference Mack, Schmid and Sesterhenn2008; Xi et al. Reference Xi, Ren and Fu2021a,Reference Xi, Ren, Wang and Fu b; Fedorov & Egorov Reference Fedorov and Egorov2022), for large sweep Mach numbers, the leading instability mode at the attachment line is inviscid in nature, while for lower sweep Mach numbers, the attachment-line instability exhibits the behaviours of viscous Tollmien–Schlichting waves. Detailed reviews for these research have been included in our previous studies (Xi et al. Reference Xi, Ren and Fu2021a,Reference Xi, Ren, Wang and Fu b) and the connection between the linear stability features of the flow and the issues discussed in this study is not particularly direct. Therefore, we will not elaborate on them here.
In fact, experimental investigations of high-speed attachment-line flow date back to 1959. Initially, Beckwith & Gallagher (Reference Beckwith and Gallagher1959) focused on the effects of sweep angles and heat flux along the attachment line in supersonic conditions. They detected the transition of attachment-line flow in their Mach 4.15 experiments, studying the effect of sweep angles over a relatively wide range. Later, Creel et al. (Reference Creel, Beckwith and Chen1986) and Chen et al. (Reference Chen, Creel and Beckwith1991) conducted experiments with a free stream Mach number of 3.5 and various sweep angles, also detecting transition along the attachment line and finding transition Reynolds numbers of approximately 650 (based on the boundary layer length scale at the leading edge). Skuratov & Fedorov (Reference Skuratov and Fedorov1991) performed similar tests to validate Creel et al.’s results. Murakami et al. (Reference Murakami, Stanewsky and Krogmann1996) studied hypersonic attachment-line flow in a Ludwieg-tube wind tunnel. In some conditions, the bypass scenario is the most possible reasons for the transition. During the experiments, without the end plates and trip wires, the attachment-line boundary layer can remain laminar along the entire attachment line.
Figure 1. Experimental model and the infrared measurements of the temperature distribution along the leading edge of the swept blunt body. Dots indicate the positions of the pressure sensors; pink represents high-temperature regions, while blue indicates low-temperature regions.
1.1. Recent experimental tests
Recently, experimental tests were performed over a swept blunt leading edge in the China Academy of Aerospace Aerodynamics (CAAA) to explore the potential leading-edge transition phenomena. The experimental model employed is the leading-edge section of a delta wing with a sweep angle of 45 o, as illustrated in the accompanying figure 1. The model has a width of 300 mm, and the length of the leading edge along the swept direction is approximately 425 mm. Due to the symmetrical nature of the model, only a half-model was used in the experimental investigation. The main wing structure is supported by steel along with a steel rear support bracket for stabilisation. The leading-edge section, specifically the blunt leading edge, is constructed from polyether ether ketone (PEEK) material, known for its high emissivity and low thermal conductivity. The use of PEEK material is specifically chosen to meet the requirements of infrared thermographic imaging. The thickness of the blunt leading edge is
$40\, \text {mm}$
, comprising two semicircular arcs of
$17.5 \,\text {mm}$
radius combined with a 5 mm thick flat plate in the middle.
The wind tunnel experiment was conducted in the conventional hypersonic wind tunnel at the CAAA. This facility operates as a blowdown free jet wind tunnel. The nozzle exit diameter is 0.5 m, with a simulated Mach number range of 4–8. The test section features a 350 mm diameter observation window, facilitating infrared thermographic imaging. The observation window is made of silicon glass coated with anti-reflective and infrared-transmitting films. The infrared thermographic camera uses a refrigerated detector with a spectral range of 3.7–4.8
$\mu \mathrm{m}$
, a thermal sensitivity of
$\lt 25 \,\text {mK}$
and a measurement temperature accuracy of
$\pm 1\,^{\circ }$
C. The system had been successfully used in a previous study (Sha et al. Reference Sha, Yue, Feng, Xiangjiang and Qing2020). Eight pressure sensors are placed along the attachment line of the blunt leading edge, uniformly distributed, in a similar manner as in a previous study (Yao et al. Reference Yao, Min, Xie, Li, Lin, Yang and Duan2018). The operation conditions for the wind tunnel test are set at a Mach number of 6.00, a total pressure of 1.00 MPa and a Reynolds number (based on a characteristic length of 1 m) of
$1.80\times 10^7$
. The temperature of the experimental gas (air) is maintained at 56.58 K. Prior to the wind tunnel test, the model’s surface temperature is at room temperature. Once the wind tunnel is activated, the infrared camera begins capturing data. After establishing a stable flow field, the model is quickly inserted into the hypersonic flow field using an insertion mechanism, and the surface temperature changes are recorded in real time. After stabilising for 5 s, the model is withdrawn from the flow field, concluding the experiment. The temperature distribution of the model surface at the 5-second mark within the flow field is recorded and analysed.
The primary finding from the experiment is the high sensitivity of the attachment-line boundary layer flow to wall perturbations, especially near the pressure sensor installation sites. When pressure sensors are mounted at the attachment-line position on the leading edge of a swept blunt body model, the model’s surface is no longer smooth. Due to unavoidable errors during the experiment, effective roughness elements – such as small protrusions or depressions – are formed along the attachment line. Experiments conducted with this configuration have demonstrated that perturbations induced by these roughness elements can effectively trigger the transition of the attachment-line boundary layer to turbulence, as illustrated in figure 1.
1.2. Motivation
Previous studies indicate that the phenomenon of subcritical transition is highly significant in the context of attachment-line flows. However, most of these studies have been limited to incompressible flows or confined to linear analysis, leaving a significant gap in the understanding of subcritical compressible flows. Furthermore, even in incompressible flow scenarios, existing computational analyses have often employed simplified models or introduced artificial perturbations to facilitate numerical studies, raising questions about their validity under actual conditions. Therefore, it is imperative to conduct numerical investigations of the three-dimensional boundary layer at the leading edge of compressible blunt bodies under relatively realistic conditions. In this study, we perform numerical simulations of transitional high-speed attachment-line boundary layers that develop from finite amplitude initial perturbations caused by roughness elements with different heights. Using LST, we find that the attachment-line boundary layer is subcritical with respect to small initial perturbations. When considering the height of the roughness elements, being subcritical implies that the roughness height is insufficient to independently induce a transition. Both subcritical and supercritical roughness elements are employed in this study. These simulations correspond to experimental observations of roughness-induced transition over a real blunt configuration, without assuming an infinite span. Unlike typical transitional studies, we have calculated the complete transition to turbulence over a real configuration.
Our primary aim in this paper (Part 1) is to investigate the physical mechanisms behind the roughness-induced transition in three-dimensional attachment-line boundary layers at the leading edge. In this study, we focus specifically on the transitional processes, while a companion paper (Part 2) delves into the turbulence characteristics further downstream. The structure of this paper is as follows. In § 2, the governing equations and numerical methodologies are introduced in detail. The core findings and detailed analysis for transitional three-dimensional boundary layers are presented in the following sections. Section 3 describes the general features of the flow fields, with a particular focus on the mean flow characteristics. In § 4, the transition mechanisms are discussed in detail and the conclusions as well as some discussions are given in § 5.
2. Methodology
2.1. Governing equations
The governing equations for all simulations in this work are the dimensionless compressible Navier–Stokes (NS) equations for a Newtonian fluid, which can be written as
\begin{align} \frac {\partial Q}{\partial t}+\frac {\partial F_j}{\partial x_j} + \frac {\partial F_{j}^{v}}{\partial x_j}=0,\qquad\qquad \end{align}
\begin{align} Q = \left [ \rho ,{\rho {u_1}},{\rho {u_2}},{\rho {u_3}},{{E_t}} \right ]^{T}\!\!,\qquad\quad \end{align}
\begin{align} {F_j} = \left [ {\begin{array}{c} {\rho {u_j}} \\ {\rho {u_1}{u_j} + p{\delta _{1j}}} \\ {\rho {u_2}{u_j} + p{\delta _{2j}}} \\ {\rho {u_3}{u_j} + p{\delta _{3j}}} \\ {\left ( {{E_t} + p} \right ){u_j}} \end{array}} \right ],\quad {F_{j}^{v}} = \left [ {\begin{array}{c} 0 \\ {{\tau _{1j}}} \\ {{\tau _{2j}}} \\ {{\tau _{3j}}} \\ {{\tau _{jk}}{u_k} - {q_j}} \end{array}} \right ].\\[6pt] \nonumber \end{align}
Throughout this work, the coordinates
$x_i, (i =1, 2, 3)$
are referred to as
$x, y, z$
, respectively, with corresponding velocity components
$u_1 = u, u_2 = v, u_3 = w$
. Here,
$F_j$
and
$F_j^{v}$
stand for the inviscid and viscous flux. The total energy
$E_t$
and the viscous stress
$\tau _{ij}$
are given respectively as
\begin{equation} \begin{aligned} E_{t}&=\rho \left (\frac {T}{\gamma (\gamma -1) M^{2}_{\infty }}+\frac {u_{k} u_{k}}{2}\right ), \\ \tau _{i j}&=\frac {\mu }{Re_{\infty }}\left (\frac {\partial u_{i}}{\partial x_{j}}+\frac {\partial u_{j}}{\partial x_{i}}-\frac {2}{3} \delta _{i j} \frac {\partial u_{k}}{\partial x_{k}}\right ). \end{aligned} \end{equation}
The pressure
$p$
and heat flux
$q_i$
are obtained from
\begin{equation} p=\frac {\rho T}{\gamma M_{\infty }^{2}}, \quad q_{i}=-\frac {\mu }{(\gamma -1) M_{\infty }^{2} Re_{\infty } Pr} \frac {\partial T}{\partial x_{i}}. \end{equation}
The viscosity is calculated using the Sutherland law,
\begin{equation} \mu = T^{3/2}\frac {T_{\infty } + C}{T\cdot T_{\infty } + C}, \end{equation}
with
$C = 110.4\,\mathrm{K}$
. The free stream Reynolds number
$Re_{\infty }$
, Mach number
$M_{\infty }$
and Prandtl number
$Pr$
are defined as
\begin{equation} Re_{\infty } = \frac {\rho _{\infty }^* U_{\infty }^* l_0^*}{\mu _{\infty }^*}, \quad M_{\infty } = \frac {U_{\infty }^*}{\sqrt {\gamma R_g^* T_{\infty }^*}}, \quad Pr = 0.72, \end{equation}
where
$\rho _{\infty }^*$
,
$U_{\infty }^*$
,
$T_{\infty }^*$
and
$\mu _{\infty }^*$
stand for the free stream density, velocity, temperature and viscosity, respectively. Here,
$R_g^* = 287 \text {J}/(\text {K}\cdot \text {Kg})$
represents the gas constant and
$\gamma$
stands for the ratio of specific heat. The length scale
$l_0^*$
is chosen as
$1$
millimetre in this research. The
$*$
denotes dimensional flow parameters.
2.2. Numerical method
Two solvers have been employed in this study. The first code we use to perform computations is the high-order finite difference code developed recently at Tsinghua University. A shock-fitting (S-F) method (Zhong Reference Zhong1998) is used to compute steady hypersonic viscous flow together with the high-order accurate non-compact finite differences methods. The fifth-order upwind scheme (for inviscid flux
$F_{j}$
) and the sixth-order centre scheme (for viscous flux
$F_j^v$
), are used to compute the flow field. A fourth-order Runge–Kutta method is applied for the time integration, and the simulations are performed until the maximum residual reaches a small value of the order of
$10^{-15}$
. A full implicit scheme can also be used for fast convergence. Validations of the code and some applications for calorically perfect gas and thermal-chemical non-equilibrium flow can be found in our previous studies (Xi et al. Reference Xi, Ren and Fu2021a,Reference Xi, Ren, Wang and Fu
b; Chen et al. Reference Chen, Xi, Ren and Fu2022). The solver is used mainly to determine the location of the leading shock and give a high-quality initial field.
The second code, used in this study, is a well-validated fluid dynamic shock capture (S-C) solver OPENCFD, developed by Li et al. (Reference Li, Fu and Ma2008), which is mainly used to simulate the whole transition/turbulent processes. The code has been validated and verified in previous studies (Li et al. (Reference Li, Fu and Ma2008); Li et al. (Reference Li, Fu and Ma2010); Liang et al. (Reference Liang, Li, Fu and Ma2010)). For three-dimensional calculations presented in this study, a hybrid high-order finite difference scheme, including the seventh-order upwind scheme, fifth-order and seventh-order WENO schemes (Jiang & Shu Reference Jiang and Shu1996), together with a shock sensor (Dang et al. Reference Dang, Liu, Guo, Duan and Li2022), is used for the inviscid flux in the characteristic form. Based on that formula, during the calculation, more than 98 % of the regions use the linear seventh-order upwind scheme and only a few regions corresponding to discontinuities use the nonlinear WENO schemes, which greatly increase the calculation efficiency. A standard sixth-order central difference scheme is used for viscous flux. Time advancement is achieved by the explicit third-order total variation diminishing Runge–Kutta scheme.
2.3. Models
Figure 2. Schematics of the swept blunt leading edge used for numerical simulations. The blue lines indicate the leading shock.
Our computational model is designed based on the experimental configuration introduced in § 1 and is shown in figure 2. The model is the front part of a delta wing with a swept angle
$\Lambda$
of 45 o. The thickness of the wing is
$2L_4 = 40$
. The computational model can be likened to a sandwich-like configuration, where the top and bottom layers consist of semicircles with a radius of
$R_1 = 17.5$
and the intermediate layer is a flat plate with a width of
$L_3 = 5$
. Together, these three layers form the complete swept blunt body configuration. The roughness element is located at
$z = L_2 = 40$
, at the centre of the leading plate. The radius of the roughness is
$R_2 = 1$
. The length of the whole model is designed as
$L_1 = 200$
. Based on experiments, the surface temperature is set to
$T_w^* = 370\text {K}$
, other relative flow parameters are listed in table 1. In this study, we have established a coordinate system, as in figure 2 wherein the
$z$
-axis aligns with the leading edge of the swept blunt model, coinciding with the attachment line and extending in the corresponding spanwise direction. The normal direction on the corresponding attachment line and swept blunt body is defined as the
$y$
-axis. Finally, the
$x$
-axis is defined to complete the typical Cartesian coordinate system in conjunction with these two axes. As usual, a body-fitted coordinate
$(\xi ,\eta ,z)$
is also established with the same spanwise direction as the Cartesian coordinate system, while the
$\xi$
-axis is defined along the chordwise direction and the
$\eta$
-axis is defined along the surface normal directions.
Table 1. Basic parameters for flow and roughness at basic grid.
$N_k$
is the number of wall normal points for
$0 \leqslant y \leqslant k_h$
.
$\delta _{bl}^* = 0.2\, \text {mm}$
is the thickness of the laminar boundary layer at the attachment-line boundary layer.
As previous analysis around the attachment-line boundary layers, we define the sweep Mach number
$M_s$
and the sweep Reynolds number
$Re_s$
as
\begin{equation} Re_s = \frac {w^*_{\infty } \delta _{al}^* }{\nu ^*_r} \approx 714, \quad M_s = \frac {w^*_{\infty }}{c^*_s} \approx 2, \end{equation}
based on the length scale
$\delta _{al}^* = \sqrt {\nu ^*_r {\partial u^*_e}/{\partial x^*}}$
at the exact attachment line
$x^*=0$
and the variables outside of the attachment-line boundary layer. Here,
$c^*_s$
is the sound speed after the leading shock,
$\nu ^*_r$
stands for the kinematic viscosity at recovery temperature
$T^*_r\approx 433\,\mathrm{K}$
and the temperature at the edge of the leading attachment-line boundary layer is
$T_{at, e}^*\approx 260\,\mathrm{K}$
. The value
${\partial u^*_e}/{\partial x^*}$
at the exact attachment line
$x^*=0$
is not known a priori for the present case, the potential flow around a circular cylinder with equivalent radius
$L_4$
is thus used to evaluate the derivative. By using the linear stability theory over two-dimensional domains, the neutral surface of the most dangerous discrete mode are presented in figure 3, over a
$Re_s {-} M_s {-} \beta$
coordinate. The base states for the stability analysis are computed using a combination of solving the Euler equations and the three-dimensional boundary layer equations. Detailed settings and calculations can be found in many previous studies (Theofilis et al. Reference Theofilis, Fedorov and Collis2006; Gennaro et al. Reference Gennaro, Rodríguez, Medeiros and Theofilis2013; Xi et al. Reference Xi, Ren and Fu2021a). Here,
$\beta$
is the normal spanwise wavenumbers. The red line in figure 3 indicates the case on which we focus in this study. It is found that the critical Reynolds number increases with increasing sweep Mach number, which indicates that the leading viscous mode (Görtler–Hämmerlin mode) is supressed by the compressible effects. Also, it is clearly shown that the present case (the red line) is located at the stable or subcritical region, which means that the transition to turbulence at the present case is not triggered by a linear instability.
Figure 3. Neutral surface of the most dangerous discrete temporal mode over the
$Re_s{-}M_s{-}\beta$
plane. The growth rate space is divided into stable and unstable regions by the neutral surface.
2.4. Roughness elements
Figure 4.
$(a)$
Grid distributions arount the roughness with the roughness height
$k_h = 0.1$
mm in full resolution.
$(b)$
Shape of roughness in two cross-sections.
In this study, the roughness element, depicted in figure 4, is designed to simulate a pressure sensor with a circular disk configuration. The geometry of this element is described in polar coordinates
$(r, \phi )$
, where the shape of the shoulders is defined using a hyperbolic tangent function, similar to the methods employed in previous studies (Kurz & Kloker Reference Kurz and Kloker2014, Reference Kurz and Kloker2016). The function is defined as
\begin{equation} h(r, \phi ) =\frac {k_h}{2} + \frac {k_h}{2}\tanh \left [\frac {S_r}{k_h}\left (\frac {d}{2} - r\right )\right ], \end{equation}
with
$k_h$
and
$d$
being the height and diameter of the roughness, respectively. The slope factor
$S_r$
is set to 1.0 for all cases in the present study. In general, the centre of the roughness is located at the points
$(x, z) = (0, 40)$
, the diameter is
$d = 2R_2 = 2$
.
Another important parameter for roughness-induced transition is the roughness reynolds number
$Re_{kk}$
, characterised based on the height
$(k_h)$
and the velocity
$(w)$
in the undisturbed laminar flow with respect to the position of the roughness. This roughness Reynolds number is defined as a function of
\begin{equation} Re_{kk} = \frac {\rho (k_h) w(k_h) k_h }{\mu (k_h)}, \end{equation}
and listed in table 1 based on the laminar boundary layer.
2.5. Simulation strategy
Figure 5. Outline of calculation processes. S-F, shock fitting; S-C, shock capture.
The calculation process for this kind of problem is divided into three steps, with two kinds of compressible solvers (Li et al. Reference Li, Fu and Ma2008; Xi et al. Reference Xi, Ren, Wang and Fu2021 b) as shown in figure 5.
-
Step 1. Assuming the incoming flow has reached an asymptotic state, a two-dimensional calculation with infinity span assumption
$(\partial /\partial z = 0)$
is performed first using a shock-fitting solver. With the exact location of the shock revealed, alignment and clustering of the mesh along the bow shockwave can be easily achieved for the following shock-capture calculation. -
Step 2. To diminish the numerical perturbations between the two solvers, a three-dimensional domain is further designed for pre-calculation with a periodic boundary condition along the spanwise direction, with the newly built grid and the initial field from the fitting solver.
-
Step 3. When the calculation is converged, the solution from the middle slice of the periodic three-dimensional domain is used for the fully three-dimensional calculation. As the boundary layer develops along the attachment line and the chordwise direction, non-reflection outlet boundary conditions are used further downstream along the attachment line and the chordwise direction. Away from the surface, as the shock is embeded in the computational domain, free stream boundary conditions are used at the outside.
To closely mimic the conditions observed during experimental tests, the generation of unsteady perturbations is implemented in two distinct phases. First, the free stream noise level in the wind tunnel experiments is approximately 1.5 %–3 % based on pressure fluctuations. Second, it was clearly observed in the experiments that disturbances were generated near the upstream pressure sensors, which did not directly trigger transition. However, these disturbances have the potential to influence the flow state downstream. Therefore, in the first phase, random velocity perturbations, with maximum amplitude constituting approximately
$2\,\%$
of the free stream velocity, are introduced upstream of the leading shock waves. This procedure aims to replicate the perturbations measured in wind tunnel experiments. In the second phase, to simulate the perturbations inherent to upstream boundary layers along the attachment line, random wall normal blowing and suction are executed via a hole on the wall. These perturbations also possess a maximum amplitude of roughly
$2\,\%$
of the free stream velocity. The centre of the specified hole is positioned at coordinates
$(x, z) = (0, 30)$
and the hole is defined with a radius of
$2$
.
In the computational analyses conducted within the scope of this study, two distinct cases were examined. In the first scenario, characterised by a roughness element height of 0.1, unsteady perturbations were deliberately introduced to facilitate the onset of transition. Otherwise, the transition would not occur within the wake flow induced by the roughness element. Conversely, the scenario involving a roughness element height of 0.2 presented a fundamentally different dynamic. The inherent absolute instability associated with this configuration amplified numerical perturbations during the computation. This natural amplification process effectively initiates the transition, obviating the need for the introduction of external perturbations.
2.6. Computational grids
2.6.1. Grid distributions for shock-fitting solver
The basic grid number used for shock-fitting simulations is
$N_{\xi } \times N_{\eta } = 801 \times 401$
(for Step 1). The distribution of the grid points in the wall normal direction is controlled through a function that provides clustering towards the wall, with two parameters
$h_{im}$
and
$\sigma _s$
. The distribution function, which maps
$\eta$
to
$h$
, can be expressed as
\begin{equation} \begin{aligned} h = H_{shk}\frac {a_y (1 + Y)}{(b_y - Y)},\!\quad b_y =1 + 2 a_y,\!\quad a_y =\frac {h_{im}}{1 - 2 h_{im}},\!\quad Y =2 \frac {\left [1 - \tanh {(\sigma _s)}\right ]\frac {1+\eta }{2}}{ 1 - \tanh {\left (\sigma _s\frac {1+\eta }{2}\right )}} - 1, \end{aligned} \end{equation}
where
$H_{shk}$
is the local shock height and is solved as a dependent variable with the flow field in shock-fitting methods;
$\eta$
is a uniform grid distribution along the region
$ [-1.0,1.0 ]$
;
$h$
is the actual wall normal grid distributions. The values of
$h_{im}$
and
$\sigma _s$
are chosen to be
$0.3$
and
$0.95$
for the fitting simulations presented in this paper. Along the chordwise direction, at the wall surface, the surface grid
$s(x,y)$
is clustered at the round head with the function
\begin{equation} \frac {s}{S} = \frac {a_{\xi }(1 - \xi )}{b_{\xi } - \xi },\quad b_{\xi }= 1 + 2 a_{\xi },\quad a_{\xi }=\frac {s_{im}}{1 - 2s_{im}}, \end{equation}
where
$s_{im}$
is chosen to be
$0.2$
for the simulations based on shock-fitting methods,
$s$
is the local surface curve length along the model surface and
$S$
is the total model surface curve length. Note that the surface grid does not change during the calculation.
2.6.2. Grid distributions for shock-capture solver
The basic grid numbers used for shock-capture simulations are
$N_{\xi } \times N_{\eta } \times N_{z} = 801 \times 401 \times 8$
(for Step 2) and
$N_{\xi } \times N_{\eta } \times N_{z} =2401 \times 401 \times 4401$
(for Step 3). More grids points with
$N_{\xi } \times N_{\eta } \times N_{z} =3001 \times 601 \times 6601$
(for Step 3) are used to validate the independence of the solutions, as shown in figure 6. Unlike the fitting solver, for the capture solver, the wall normal grid needs clustering towards the wall as well as the shock region. Therefore, the wall normal grid is divided into three parts. The first part is at the region
$ [0, 0.105H_{shk} ]$
, where
$H_{shk}$
is the local shock height. The second part is a transition region which connects the near-wall region and the shock region, and is located at the region
$ [0.105H_{shk}, 0.945H_{shk} ]$
. The final region stands for the part used to capture the shock and is located at the region
$ [0.945H_{shk},1.05H_{shk} ]$
. As the local shock height
$H_{shk}$
is solved by the fitting solver in advance, the grid lines are adapted well to both the body and the shock shape. Unlike the usual calculations of the transitional/turbulent boundary layer, the shock regions in the present simulation have been taken into account with enough accuracy.
Figure 6. Regions of flow separations in a cross-cut plane through the centre of a roughness for case H0100, marked by contours of
$w=-0.001$
, for basic and high resolutions. The axis is stretched for clarity.
For the basic grid, the first region has 201 grid points, with the distributions
$h_1$
as
\begin{equation} h_1 = \frac {21}{1600}\frac {1 + \eta _1}{3.5 -\eta _1}H_{shk}, \end{equation}
where
$\eta _1$
is the uniform distributions in region
$ [-1,1 ]$
. The third region has 101 grid points with a uniform distribution. The second region has 101 grid points and is used to link grids with different grid spacings using a Hermite function obtained by imposing C3 continuity of the resulting stretching function. The grid distributions and relative grid spacing is shown in figure 7. For the S
tep 2 calculation of the periodic box, eight grids are used to cover a spanwise region of length
$20$
. For the final full three-dimensional calculation, 4401 grids are distributed uniformly along the spanwise direction. The roughness element and unsteady hole are captured with 101 grids points.
By using the above approach, for the basic grid numbers, around the attachment-line boundary layer,
$125$
grid points are located inside the boundary layer(based on laminar boundary layer thickness) in the wall normal
$\eta$
-axis direction. Approximately
$230$
grid points are located at the leading-edge plate along the
$x$
-axis. As the boundary layers develop along the chordwise direction as well,
$136$
grid points are located inside the boundary layer at the chordwise outlet boundary. Detailed grid numbers that are used to resolve roughness are shown in table 1. Comparing with the existing research (Mayer et al. Reference Mayer, von Terzi and Fasel2011; De Tullio et al. Reference De Tullio, Paredes, Sandham and Theofilis2013; De Tullio & Sandham Reference De Tullio and Sandham2015; Groskopf & Kloker Reference Groskopf and Kloker2016; Di Giovanni & Stemmer Reference Di Giovanni and Stemmer2018; Hader & Fasel Reference Hader and Fasel2019; Shrestha & Candler Reference Shrestha and Candler2019) on numerical simulations of high-speed boundary layer transitions, these grid numbers and distributions can reveal the detailed behaviour of linear and nonlinear waves. To further identify the possible features of the turbulence further downstream, grids points along the wall-normal directions are increased to 601 in the final simulations presented in this paper. The grid sizes in wall units in the turbulent region are shown in figure 8 and listed in table 2. Based on that grid size, as there are no additional stress and heat flux terms in the present simulations, resolutions for an implicit large eddy simulation (ILES) can be recognised for the possible turbulent region.
Figure 7. Basic grid distributions with 401 grid points along the wall normal directions.
$(a)$
Scaled grid distributions.
$(b)$
Scaled grid spacing along the wall normal direction.
Figure 8. Distributions of grid sizes in wall units at the turbulent boundary region.
Table 2. Grid points and maximum grid sizes in wall units at the turbulent boundary region.
3. General features of the flow fields
3.1. Verification of the asymptotic assumption
In the current study, as the three-dimensional boundary layer develops along the spanwise
$z$
-direction, the commonly employed asymptotic assumption (
$\partial /\partial z = 0$
) for attachment-line and cross-flow analyses can be verified by examining the profiles upstream of the roughness. Taking the case H0100 as an example, the profiles of major variables such as density
$\rho$
, spanwise velocity
$w$
and temperature
$T$
in the attachment-line region
$x\in [-2.5,2.5]$
are depicted in figures 9 and 10. A perfect alignment of these profiles is observed, which confirms two critical points. First, the exact correspondence of results obtained from the S-C and S-F solvers corroborates the precision of the solvers used in this research. Second, and more significantly, because of the agreements of the profiles at inlet
$z=0$
and at
$z=18.2$
, the asymptotic assumptions are still applicable to a spatially developing attachment-line boundary layer.
Figure 9. Profiles of major variables at the
$x=0$
plane, for case H0100.
Figure 10. Profiles of major variables at the
$x=2.5$
plane, for case H0100.
Upon broadening our perspective further, it becomes apparent that even in regions where cross-flow effects are pronounced, this assumption remains valid, as shown in figure 11. However, it is crucial to emphasise that the upstream flow has reached an asymptotic state. This assumption is applicable to fully developed laminar boundary layers, provided that the spanwise extent is sufficiently large.
Figure 11. Profiles of major variables at three locations (
$30^{\circ}, 60^{\circ}, 90^{\circ}$
) over a cylinder surface, for case H0100.
We further compare the profiles at
$x=2.5$
, a little bit away from the attachment line at
$x=0$
, as presented in figure 12. As the flow develops along the chordwise direction, the basic feature of the variables are kept the same, but the boundary layer has become slightly thinner.
Figure 12. Comparison of variables at the
$x=0$
and
$x=2.5$
planes, for case H0100.
Figure 13. Instantaneous iso-surface of
$\lambda _2 = -0.035$
, colour indicates
$w$
, for the case H0100.
Figure 14. Instantaneous iso-surface of
$\lambda _2 = -0.035$
, colour indicates
$w$
, for the case H0200.
Figure 15. Surface heat fluxes
${\theta }_{tw}$
distributions for two cases:
$(a)$
for H0100;
$(b)$
for H0200.
3.2. General features of the instantaneous flow fields
Figures 13 and 14 illustrate the general characteristics of the flow field, visualised through the iso-surfaces of
$\lambda _2=-0.035$
in the instantaneous fields, highlighting the underlying vortex structures. The iso-surfaces are coloured with spanwise velocity
$w$
. The whole flow fields can be divided into three parts. The first part is the roughness region, where the initial laminar flow is disturbed by the surface roughness element. Vortex structures develop around the roughness, displaying longitudinal streaks that roll up in the wake of the roughness element and form hairpin-like vortices. These hairpin vortices eventually break down into smaller vortices. The breakdown of these vortex structures marks the onset of turbulence, with smaller vortex structures becoming apparent along the attachment line. The second part is the transitional region, in which the initial turbulences at the attachment line develop along the spanwise direction as well as the chordwise direction. At the very beginning, the turbulences are located around the attachment line and the turbulent region expands along the chordwise direction slightly. Then, as the flows develop further downstream, the turbulent structures are flushed from the attachment-line region to the chordwise outlet. The final part is the fully turbulent region, where the fully developed turbulences cover the whole region of the leading-blunt body, ranging from the attachment line to the chordwise outflow.
When the height of the roughness element is increased from 0.1 to 0.2, there are significant differences in the vortex structures formed behind the roughness, as shown in figures 13 and 14. For higher roughness element, the horseshoe vortex ahead of the roughness is much more pronounced and the relative induced streaks can directly induce transition downstream, at the side region. Close to the symmetry plane and downstream of the roughness, one can notice the presence of large streamwise streaks. Further downstream, these streaks induce the formation of hairpin-shaped structures and also lead to the transition. To provide a more intuitive and visual representation of the aforementioned process, several animations are included in the supplementary material. Additionally, these two different characteristics can be seen more directly in the contours of surface average heat fluxes
$\theta _{tw}$
and skin frictions
$\overline {\tau }_w$
, as depicted in figures 15 and 16. These metrics essentially serve as ‘footprints’ of the boundary layer dynamics, providing insights into the complex interactions and flow structures present within the boundary layers. Here, as the usual boundary layers previously, we define the velocity
$\overline {u}^{+}$
, based on inner scale, as
\begin{equation} \left . \begin{aligned} h^+ &= \frac {\overline {\rho }_w \overline {u}_{\tau } h}{\overline {\mu }_w}, \quad \overline {u}^+ = \frac {|\overline {u}_p|}{\overline {u}_\tau }, \quad |\overline {u}_p| = \sqrt {\overline {u}_{\xi }^2 + \overline {w}^2}, \\ \overline {u}_{\tau } &= \sqrt {\frac {\overline {\tau }_w}{\overline {\rho }_w}}, \quad \overline {\tau }_w = \frac {\overline {\mu }}{Re}\left . \frac {\partial \overline {u}_p}{\partial h} \right |_{h=0}, \end{aligned} \right \} \end{equation}
where
$\overline {u}_{p}$
is the velocity parallel to the surface. The skin-friction coefficient
$\overline {C_f}$
and surface heat-flux
${\theta }_{tw}$
for this kind of flow are defined as
\begin{equation} \overline {C_f} = \frac {2 \overline {\mu }_w^*}{\rho _{\infty }^* U_{\infty }^{*2}} = \frac {2 \overline {\mu }_w}{Re}\frac {\partial \overline {u}_{p}}{\partial h} = 2 \overline {\tau }_w,\ {\theta }_{tw} = - \left |\kappa \nabla \overline {T} \cdot \boldsymbol {n}\right |. \end{equation}
The derivatives of surface normals, denoted as
${\partial }/{\partial h}$
, for arbitrary variables
$f_{\psi }$
, are determined through a two-step process. Initially, the gradients of the variables
$f_{\psi }$
are computed using the identical scheme adopted for the calculation of viscous fluxes during the simulations. Subsequently, the derivatives of the surface normals
${\partial f_{\psi }}/{\partial h}$
are obtained by projecting the calculated gradients
$\nabla f_{\psi }$
onto the surface normal vectors
$\boldsymbol {n}$
.
Figure 16. Surface skin friction
$\overline {\tau }_w$
distributions for two cases:
$(a)$
for H0100;
$(b)$
for H0200.
Figure 17. Density gradient magnitude contours of the case H0200, at attachment-line plane
$x=0$
. The red line stands for the computational domain.
3.3. General influences of the roughness height on mean flow
The magnitude contour of average density gradients for case H0200, at the attachment-line plane, is depicted in figure 17. This illustration provides a comprehensive view of the general flow field characteristics for both cases. The presence of surface roughness induces a shock slightly ahead of the roughness. As this shock evolves away from the surface and progresses downstream, it shapes into a curved shock surface under the influence of the incoming flow. The interaction of this induced shock with the leading shock of the blunt body, followed by its reflection back into the boundary layer downstream, results in a noticeable deformation of the leading shock. Subsequent to the roughness-induced shock, the compressed fluids undergo expansion and acceleration, leading to the formation of a recompression shock at the roughness’s tail along the
$z$
-direction. Meanwhile, a shear layer develops behind the roughness and the recompression shock once again impinges on the leading shock, reflecting back into the boundary layers. When the flow evolves further downstream, the high-shear region at the outside of the boundary layer becomes much weaker, as reflected as the decrease of magnitudes for the density gradient
$\left | \nabla \rho \right |$
.
Figure 18. Line plots of average spanwise velocity
$\overline {w}$
around the roughness for the
$(a)$
H0100 and
$(b)$
H0200 cases. The red and blue lines stand for the wall surfaces and seperation bubbles, respectively. The shape of the separation bubble is approximated by the
$\overline {w}$
velocity contour of
$-0.0001$
. The density gradient magnitude
$|\nabla \overline {\rho }|$
contours are shown against the background.
Figures 18 and 19 show the distributions of mean velocity and temperature along the attachment line in the wall-normal direction, and the density gradient magnitude
$|\nabla \overline {\rho }|$
contours are shown against the background. The roughness-induced oblique shocks as well as the high-shear region at the edges of the boundary layers can be seen clearly through these shaded contours. Additionally, the size and specific location of the corresponding separation bubbles are indicated by blue lines in the figures and are clearly compared in figure 20. It is evident from the figure that the separation bubbles induced by small roughness elements are lower in height compared with those induced by large roughness elements, but extend farther downstream. Conversely, the separation bubbles induced by larger roughness elements extend farther upstream. The corresponding velocity and temperature profiles highlight the approximate location of the shear layer and illustrate the process through which low-speed fluid, due to the wakes induced by the roughness elements as well as the lift-up effect, is elevated away from the wall.
Figure 19. Line plots of average temperature
$(\overline {T} - T_w)/(T_{\infty } - T_w)$
at the attachment line around the roughness for the
$(a)$
H0100 and
$(b)$
H0200 cases. The red and blue lines stand for the wall surfaces and seperation bubbles, respectively. The shape of the separation bubble is approximated by the
$\overline {w}$
velocity contour of
$-0.0001$
. The density gradient magnitude
$|\nabla \overline {\rho }|$
contours are shown against the background.
Figure 20. Comparison of separation bubble shapes corresponding to different heights of roughness elements. The red and blue lines stand for the wall surfaces and seperation bubbles, respectively. The shapes of the separation bubble are approximated by the
$\overline {w}$
velocity contour of
$-0.0001$
.
Figure 21. Contours for average spanwise velocity
$\overline {w}$
at different spanwise locations (
$z$
) for the H0100 case. The thick black lines indicate the contours of the roughness element.
Figure 22. Contours for average spanwise velocity
$\overline {w}$
at different spanwise locations (
$z$
) for the H0200 case. The thick black lines indicate the contours of the roughness element.
The contours of the average spanwise velocity
$\overline {w}$
are shown in figures 21 and 22 for better understanding of the general features of the flow fields. It is clearly observed that regardless of the height of the roughness element, the wake exhibits two consistent patterns. First, there is a noticeable convergence of flow towards the centre (
$x=0$
) behind the roughness element. Second, the vortices in the regions adjacent to the roughness element (
$x\lt -1$
or
$x\gt 1$
) move in their respective chordwise directions. For the H0100 case, due to the smaller size of the roughness element, the vortex formed behind it is not strong, and thus the flow is predominantly influenced by the vortices generated at the sides of the roughness element, shown in figure 21. In the H0200 case, the induced streamwise vortices lift up low-momentum fluid from the near-wall region and give rise to a low-velocity streak away from the wall, as shown in figure 22. Once these low-velocity streaks form, they significantly alter the flow field distribution across the section, making the flow more unstable. Additionally, the induced stronger horseshoe vortices by the roughness element increase the three-dimensionality of the flow downstream near the edges of the roughness element (in the regions where
$x\lt -1$
or
$x\gt 1$
). The velocity profiles in the corresponding regions exhibit certain characteristics very close to the roughness elements. These characteristics are quite similar to the classic co-rotating rollover (or half-mushroom-like) structures formed by saturated cross-flow vortices (Saric et al. Reference Saric, Reed and White2003; Chen et al. Reference Chen, Xi, Ren and Fu2022). This may lead to earlier transition in the associated flow regions compared with the wake region.
Figure 23. Iso-surface of spanwise velocity
$\overline {w} = 0$
for the
$(a)$
H0100 and
$(b)$
H0200 cases. Limiting streamlines along the surfaces for the
$(c)$
H0100 and
$(d)$
H0200 cases. Spatial streamlines around the roughness for the
$(e)$
H0100 and
$(f)$
H0200 cases. The colours of the lines are used to distinguish the different height of the streamlines upstream. The seeds of the white lines are located at
$h=0.02$
, while those of the light blue lines are located at
$h=0.05$
.
To capture the fundamental three-dimensional nature of the flow, we present the distribution of limiting streamlines near the roughness elements and the iso-surfaces (figure 23) where the average spanwise velocity
$\overline {w}$
is zero, which are used to characterise the three-dimensional flow properties. Upstream of the roughness elements, a distinct recirculation region is observed. For lower roughness elements, this recirculation region forms a relatively complete separation bubble, as evidenced by the upstream recirculation streamlines. For higher roughness elements, in addition to a more pronounced primary separation region, a secondary separation line is clearly visible in the limiting streamlines near the roughness elements. As the flow progresses downstream, different degrees of flow convergence and vortex characteristics appear behind the roughness elements. The downstream streamlines indicate that the wake of the lower roughness element quickly returns to a more orderly state, which is also reflected by the absence of distinct vortex structures in the wake. The corresponding spatial streamlines in figure 23
$(e)$
become very orderly, with no noticeable twists or entanglements. In contrast, the higher roughness element exhibits completely different characteristics, with the formation of a pair of counter-rotating vortices in the wake, clearly reflected in the corresponding separation lines on the wall in the limiting streamlines in figure 23
$(d)$
. Also, noticeable twists or entanglements of the spatial streamlines can be found in figure 23
$(f)$
.
4. Mechanisms of the roughness-induced transition
In the subsequent sections, we further investigate the transition mechanisms induced by roughness elements using a robust statistical analysis framework. This approach is designed to systematically reveal the underlying flow dynamics and characteristics. The spectra analysis based on the Fourier transform and bispectral analysis (Kim & Powers Reference Kim and Powers1979; Nikias & Raghuveer Reference Nikias and Raghuveer1987; Nikias & Mendel Reference Nikias and Mendel1993; Collis et al. Reference Collis, White and Hammond1998) are used for analysing single-point signals. To effectively extract characteristic structures and understand nonlinear interactions across three-dimensional flow fields, spectral proper orthogonal decomposition (SPOD) (Towne et al. Reference Towne, Schmidt and Colonius2018; Yeung & Schmidt Reference Yeung and Schmidt2024), dynamic mode decomposition (DMD) (Rowley et al. Reference Rowley, Mezić, Bagheri, Schlatter and Henningson2009; Schmid Reference Schmid2010; Mezić Reference Mezić2013; Schmid Reference Schmid2022) and bispectral mode decomposition (BMD) (Schmidt Reference Schmidt2020) are used. The selected surface sampling points for the H0100 and H0200 cases are shown in figure 24. The selected surface points can roughly be divided into three groups. The first group is located along the exact attachment line, extending from upstream to downstream of the roughness element, labelled sequentially as
$s_1, s_2, \cdots , s_{14}$
. The second group is located on the side of the roughness element, also extending from upstream to downstream, labelled as
$s_{15}, s_{16}, \cdots , s_{24}$
. The third group, labelled as
$s_{25}, s_{26}, \cdots , s_{32}$
, is also on the side of the roughness element, but is further from the second group of detection points. The first two groups are the same for both cases, while the third group appears only in the H0200 case to track the evolution of the corresponding horseshoe vortex. More details about the statistical sampling strategy, analysis methods and some supplementary analysis are described in the supplementary material to enhance transparency and reproducibility. In the following, the key transition mechanisms and the findings are introduced.
Figure 24. Selected points of surface sampling data for two cases:
$(a)$
for H0100;
$(b)$
for H0200. The corresponding points are sequentially recorded as
$s_1,s_2,\cdots ,s_{32}$
along the
$z$
-axis from upstream to downstream, starting from the attachment line to chordwise downstream. The subscripts
$h_1$
and
$h_2$
are used to distinguish the points in different cases. The average surface skin friction are used to illustrate the relative transition positions.
4.1. Evolutions of the perturbations along the surface
The evolution of perturbation amplitudes of different frequencies, measured along the spanwise (
$z$
-)direction under two operating conditions, is illustrated in figures 25, 26, 27 and 28. In figures 25 and 27, all resolved frequencies are depicted as semi-transparent grey lines, while the evolution of perturbations at several representative frequencies is prominently highlighted. These surface lines along the spanwise
$z$
-direction correspond to the selected group points shown in figure 24. Figures 26 and 28 highlight the features of the selected surface points with longer sampling durations and higher sampling resolutions.
Figure 25. Spanwise evolutions of surface perturbations
$|\hat {\rho }^{\prime }|$
with different frequencies, around the roughness element for the case H0100.
$(a)$
and
$(b)$
Perturbations along specific surface lines, which coincide with the selected first and second groups of points in figure 24
$(a)$
, respectively. The centre of roughness is locate at
$z=40$
.
4.1.1. H0100 case
First, we analyse the evolution of disturbances at the centreline position under the H0100 condition, as illustrated in figure 25
$(a)$
. From this figure, it is evident that as the flow passes over the roughness element, the amplitudes of low-frequency disturbances (e.g. 46.74 kHz) and high-frequency disturbances (e.g. 283.46 kHz, 675.48 kHz and 740.32 kHz) remain largely unaffected. This indicates that the flow field near the roughness element exhibits minimal selective amplification for both low- and high-frequency components.
In contrast, disturbances at intermediate frequencies, such as 93.48 kHz and 141.73 kHz, exhibit notable growth downstream of the roughness element. This suggests that these frequency components are preferentially amplified by the flow field induced around the roughness element. The trend is consistent with the single-point analysis provided in the supplementary material.
As the flow continues to develop further downstream, higher-frequency disturbances are incrementally excited. Eventually, disturbances across all frequency bands saturate, marking the onset of the transition to turbulence. The spanwise evolution of disturbances near the edge of the roughness element shows a similar pattern to that observed at the centreline, as depicted in figure 25
$(b)$
. However, the downstream distance required for disturbances to reach saturation at the edge is significantly greater than at the centreline. This observation aligns well with the features previously discussed.
Figure 26. Spectra
$E_{\rho ^{\prime }}$
of perturbations density
$\rho ^{\prime }$
at the selected points in figure 24
$(a)$
.
$(a)$
and
$(b)$
Two groups of selected points in the H0100 case.
Figure 27. Spanwise evolutions of perturbations
$|\hat {\rho }^{\prime }|$
with different frequencies, around the roughness element for the case H0200.
$(a), (b)$
and
$(c)$
Perturbations along the surface lines of the first, second and third groups of selected points in figure 24
$(b)$
, respectively. The centre of roughness is located at
$z=40$
.
Figure 28. Spectra
$E_{\rho ^{\prime }}$
of perturbations density
$\rho ^{\prime }$
at the selected points in figure 24
$(b)$
.
$(a), (b)$
and
$(c)$
Three groups of selected points in the H0200 case.
These features can be seen more directly and accurately by examining the signal at the selected surface sampling points (shown in figure 24). This is because, for single-point data, it is feasible to use a longer sampling duration and achieve higher frequency resolution, making the spectral characteristics more apparent.
As shown in figure 26, at point
$P_{s_1,h_1}$
, the disturbances exhibit a typical broadband distribution, with a peak around 70 kHz. As the flow approaches the roughness element (from point
$P_{s_1,h_1}$
to
$P_{s_3,h_1}$
), the fundamental amplitude-frequency characteristics remain unchanged, although the amplitude of the perturbations slightly increases. Upon reaching point
$P_{s_4,h_1}$
, the flow enters the leading-edge separation zone induced by the roughness element, where a noticeable rise in the amplitude of low-frequency disturbances occurs. However, due to the strong adverse pressure gradient at the roughness element (at points
$P_{s_5,h_1}$
and
$P_{s_6,h_1}$
), these low-frequency components are partially suppressed, as reflected in the dashed and dotted lines in figure 26
$(a)$
. As the flow passes over the roughness element and enters the wake region (points
$P_{s_7,h_1}$
,
$P_{s_8,h_1}$
and
$P_{s_9,h_1}$
), the perturbation amplitude-frequency characteristics stabilise, nearly coinciding with those at upstream locations (point
$P_{s_6,h_1}$
). Interestingly, at a frequency of approximately 130 kHz, there is a trend of increasing amplitude, which ensures that some of the perturbations are selected by the flow fields. Further downstream (from points
$P_{s_{10},h_1}$
to
$P_{s_{12},h_1}$
), these components grow first and become dominant, accompanied by noticeable growth in low-frequency disturbances. Eventually, at points
$P_{s_{13},h_1}$
and
$P_{s_{14},h_1}$
, perturbations across all frequencies begin to grow, gradually exhibiting spectral characteristics indicative of turbulent flow. The evolution of disturbances at the edge of the roughness element, as shown in figure 26
$(b)$
, exhibits similar characteristics. However, the corresponding disturbance amplitudes remain consistently lower, which indirectly explains why the transition does not initially occur at the edges.
4.1.2. H0200 case
Next, we examine the evolution characteristics of disturbances under the H0200 condition, which is characterised by a higher roughness element, as shown in figure 27. As the flow progresses towards the roughness element, disturbances across nearly all frequencies are amplified, particularly in the region where
$z= [40,42 ]$
. Further downstream, low-frequency components (e.g. 10.55 kHz, 21.11 kHz and 42.22 kHz) rapidly grow to saturation. Subsequently, high-frequency components begin to grow only after the low-frequency disturbances have remained saturated for some time. Eventually, disturbances at all frequencies reach a nearly saturated state, signalling the transition to turbulence.
A noteworthy phenomenon is the variation in dominant low-frequency disturbances at different positions. At the centreline, as shown in figure 27
$(a)$
, the 21.11 kHz disturbance dominates. In contrast, at the edge of the roughness element and in the trailing vortex region, as depicted in figures 27
$(b)$
and 27
$(c)$
, the 10.55 kHz disturbance becomes the dominant low-frequency component. This variation highlights the three-dimensional nature of the disturbance characteristic structure and we will clarify this in § 4.3.1.
Another critical observation concerns the evolution of disturbances in the horseshoe vortex region and the edge of the roughness element, as shown in figures 27
$(b)$
and 27
$(c)$
. In the horseshoe vortex region (notably in the
$z= [39,40 ]$
interval), low-frequency disturbances remain largely unaffected, exhibiting nearly linear growth in amplitude. In contrast, high-frequency components experience significant modulation, reflected in the synchronous behaviour of their amplitudes. Interestingly, in the edge region of the roughness element and the trailing vortex region, the low-frequency disturbances saturate earlier than those at the centreline. This phenomenon aligns well with the trends observed in the mean flow analysis discussed previously.
Similar to the previous analysis, we further validate the aforementioned process by using signals with extended sampling durations and higher resolutions obtained at the surface characteristic points. Observing the evolution patterns along the centreline (figure 28
$a$
), the absence of upstream disturbances results in a significantly quieter flow field ahead of the roughness element compared with the H0100 scenario. Overall, disturbance amplitude levels are relatively low, but increase as the flow approaches the roughness element (from points
$P_{s_1, h_2}$
to
$P_{s_4, h_2}$
).
As the flow reaches the roughness element, disturbances are suppressed due to the adverse pressure gradient, particularly at point
$P_{s_5, h_2}$
. However, upon entering the region behind the roughness element (point
$P_{s_6, h_2}$
), disturbances in the low-frequency range (around 20 kHz) exhibit a significant increase, with amplitudes rising by nearly three orders of magnitude. Further downstream, high-frequency disturbances begin to grow, and by point
$P_{s_8, h_2}$
, disturbances at frequencies around 450 kHz grow preferentially. With continued downstream progression (from points
$P_{s_{9}, h_2}$
to
$P_{s_{14}, h_2}$
), disturbances across various frequency ranges reach saturation and the flow transitions into the turbulent phase.
In the H0200 scenario, the increased height of the roughness element results in a more pronounced horseshoe vortex (figure 22). We further investigate the evolution of disturbances in its wake region, corresponding to the third group of points (figure 28
$c$
). Influenced by the horseshoe vortex, the disturbance spectrum at point
$P_{s_{25},h_2}$
shows considerable amplitude in the low-frequency region. As the flow progresses downstream, the disturbance amplitude increases and at point
$P_{s_{27},h_2}$
, disturbances at approximately 1000 kHz exhibit marked preferential growth. This component becomes dominant in the high-frequency range, reaching saturation first at point
$P_{s_{29},h_2}$
.
By comparing the amplitudes at the same spanwise location in the centre wake region (point
$P_{s_{9}, h_2}$
, figure 28
$a$
), it is evident that high-frequency disturbances are triggered earlier and reach saturation preferentially, leading to the transition in the edge region occurring prior to that along the centreline. This is also clearly reflected in the transition visualisation, as discussed previously.
Finally, we examine the second set of points located at the edge of the roughness element (figure 28
$b$
). The basic spectral characteristics are similar to those of the third set, although the growth along the mainstream direction and the corresponding amplitudes are smaller. This observation partially explains why, in the H0200 scenario, transition is preferentially observed at the sides and along the centreline.
4.2. Nonlinear interactions of the transitional flow fields
To further understand the potential nonlinear relationships among the flow fields, BMD methods have been used to analyse the possible interactions in a global point of view. As a complementary effort, traditional bispectral analysis was also conducted and the detailed results can be found in the supplementary material. The magnitudes of mode bispectrum
$\log {|\lambda _1|}$
for both cases are shown in figures 29 and 30. For the H0100 case, the primary nonlinear effect occurs between the unsteady disturbances and the mean-flow field. This is evident as the intervals with significant mode bispectrum amplitude are located along the
$\omega _2^*=0$
horizontal axis in the figure. This nonlinear effect indicates that the mean-flow field is nonlinear coupled with the dominant unsteady disturbances. The dominant disturbance frequencies extracted using the BMD method correspond to the characteristic frequencies identified by the three-dimensional SPOD mode analysis provided in detail in the supplementary material. Additionally, there exist relatively weaker nonlinear effects among these dominant structures, as also presented in the figure. It is important to note the significant differences from the single-point bispectrum analysis shown in supplementary figure 4. This discrepancy arises because the projections of nonlinear interactions among three-dimensional spatial structures at a single point are often inaccurate.
Figure 29. Magnitude of mode bispectrum
$\log {|\lambda _1|}$
of spanwise velocity perturbation
$w^{\prime }$
for the H0100 case.
Figure 30. Magnitude of mode bispectrum
$\log {|\lambda _1|}$
of spanwise velocity perturbation
$w^{\prime }$
for the H0200 case:
$(a)$
low-frequency region;
$(b)$
high-frequency region.
For the H0200 case, the BMD analysis reveals notably different structures. The primary nonlinear interactions concentrate in the low-frequency range, predominantly featuring the self-interaction of the 10 kHz disturbances, peaking at the point
$(\omega _1^*, \omega _2^*) \approx (10\,\mathrm{ kHz}, 10 \,\mathrm{kHz})$
in figure 30
$(a)$
. Additionally, there are significant nonlinear interactions both between 20 kHz and 10 kHz disturbances, and within themselves. However, unlike the H0100 case, the interactions with the mean-flow field are weaker. In other regions, nonlinear interactions span a broad frequency range. These characteristics partially explain the earlier transition observed in the H0200 scenario. In figure 30
$(b)$
, the mode bispectrum displays prominent distributions at several points, indicating nonlinear interactions between the low-frequency 10 kHz and 20 kHz disturbances and higher frequency disturbances. Coupled with modal energy(eigenvalues) and the modal structures from SPOD analyses provided in the supplementary material, these high-frequency disturbances are likely the result of secondary instabilities occurring after the low-frequency unsteady streak structures have fully developed.
4.3. Clarification of the absolute instability for H0200 case
In previous analyses, the primary flow mechanisms for two distinct scenarios were successfully identified. However, for the H0200 case, two fundamental questions remain unanswered.
-
1. What causes the inconsistencies among the dominant low-frequency modes extracted from spectral analyses conducted at different spatial locations? Which frequency truly represents the global dominance in the system?
-
2. What is the origin of the self-sustained low-frequency disturbances observed in the flow?
To answer the first question, we propose applying SPOD to the three-dimensional transitional flow field. This technique allows us to identify the characteristic frequency associated with the most dominant global three-dimensional structure directly. Addressing the second question involves considering two potential mechanisms that could drive these disturbances.
-
1. Local absolute instabilities, leading to the spontaneous oscillation of flow separation bubbles.
-
2. Feedback loops, resulting from flow separation dynamics and their downstream influence.
In this study, we hypothesise that the first mechanism – local absolute instabilities – is the more likely explanation. This hypothesis is supported by observations of spanwise upstream propagation in the dominant SPOD modes for the H0200 case, as revealed in animations of the flow field. Specifically, we argue that the local absolute instabilities within the separation bubbles forming downstream of roughness elements lead to the observed self-sustained low-frequency disturbances.
4.3.1. Clarification of the dominant frequency
Figure 31. SPOD modes of the spanwise velocity perturbation
$w^{\prime }$
around the three-dimensional transition flow field for the H0200 case:
$(a)$
leading four SPOD spectra (
$\lambda _{f_{k,1}}, \lambda _{f_{k,2}}, \lambda _{f_{k,3}}$
and
$\lambda _{f_{k,4}}$
);
$(b) {-} (e)$
real parts of the spanwise perturbation velocity fields
${\hat{w}}^{\prime}$
of the first SPOD modes
$\lambda _{f_{k,1}}$
at the four indicated frequencies in panel
$(a)$
. These modes are depicted by isocontours of the spanwise velocity perturbation
${\hat{w}}^{\prime }$
. The contour levels depict
$\pm$
10 % of the mode’s maximum spanwise velocity perturbation. The corresponding light purple isosurface in panel
$(b)$
represents the separated bubble in the corresponding average flow field, approximated using the isosurface of
$\overline {w} = -0.0001$
.
The SPOD results are presented in figure 31. The spectra in figure 31
$(a)$
identify a dominant low-frequency component at approximately
$11$
kHz, consistent with the previously observed frequency of
$10.55$
kHz. The corresponding three-dimensional structure, shown in figure 31
$(b)$
, reveals streak-like features located near the separation bubble downstream of the roughness element. In addition, a high-frequency structure at
$486$
kHz is identified, which aligns well with the BMD mode bispectrum shown in figure 30
$(b)$
and the spectra reported previously. The spatial distribution of this mode, illustrated in figure 31
$(d)$
, overlaps with the low-frequency streak-like structure described earlier. Collectively, these findings strongly suggest that the identified high-frequency mode is most likely driven by the secondary instability of the low-frequency streaks.
Two additional tonal signals at
$61$
kHz and
$156$
kHz are identified through SPOD analysis. The corresponding mode structures are not confined solely to the central region behind the roughness element, but extend into the lateral regions. This broader spatial distribution may explain the significantly larger transition region observed in this case. As discussed earlier, transition onset occurs earlier at the edge of the roughness element; however, the spatial structures of the tonal signals do not exhibit a similar trend. This discrepancy arises because the localised flow features near the roughness edge are not globally dominant in the three-dimensional transitional flow field.
In the H0200 scenario, three-dimensional SPOD analysis identifies a dominant low-frequency mode at approximately 11 kHz. In contrast, spectral analyses along the centreline reveal a dominant frequency of approximately 21 kHz – roughly double the frequency identified in the three-dimensional analysis. This discrepancy likely arises from the inherently three-dimensional nature of the unsteady flow behaviour. To further understand this discrepancy, the original physical signal can be examined, as illustrated in supplementary figure 3. Observations of the chordwise perturbation velocity component at a point perpendicular to the attachment line reveal a substantially longer oscillation period compared with other physical variables, nearly twice as long. This observation aligns with the 11 kHz dominant mode identified in the three-dimensional SPOD analysis and suggests that the chordwise low-frequency unsteady evolution mode governs the dynamics of the three-dimensional transitional flow field, rather than unsteady evolution occurring along the mainstream direction. The relationship between the two dominant frequencies, 11 kHz and 21 kHz, is attributed to harmonic interactions arising from the intrinsic coupling of the unsteady behaviours in the flow. These findings highlight the necessity of integrating both high-dimensional spectral analyses and localised investigations to comprehensively characterise the dominant structures and dynamics in transitional flow fields.
The length and velocity scales are introduced to give an estimation of the dimensionless frequency. The Strouhal number (
$St$
) can be expressed in terms of the length and velocity scales as follows:
\begin{equation} St = \frac {\omega ^* \cdot L_{r_2}^*}{U_{r_2}}. \end{equation}
The length scale
$L_{r_2}^*$
is defined directly as the diameter of the roughness element. The characteristic velocity
$U_{r_2}^*$
is defined as
\begin{equation} U_{r_2}^* = \frac {1}{\nu ^*_w}\int _0^{k_h^*}{W^{*2}}\mathrm {d}h, \end{equation}
to consider the local kinetic energy contributions while accounting for viscous effects, providing a physically meaningful measure of the flow’s energetic state near the wall. The spanwise velocity
$W^*$
and kinematic viscosity
$\nu ^*_w$
at the wall surface are obtained at the attachment line in a smooth configuration without the roughness element, and
$k_h^*$
represents the height of the roughness element. This definition is specifically designed to decouple the local influence of roughness from free stream conditions, making it ideal for analysing roughness-induced oscillation. Based on these scales, the Strouhal number for the dominant low-frequency perturbation (
$10.55$
kHz) is found to be approximately 0.75 for the present H0200 case.
4.3.2. Origin of the self-sustained perturbations
As mentioned before, the local absolute instabilities are hypothesised to be responsible for the self-sustained perturbations. To rigorously investigate this hypothesis, DMD is used to distinguish the respective modes and quantify their growth rates directly linked to such phenomena. Absolute instability implies spontaneously growing disturbances, characterised by dominant DMD modes concentrated in specific regions with significantly high growth rates, indicating robust local growth. In contrast, feedback-loop-induced global instabilities might exhibit coupled modes extending across broader spatial regions. For the feedback loop mechanism, the reliance on flow feedback often leads to a balanced energy transfer, resulting in slower mode growth, as indicated by relatively small growth rates. Under certain critical conditions, when the system reaches a state of sustained oscillation or coherent behaviour, the growth rate may approach zero.
The results of DMD are shown in figure 32. Figures 32
$(d)$
and 32
$(e)$
illustrate the variation in growth rates of DMD modes with frequency, identifying three unstable modes, also shown in figures 32
$(a)$
and 32
$(b)$
. One of these is the mean-flow mode, characterised by
$\text {Im}(\mu _i) = 0$
, or a perturbation frequency of
$\omega ^*=0$
kHz, referred to as Mode 0. The other two unstable modes are named Mode 1 and Mode 2, corresponding to frequencies of 11 kHz and 131 kHz, respectively, and both exhibit substantial growth rates. Figures 32
$(f)$
and 32
$(g)$
depict the three-dimensional spatial distribution of these unsteady modes. Mode 1, at 11 kHz, concentrates near the separation bubble at the trailing edge of the roughness element. This mode closely resembles the previous SPOD analysis results, displaying a streak-like structure, with a frequency matching those identified as dominant in SPOD. Unlike the localised distribution of Mode 1, Mode 2 (131 kHz) is spatially broader, extending downstream of the roughness element with a distinct three-stranded distribution. This characteristic aligns with previously general features of the transition for the H0200 case. Therefore, building on previous analyses, we conclude that under the H0200 condition, the low-frequency disturbances (approximately 10 kHz) originate from local absolute instabilities in the separation bubble region behind the roughness elements. This finding also explains why, under H0200 conditions, transition can be effectively triggered solely by inherent numerical disturbances in the flow field, without the need for additional perturbations.
Figure 32. DMD modes for spanwise velocity perturbations
$w^{\prime }$
:
$(a)$
eigenvalues
$\mu _i$
(black circles) for DMD. The red points stand for selected modes;
$(b)$
,
$(c)$
eigenvalues (black circles) near selected modes (red points) in panel
$(a)$
with blue dashed lines representing neutral lines.
$(d)$
Growth rate of DMD modes relative to frequency;
$(e)$
growth rate near selected modes, with blue lines indicating zero-growth neutral lines.
$(f)$
,
$(g)$
Real parts of selected unsteady DMD modes, shown as isocontours of spanwise velocity perturbation, at
$\pm$
10 % of the maximum spanwise velocity perturbation. In panel
$(f)$
, the light purple isosurface denotes the separated bubble in the average flow field, approximated as
$\overline {w} = -0.0001$
.
4.4. Effects of the roughness height on boundary layer transition
Based on the previous results, we can summarise and discuss the influence of roughness element height on the transition characteristics of the present flows.
First, when the height of the roughness elements is below the critical roughness height, their impact on the flow field is limited. In this regime, the roughness elements merely excite steady (time-independent) flow structures and are unable to directly induce unstable disturbances. However, the three-dimensional flow field induced by the roughness element acts as a selector and amplifier of unsteady disturbances from the upstream flow, preferentially selecting certain disturbances and promoting their growth. After the saturation of these disturbances, high-frequency perturbations are subsequently excited, eventually leading to the final transition.
When the roughness element height exceeds the critical value, local regions of absolute instability emerge within the three-dimensional flow field induced by the roughness elements. These regions are located in the downstream three-dimensional separation bubble. The absolute instability amplifies disturbances continuously, leading to a spontaneous transition process. The characteristic frequency of the disturbances associated with this absolute instability directly corresponds to the low-frequency movement in the chordwise direction of the separation bubble, rather than the streamwise direction. Subsequently, the secondary instability of the streaks induced by these low-frequency disturbances further drives the transition in the wake region downstream of the roughness elements.
Moreover, as the height of the roughness elements increases, stronger horseshoe vortices are induced, which rapidly generate saturated streak structures in the downstream region. These saturated streak structures, resembling saturated cross-flow vortices, undergo high-frequency secondary instability that leads to transition.
5. Conclusions and discussions
5.1. Conclusions
In summary, this study presents simulations of the complete transition process from laminar flow to turbulence in a hypersonic three-dimensional swept leading-edge boundary layer over an experimental configuration. Based on the results from linear stability theory and the experimental observations, the subcritical state of the possible transition along the attachment line is confirmed. Two roughnesses based on the experimental tests are modelled and designed to trigger transitions. The main findings of this study can be summarised as the following points.
-
1. In the flow over a swept blunt body discussed in this paper, even without the assumption of infinite sweep, if the incoming boundary layer reaches an asymptotic state (for laminar flow), the subsequent flow state also satisfies the infinite sweep assumption – the boundary layer is homogeneous along the swept direction.
-
2. The two different heights of roughness elements in the configuration studied in this paper result in completely different transition characteristics. For a lower-height roughness element, the element alone cannot directly induce the corresponding boundary layer transition. Certain random perturbations need to be introduced during the simulation. For a higher-height roughness element, it can directly induce boundary layer transition by itself without the need for additional perturbations. The flow near the roughness element resembles that of a flat-plate boundary layer, where vortex structures are triggered in its vicinity. In the transition phenomena induced by external perturbations and lower-height roughness elements, the transition mainly occurs directly downstream of the wake of the roughness element, but the wake vortices induced by the roughness element do not directly destabilise and lead to the final transition. In the transition simulation with higher-height roughness element, the horseshoe vortices generated by the roughness element form corresponding streaks, which preferentially destabilise and lead to the transition to turbulence, while the wake vortices directly behind the roughness element destabilise and transition further downstream.
-
3. Using single-point, surface and full three-dimensional statistical data combined with spectral, bispectral, SPOD, DMD and BMD analyses, the complete transition processes are elucidated. Fora small roughness element, the wake flow behind the roughness element acts as a perturbation selector and amplifier by selecting and amplifying incoming perturbations from upstream, particularly those with frequencies between 100 kHz and 200 kHz, leading to the final transition. In contrast, fora larger roughness element, a low-frequency local absolute instability exists in the separation bubble induced by the roughness element. This low-frequency perturbation, approximately 10 kHz, generates corresponding low-frequency streaks. The high-frequency secondary instability of these low-frequency streaks is the primary mechanism driving transition.
5.2. Further discussions
In this study, we demonstrate that combining traditional local statistical analysis with high-dimensional (two-dimensional (2-D)/3-D) statistical approaches is crucial for a comprehensive understanding of transitional flow fields. High-dimensional analysis provides direct insights into the globally significant structures in physical space, offering a holistic perspective on the flow dynamics. However, it is important to note that certain local features may be masked or smoothed out in such global analyses. Conversely, traditional single-point statistical analysis captures all signal characteristics at specific locations, but it neglects the global context and the interconnections of these features across the flow field. Therefore, we argue that integrating these two approaches is essential for effectively analysing complex three-dimensional flow fields.
In this study, under the conditions, even when the height of the roughness element reaches half of the local boundary layer thickness, significantly exceeding the critical height in normal flat-plate boundary layers (Bernardini et al. Reference Bernardini, Pirozzoli and Orlandi2012, Reference Bernardini, Pirozzoli, Orlandi and Lele2014; Estruch-Samper et al. Reference Estruch-Samper, Hillier, Vanstone and Ganapathisubramani2017), vortical structures formed downstream of the roughness element rapidly decay and do not lead to transition in the absence of additional disturbances. However, the presence of small disturbances can potentially induce transition in the wake. This phenomenon is not explored in detail in this study, but possible explanations include the significant differences between attachment-line boundary layer flow and conventional flat-plate boundary layer flow. In the external region of a flat-plate boundary layer (where the major flow direction is along the
$z$
-axis), both the
$u$
velocity and its gradient
$\partial u/\partial x$
are approximately zero. In contrast, within the attachment-line boundary layer, although the
$u$
velocity is small at the region, the gradient
$\partial u/\partial x$
is relatively large. This tendency or acceleration from the centre towards the chordwise direction might inhibit the concentration of disturbances within the wake of the roughness element. Additionally, in the attachment-line boundary layer, the flow exhibits a wall-normal impingement
$(v\lt 0)$
, whereas, in normal flat-plate boundary layers, a positive wall-normal velocity
$(v\gt 0)$
is often observed. This results in a thinner boundary layer and suppression of disturbance development. These differences may collectively contribute to the variation in transition behaviour.
Another notable point is that under the H0200 condition, the transition near the edges of the roughness element occurs earlier than in the wake region. In fact, this phenomenon is also observed in roughness-induced transitions of hypersonic flat-plate boundary layers (Lefieux et al. Reference Lefieux, Garnier and Sandham2019). These findings indicate significant differences in transition behaviour between the edge and wake regions induced by the roughness element. Therefore, when addressing such problems in the future, particular attention should be paid to potential transitions occurring near the edge regions.
Acknowledgements.
We acknowledge Yancheng MetaStone Tech. Co. for providing us with the computational resources required by this work. Useful discussions with Professor Qibing Li of Tsinghua University and Professor Jie Ren of Beijing Institute of Technology are gratefully acknowledged.
Funding.
This work received support from the NSFC Grants 12388101, 12202242 and 12172195. The authors are also grateful for the support from the National Key Research and Development Plan of China through project no. 2019YFA0405201, the National Key Project GJXM92579 and the Grant 20231001 for the Supercomputing Center in Yancheng.
Declaration of interests.
The authors report no conflict of interest.
Data availability statement.
The full data set of the simulations is of the order of 80 thousands of gigabytes. By contacting the authors, a smaller subset can be made available.
Author contributions.
Youcheng Xi: funding acquisition, conceptualisation, data curation, formal analysis, coding, investigation, methodology, validation, writing-original draft. Bowen Yan and Guangwen Yang: funding acquisition, software and optimisation on high performance computing, computational resources. Xinguo Sha and Dehua Zhu: model designing, wind tunnel testing, experimental investigation. Song Fu: funding acquisition, supervision, resources, writing-review and editing.
Supplementary materials and movies.
Supplementary materials and movies are available at https://doi.org/10.1017/jfm.2025.269
Appendix A. Comparison of the laminar states of the present configuration and a swept circular cylinder
The laminar flow for a swept circular cylinder is presented for comparison. The basic comparison of the surface shapes are shown in figure 33 and the laminar profiles at the exact attachment line are shown in figure 34. Due to the relatively large equivalent radius of curvature, the boundary layer formed on the present configuration is slightly thicker than that on a circular cylinder. However, the fundamental distribution patterns of the physical quantities remain consistent.
Figure 33. The basic shapes of the leading edge for circular cylinder and the present configuration.
Figure 34. Profiles of major variables at the
$x=0$
plane.
$h$
stands for the distance away from the wall surface.
Open access
