3.1. Global flow pattern of double-cone flows

The schlieren technique effectively visualises shock wave structures and flow separation according to the density changes in hypersonic flow fields. Figure 5 displays schlieren images captured experimentally at a sampling rate of 75 kfps, along with corresponding numerical schlieren results for the three scenarios. In the numerical simulations, the transition locations were fixed at $x_{t}/L$ = 0.86, 0.85 and 0.59 for cases RN0, RN4 and RN8, respectively (details on the estimation of $x_{t}$ can be found in Appendix A), yielding the closest alignment with experimental data. The schlieren results depicted in figure 5 distinctly reveal the presence of a separation bubble at the cone junction, alongside the complex shock wave system generated by the flow interaction. Both features are well reproduced by the numerical simulation, closely matching the experimental phenomena.

Figure 5.

The experimental and numerical schlieren images of cases ( $a$ ) RN0, ( $b$ ) RN4 and ( $c$ ) RN8, with flows structures annotated, including finescale structures (FSS), nose-induced oblique shock wave (NOSW), separation shock waves (SSW), transmitted oblique shock wave on the first cone (OSW1), bow shock wave (BSW), transmitted shock wave (TSW), reattachment shock wave (RESW), Mach stem (MS), sonic line (SOL), triple point (TP) and a relatively weak oblique shock wave on the second cone (OSW2). The white markers indicate the positions of the four Kulite sensors.

The schlieren results indicate that the separation point gradually moves upstream as the nose bluntness increases. The separation bubbles in cases RN0 and RN4 are very similar, exhibiting typical laminar boundary-layer separation morphology on the first cone. When the nose radius increases to 8 mm, the thick entropy layer downstream of the nose remains unswallowed by the boundary layer within a short streamwise distance, and its influence persists. The interaction between the shock wave and the boundary layer, along with the entropy layer, leads to a significant morphological alteration of the separation bubble compared with cases RN0 and RN4. Notably, case RN8 exhibits the emergence of finescale, unsteady vortical structures on the first cone, as shown by the two snapshots in figure 5( $c$ ). Nevertheless, a common feature across all three cases is a prominent Type V SSI. In each case, the separation shock wave first interacts with the nose-induced shock wave in a Type VI interaction, generating the transmitted oblique shock wave OSW1. The OSW1 subsequently intersects the detached bow shock ahead of the second cone, forming a distinct TSW and a TP (TP1) characteristic of a Type V pattern. In the flow field above TP1, the flow is subsonic immediately downstream of the bow shock, with the corresponding region enclosed by the sonic line (the purple solid line) indicated in the numerical schlieren images. Conversely, the airflow below TP1 remains supersonic as it passes through the TSW. A detailed inspection reveals two types of shock wave systems below TP1 in the three cases. For cases RN0 and RN4, the separation bubble is relatively small, and the compression waves rapidly coalesce into an oblique shock wave during the reattachment process. This shock wave intersects with TSW, resulting in a Mach shock wave reflection typified by a discernible Mach stem and downstream slip lines, observable in the close-up images of both experimental and numerical schlieren results in figure 5. In contrast, the extended separation region in case RN8 inhibits full convergence of compression waves, yielding a weaker reattachment shock wave that engages in regular shock wave reflection with the TSW. Consequently, due to the two different wave structures below TP1, the downstream supersonic jet exhibits two patterns: RN0 and RN4 develop a three-layer jet stratified by slip lines emanating from TPs TP2 and TP3, whereas RN8 forms a single-layer structure. It is worth noting that in the experimental schlieren images, both the three-layer and single-layer structures of the supersonic jet maintain only a short streamwise distance before large-scale vortex structures emerge and disrupt the stratified flow. However, these large-scale vortical structures within the supersonic jet are not successfully reproduced in steady RANS simulations due to their inherent limitations in resolving unsteady vortical structures.

In fact, the detached bow shock wave undergoes progressive attenuation along the streamwise direction, transitioning into an oblique shock wave while maintaining supersonic flow conditions downstream. In cases RN0 and RN4, schlieren images indicate that the supersonic airflow traversing the middle region of the second cone decelerates again through a relatively weak oblique shock wave (designated OSW 2 in figure 5), aligning the airflow gradually parallel to the cone surface. This weak shock wave intersects with the outer oblique shock wave, resulting in a Type VI SSI, along with another TP (TP4) and a discernible slip line. In contrast, although the angle of the outer oblique shock wave in case RN8 shows an increasing trend near the tail region of the second cone, no typical Type VI SSI is observed. It can be inferred that this absence may be attributed to complex wave interactions wherein the outer oblique shock wave interacts with distributed weak compression waves. These compression waves, while physically present, possess insufficient density gradients for clear visualisation in schlieren images.

Figure 6.

Enlarged Mach number contours near the cone junction and the friction coefficient distributions for cases ( $a{,}b$ ) RN0, ( $c{,}d$ ) RN4 and ( $e{,}f$ ) RN8.

Regarding the separation bubble, a comprehensive analysis is carried out based on the numerical results. Figure 6( $a$ , $c$ , $e$ ) present the enlarged Mach number contours delineating flow structures near the separation zones. The streamline within the separation bubble is depicted through the yellow arrowed lines. Figure 6( $b$ , $d$ , $f$ ) display the friction coefficient ( $C_{\kern-1pt f}$ ) distributions; herein, $C_{\kern-1pt f}$ = $\tau _w$ /0.5 ${\rho _0}{u_0^2}$ , where $\tau _w$ , $\rho _0$ , $u_0$ denote the wall shear stress, free stream density and free stream velocity, respectively. The results indicate that for cases RN0 and RN4, the reflected shock wave originating from TP3 impinges on the boundary layer downstream of the reattachment of the separation bubble at the cone junction, leading to another separation bubble. The streamlines within the separation bubbles demonstrate the spatial segregation of these two recirculation regions, as further evidenced by the $C_{\kern-1pt f}$ distributions. It is worth noting that a secondary separation occurs within the first separation bubble, as shown in both the Mach number contours and $C_{\kern-1pt f}$ distributions. Conversely, for case RN8, the adverse pressure gradient caused by the geometric deflection at the cone junction combines with that induced by the impingement of TSW, resulting in a relatively large-scale coupled separation bubble without any secondary separation occurring beneath it. The axial locations of the separation and reattachment points for all three cases are listed in table 2. Results regarding the secondary separation beneath the first separation bubble are not included.

Table 2.

Axial locations of the separation and reattachment points. The symbols $S$ and $R$ in the subscripts denote the separation and reattachment points, respectively. The numbers 1 and 2 correspond to the first and second separation bubbles, respectively.

The PSP-measured pressure contours within a 120 $^{\circ }$ azimuthal region are shown in figure 7( $a$ , $c$ , $e$ ). These results demonstrate a clearly axisymmetric flow, in contrast to the significant three-dimensional structures reported in the numerical work of Hao et al. (Reference Hao, Fan, Cao and Wen2022) when the flow is fully laminar. Notably, the observed minor curvature of high-pressure regions along the azimuthal direction is attributed to the optical image distortion, and the actual cone junction is marked by the red dashed lines. The streamwise pressure coefficient distributions ( $C_p = {p_w}/0.5{\rho _0}{u_0^2}$ , where $p_w$ is the wall pressure), derived from numerical simulations, pressure sensors and PSP results along the central meridian line, are presented in figure 7( $b$ , $d$ , $f$ ), indicating good agreements between the experimental and numerical results. It is noteworthy that the PSP data along the model centreline were prioritised for optical distortion correction, yielding a streamwise displacement $\Delta x \lt 0.4$ mm. This ensures a negligible impact on the measured pressure peaks, as shown in figure 7( $b$ ). It is clear that the peak pressures in cases RN0 and RN4 are significantly higher than that in case RN8. Combining the flow structures shown in figures 6 and 7, the pressure rises before reaching the peak value in cases RN0 and RN4, and can be divided into three stages. The first is the pressure rise induced by the airflow deflection at the separation point of the first bubble; the second pressure-rise stage occurs during the reattachment process of the first separation bubble; and the third stage is induced by the impingement of TSW on the boundary layer on the second cone. The latter two stages exhibit a steeper pressure rise compared with the first. In case RN8, the coupling of the two separation regions causes the latter two pressure rises to merge into a single gradual increase, primarily driven by a series of compression waves. Furthermore, the shock wave emitted from the TP TP2 in cases RN0 and RN4, along with the reattachment shock wave in case RN8, impinge on the shear layer downstream of TP TP1. These interactions result in strong reflected expansion waves, which subsequently impinge on the cone surface, causing a rapid pressure decrease. Additionally, these expansion waves and the shock waves within the supersonic jet undergo repeated reflections between the shear layer and the cone surface, leading to fluctuations in pressure distributions on the second cone.

Figure 7.

Pressure contours and pressure distributions along the meridian lines for cases ( $a{,}b$ ) RN0, ( $c{,}d$ ) RN4 and ( $e{,}f$ ) RN8.

Figure 8.

The comparison of ( $a$ ) pressure coefficient $C_p$ and ( $b$ ) friction coefficient $C_{\kern-1pt f}$ distributions. The black dashed line indicates the location of the cone junction. The dots denote the upstream conditions prior to the interaction region.

Figure 8 presents the comparison of $C_p$ and $C_{\kern-1pt f}$ distributions for the three cases. While all cases exhibit a characteristic pressure plateau downstream of the separation point, distinct plateau magnitudes are observed. As the nose bluntness increases, the pressure plateau magnitude $p_p$ gradually decreases. As shown in the $C_p$ distributions, an over-expansion effect occurs downstream of the blunt noses in cases RN4 and RN8, leading to lower pressures compared with the sharp-nose case; correspondingly, the friction coefficients in these blunt cases are also reduced relative to case RN0. Comparing cases RN0 and RN4 reveals that although the wall pressure upstream of separation recovers to the RN0 level in case RN4, the friction coefficient remains significantly lower. According to the free interaction theory (Chapman, Kuehn & Larson Reference Chapman, Kuehn and Larson1957), $p_p$ is related to ${p_0}{C_{\kern-1pt f,0}^{1/2}}$ (where $p_0$ and $C_{\kern-1pt f,0}$ represent the wall pressure and friction coefficient upstream of the separation region). The lower $C_{\kern-1pt f,0}$ in case RN4 directly explains its reduced $p_p$ relative to case RN0. In case RN8, insufficient streamwise distance prevents the wall pressure and friction coefficient from fully recovering to RN0 levels upstream of the separation point. These lower upstream values consequently contribute to a further reduction in the $p_p$ compared with both cases RN0 and RN4. Moreover, the plateau terminations in cases RN0 and RN4 coincide with the cone junction (marked by the black dashed line), and the pressure rise downstream of this point is caused by the reattachment process. However, a notable deviation occurs in case RN8: both numerical simulation and PSP data reveal a gradual pressure rise in the region downstream of the plateau (at approximately $x/L$ = 0.7, marked by the green dashed line) and before the cone junction. This subtle pressure increase precedes the sharper rise associated with reattachment. Considering that free interaction theory predicts higher plateau pressures for shock interactions involving transitional or turbulent boundary layers compared with purely laminar interactions (Babinsky & Harvey Reference Babinsky and Harvey2011), one possible explanation for this slow pressure rise in case RN8 is the onset of transition to turbulence within the shear layer overlying the separation bubble on the first cone. This developing turbulence would enhance momentum transport, gradually increasing the pressure upstream of the reattachment.

3.2. Dynamic analysis of double-cone flows

In the previous section, the global flow patterns for the three configurations were introduced. Notably, examination of the schlieren sequences revealed the presence of shock-wave oscillations across all cases (refer to the Supplementary materials for the movies are available at https://doi.org/10.1017/jfm.2026.11443). However, the intensity of these oscillations is significantly weaker than the large-scale pulsations reported in earlier experimental and numerical studies (Jagadeesh et al. Reference Jagadeesh, Hashimoto, Naitou, Sun and Takayama2003; Hornung et al. Reference Hornung, Gollan and Jacobs2021). To delve deeper into the unsteady characteristics of these cases, we conduct a comprehensive analysis that integrates data obtained from both schlieren visualisation and pressure measurements in this section.

Unlike the two-dimensional geometries, such as compression corners and double-wedge configurations, the axisymmetric double-cone models employed in the current study offer a distinct advantage when utilising schlieren techniques. The axisymmetric features effectively minimise spatial integration effects along the spanwise direction, thereby enhancing the fidelity of the captured flow dynamics. This set-up enables spectral analysis via the SPOD method. This technique has been widely applied to investigate unsteadiness in supersonic and hypersonic flows, including applications to second-mode instabilities (Butler & Laurence Reference Butler and Laurence2021), shear-layer structures (Benitez et al. Reference Benitez, Borg, Scholten, Paredes, McDaniel and Jewell2023, Reference Benitez, Borg and Hill2024) and shock wave-separation bubble dynamics (Hoffman et al. Reference Hoffman, Rodriguez, Cottier, Combs, Bathel, Weisberger, Jones, Schmisseur and Kreth2022; Li et al. Reference Li, Li, Ye and Zhang2024; Posladek et al. Reference Posladek, Andrade, Combs and Glasby2024). The SPOD method decomposes the flow field into a hierarchy of orthogonal modes, each oscillating at distinct frequencies, which collectively capture the coherent spatiotemporal evolution of dominant flow structures. Motivated by the identification of relatively high-frequency structures (approaching 100 kHz) in the supersonic jet within a Mach 7 double-cone flow, as reported in the previous experimental study by Posladek et al. (Reference Posladek, Andrade, Combs and Glasby2024), this study carried out the spectral analysis based on the schlieren sequences at sampling rates of both 75 kfps and 300 kfps. The schlieren visualisation at the higher sampling rate of 300 kfps primarily focuses on the flow field around the second cone, targeting the potential capture of these high-frequency phenomena.

Regarding the SPOD analysis applied to schlieren sequences acquired at 75 kfps, the Welch method was employed with a Hamming window of length 512 points to mitigate spectral leakage. The schlieren sequences were segmented into 14 overlapping blocks using a 75 $\,\%$ overlap ratio. This set-up yielded a frequency resolution of 146 Hz. Normalising this resolution using the characteristic length $L$ and free stream velocity $u_0$ results in a Strouhal number resolution of $St_L$ = 0.008. Consistent with typical practice for identifying dominant flow features, the analysis focuses exclusively on the leading SPOD mode due to its highest energy content. Figure 9( $a$ ) present the modal energy distributions of the leading SPOD modes for all three cases. The results indicate that for cases RN0 and RN4, the modal energies decrease with increasing frequency, reflecting broadband low-frequency characteristics. Interestingly, when the nose radius increases to 8 mm, a significant peak emerges at around $St_L$ = 0.12 ( $f$ = 2.197 kHz) in the leading mode energy, signifying the presence of dominant low-frequency flow unsteadiness in this case. Figure 9( $b$ , $c$ , $d$ ) show the selected SPOD mode shapes at $St_L$ = 0.12 for all three cases. In cases RN0 and RN4, the broadband low-frequency unsteadiness is primarily associated with the oscillation of shock wave system, encompassing the separation shock wave, bow shock wave, transmission shock wave and the Mach shock wave reflection occurring downstream of the first separation bubble, as well as the shear layer of the separation bubble. In case RN8, although a dominant frequency in the low-frequency region is observed, the SPOD mode shape indicates that the unsteady structures corresponding to this frequency are remarkably similar to those observed in cases RN0 and RN4, predominantly characterised by wave structures and the shear layer of the separation bubble.

Figure 9.

( $a$ ) The SPOD leading mode energy distributions for all cases. Selected SPOD mode shapes at $f$ = 2197 Hz for cases ( $b$ ) RN0, ( $c$ ) RN4 and ( $d$ ) RN8, derived from the schlieren sequences at the sampling rate of 75 kfps.

Figure 10.

( $a$ ) The SPOD leading mode energy distributions for all cases. Selected SPOD mode shapes for cases ( $b$ ) RN0, ( $c$ ) RN4 and ( $d$ ) RN8, derived from the schlieren sequences at the sampling rate of 300 kfps.

As for the unsteadiness observed on the second cone, figure 10 presents the SPOD results derived from the schlieren sequences captured at the sampling rate of 300 kfps. The energy distributions of SPOD leading mode depicted in figure 10( $a$ ) reveal the presence of broadband low-frequency characteristics in cases RN0 and RN4, alongside the notable peak frequency at $St_L$ = 0.113 ( $f$ = 2.05 kHz) in case RN8. This is attributed to the inclusion of the bow shock wave, transmission shock wave and portions of the separated shock wave within the field of view. In the higher-frequency region, the analysis indicates that the leading modal energy begins to rise at approximately $St_L$ = 2 for cases RN0 and RN4, and at $St_L$ = 1 for case RN8. In cases RN0 and RN4, the peak frequencies of the leading modal energy are very close, approximately falling within the range around $St_L$ = 3.9–4.4 ( $f$ = 70–80 kHz ). In comparison, in case RN8, the peak of the leading mode energy in the high-frequency range occurs at slightly lower frequency around $St_L$ = 3.3–3.9 ( $f$ = 60–70 kHz). Figure 10( $b$ , $c$ , $d$ ) illustrate the selected modal shapes corresponding to $St_L$ = 4.1 ( $f$ = 75 kHz) for cases RN0 and RN4, and $St_L$ = 3.6 ( $f$ = 65 kHz) for case RN8. The SPOD mode shapes indicate that, in these cases, high-frequency unsteadiness is primarily associated with the flow structures along the outer boundary of the supersonic jet. Based on their location within the high-shear region and their coherent vortical appearance in schlieren snapshots (figures 5 and 10), these high-frequency structures are suggestive of an instability intrinsic to the shear layer. Their characteristics are plausibly linked to a Kelvin–Helmholtz type instability (Knisely Reference Knisely2016; Tumuklu et al. Reference Tumuklu, Levin and Theofilis2018a , Reference Tumuklu, Theofilis and Levinb ).

Consequently, the SPOD analysis of flow structures captured through schlieren visualisation reveals that the unsteady phenomena in the three cases manifest in two primary forms. The first form consists of low-frequency unsteady oscillations associated with the shock wave system and the separation bubble. The second form involves large-scale vortices along the outer boundary of supersonic jet corresponding to the high-frequency unsteadiness, which plausibly arises from the Kelvin–Helmholtz instability of the strong shear layer on the second cone. Notably, unlike the broadband low-frequency oscillations observed in cases RN0 and RN4, case RN8 exhibits a distinct dominant frequency within the low-frequency range – a phenomenon not previously documented. The underlying mechanism for this dominant frequency will be further explored in § 4.

Figure 11.

Pressure signals captured by the Kulite sensors and PSD distributions of the pressure signals for cases ( $a{,}b$ ) RN0, ( $c{,}d$ ) RN4 and ( $e{,}f$ ) RN8.

Notably, the spectral analysis of the schlieren images primarily highlights the spatial unsteady structures in the meridian plane. In this study, dynamic wall pressure measurements were also conducted using pressure sensors and fast PSP technique. Actually, the wall surface parameters are intrinsically linked to the spatial flow field structures, prompting us to perform a dynamic analysis of the pressure fluctuations on the wall surface. As mentioned in § 2.1.4, the frequency response of pressure measurements from the NS2 sensors is influenced by the steel tube, leading us to focus solely on the pressure signals obtained from the four Kulite sensors and fast PSP technique. In the current study, the Kulite sensors and PSP system were sampled at 50 kHz and 10 kHz, respectively. Consequently, the effective upper frequency for spectral analysis was governed by the corresponding Nyquist frequencies (25 kHz and 5 kHz), restricting the dynamic analysis of wall pressure to low-frequency unsteady phenomena in this section. Figure 11( $a$ , $c$ , $e$ ) present the normalised pressure fluctuation signals (relative to local mean pressure) recorded by the four Kulite sensors on the second cone. In all three cases, the first sensor is positioned near the peak pressure location, while the remaining three sensors are situated beneath the supersonic jet, as illustrated in figure 5 (the sensor positions are indicated by white markers). It is evident that the normalised fluctuation amplitude at the first sensor is significantly higher than those at the other three locations. To quantify these observations, power spectral density (PSD) analysis of $p'/p_{\textit {local}}$ was conducted using the Welch method, as shown in figure 11( $b$ , $d$ , $f$ ). The results reveal that, in the low-frequency range, the energy content at the first sensor is approximately one to two orders of magnitude greater than those of the other three sensors. This finding indicates that the supersonic jet on the second cone exhibits very weak low-frequency unsteadiness, which aligns well with the SPOD results presented in figure 9. Additionally, the PSD distributions for the first sensor in case RN8 reveal a distinct peak at around $St_L$ = 0.118 ( $f$ = 2148 Hz), whereas cases RN0 and RN4 display broadband low-frequency characteristics.

Figure 12.

Root mean square (r.m.s.) contours of wall-pressure fluctuations obtained from the PSP and the corresponding distribution along the meridian lines for cases ( $a{,}d$ ) RN0, ( $b{,}e$ ) RN4 and ( $c{,}f$ ) RN8. The red dots denote the r.m.s. values measured by the first Kulite sensor.

It is well established that the spectral analysis of schlieren sequences captures unsteady oscillations in the meridian plane. However, previous studies on shock-induced separations have demonstrated that low-frequency unsteadiness can also manifest in spanwise movements within nominally two-dimensional flows (Sidharth et al. Reference Sidharth, Dwivedi, Candler and Nichols2018; Cao et al. Reference Cao, Hao, Klioutchnikov, Olivier and Wen2021b ; Hao et al. Reference Hao, Cao, Wen and Olivier2021, Reference Hao, Fan, Cao and Wen2022). Although the time-averaged pressure contours shown in figure 7 do not reveal obvious large-scale three-dimensional flow structures, this analysis aims to further specifically investigate the potential presence of low-frequency unsteady oscillations in the azimuthal direction. Figure 12 presents the r.m.s. of the wall pressure fluctuations ( $p'_{\textit {rms}}$ ). Figure 12 $(a-c)$ show contours of $p'_{\textit {rms}}$ normalised by the maximum of mean pressure ( $p_{max}$ ), obtained from the PSP data. Figure 12 $(d{-}f)$ plot the corresponding $p'_{\textit {rms}}$ distribution along the central meridian line. For comparison, the r.m.s. values measured by the first Kulite sensor are marked by red symbols. The comparison shows that the PSP derived $p'_{\textit {rms}}$ levels are approximately 40 $\,\%$ lower than the sensor measurements. Actually, this difference aligns with the frequency-dependent attenuation inherent to the fast PSP technique. As noted in § 2.1.4, the PSP attenuation is approximately 3.3 dB at 2 kHz (amplitude reduction to approximately 70 $\,\%$ ) and 5.6 dB at 4 kHz (amplitude reduction to approximately 50 $\,\%$ ). It is important to note that the frequency-dependent attenuation is expected to be spatially uniform, a condition intrinsically linked to the homogeneity of the PSP coating itself (Funderburk & Narayanaswamy Reference Funderburk and Narayanaswamy2019). This expectation is supported by the absence of any non-physical artefacts in the time-averaged PSP pressure fields, indicating a high degree of spatial uniformity in the applied paint layer. Given this consistent PSP application and optical set-up, the data are employed in this study primarily for the qualitative observation of spatial mode shapes, effectively complementing the pointwise measurements provided by conventional sensors. A quantitative comparison of unsteady characteristics across different cases, which would require accounting for the precise frequency-dependent attenuation in amplitude, is beyond the scope of the present work. In figure 12, it is of interest that distinct dual peaks are observed within the $p'_{\textit {rms}}$ contours on the second cone for cases RN0 and RN4. In contrast, for case RN8, the first peak diminishes significantly, becoming barely discernible. For case RN0, the two high $p'_{\textit {rms}}$ regions are located at approximately $x/L$ = 0.90–0.95 and $x/L$ = 0.98–1.02, while for case RN4 they occur at approximately $x/L$ = 0.93–0.97 and $x/L$ = 1.00–1.04. As illustrated in the Mach number contours and $C_{\kern-1pt f}$ distributions in figure 6, both cases RN0 and RN4 display two individual separation bubbles. According to the axial locations of the separation and reattachment points listed in table 2, it can be observed that the first high $p'_{\textit {rms}}$ region corresponds to the area between points $R1$ and $S2$ , while the second high $p'_{\textit{rms}}$ region corresponds the region downstream of point $R2$ (here, $S$ and $R$ denote the separation and reattachment points, respectively; the numbers 1 and 2 correspond to the first and second separation bubbles, respectively). The spatial correlation indicates that the dual-peak structure in the pressure fluctuations corresponds to the two spatially distinct separation regions. Critically, in the two cases, the second, much smaller separation bubble is directly induced by the impinging shock wave emanating from the TP TP3 (as visible in figures 5 $a$ and 5 $b$ ). Consequently, the oscillation of this bubble serves as a direct manifestation for the unsteadiness of the impinging shock, which is the primary cause of the distinct second $p'_{\textit {rms}}$ peak on the cone surface. In contrast, case RN8 features a large, coupled separation bubble. The thickened boundary layer at the shock impingement point substantially attenuates the shock’s direct impact on the cone surface. Here, the shock’s influence is indirectly transmitted by driving the low-frequency, large-scale oscillation of the entire bubble. This integrated dynamic results in a single, broad $p'_{\textit {rms}}$ peak, thereby explaining the absence of the dual-peak configuration in the $p'_{\textit {rms}}$ distribution observed in the other cases.

Spectral analysis for the surface pressure field was performed using the SPOD method. The datasets were segmented into 14 blocks with 75 $\,\%$ overlap, yielding a frequency resolution of 156 Hz. As illustrated in figure 13( $a$ ), the leading SPOD modes for cases RN0 and RN4 exhibit broadband low-frequency characteristics while case RN8 displays a pronounced spectral peak at $St_L$ = 0.112 ( $f$ = 2031 Hz). The selected SPOD mode shapes at $St_L$ = 0.112 are also shown in figure 13( $b$ , $c$ , $d$ ). Analysis of these modes reveals two key observations: no significant large-scale low-frequency three-dimensional unsteady structures are evident in any case, and the observed low-frequency unsteadiness primarily manifests as axial-direction oscillations, i.e. the axisymmetric oscillation mode. Moreover, the spatial distribution of modal energy demonstrates distinct patterns: in cases RN0 and RN4, the most energetic components are predominantly localised near the high $p'_{\textit{rms}}$ regions identified previously, with other regions exhibiting less distinct modal energy. The corresponding SPOD mode on the first cone exhibits a diffuse and less organised spatial structure, manifesting as a weak and indistinct alternating pattern. This stems from the relatively small separation bubble, which is characterised by inherently broadband, low-amplitude unsteadiness. In comparison, case RN8 is dominated by a strong, narrowband unsteadiness associated with a large separation bubble. This bubble effectively couples and amplifies the unsteady signal, transmitting vigorous dynamics from the downstream impinging shock upstream. The result is a highly coherent and spatially organised mode at a dominant frequency, with primary energy concentrated on the second cone and a visible pattern linked to bubble breathing on the first cone. This mode is clearly visualised as a strong red/blue alternation. Differences are also evident in the downstream region beneath the supersonic jet. In cases RN0 and RN4, the jet has a complex three-layer structure (as shown in figures 5 $a$ and 5 $b$ ). Reflections of upstream unsteady shock waves within this irregular environment are highly disturbed. This complexity, combined with the inherent broadband oscillations in this region, yields a diffuse SPOD mode with indistinct spatial organisation. Conversely, in case RN8, the jet possesses a simpler single-layer structure (figure 5 $c$ ), which supports clean and regular shock reflections that form a distinct wave pattern. Consequently, the wall-pressure mode shapes in figure 13( $d$ ), particularly the clear alternating red/blue bands in case RN8, serve as a direct spatial imprint of the regular reflection of upstream unsteady shock waves on the second cone surface.

Figure 13.

$(a)$ The SPOD leading-mode energy distributions. Selected SPOD mode shapes at $f$ = 2031 Hz for cases ( $b$ ) RN0, ( $c$ ) RN4 and ( $d$ ) RN8, derived from the pressure field data obtained using PSP technique.

4.1. Global stability analysis

In the GSA method, the flow field is assumed can be decomposed into two components: the steady base flow and small-amplitude unsteady perturbations. This decomposition is formally expressed as

(4.1) \begin{align} \begin{aligned} \boldsymbol{U}=\boldsymbol{\bar {U}+{U}^{\prime}} \end{aligned} , \end{align}

where $\boldsymbol{\bar {U}}$ and $\boldsymbol{{U}^{\prime}}$ denote the steady base flow and unsteady perturbations, respectively. Since the steady base flow satisfies the governing (2.1), substituting the decomposition (4.1) into (2.1) and neglecting higher-order terms yields the linearised perturbation equation,

(4.2) \begin{align} \begin{aligned} \frac {\partial \boldsymbol{{U}^{\prime}}}{\partial t}=\mathcal{J}\boldsymbol{{U}^{\prime}} \end{aligned} , \end{align}

where $\mathcal{J}$ is the Jacobian matrix, which corresponds to the linearisation of the RANS operator and can be evaluated based on the steady base flow. The perturbation $\boldsymbol{{U}^{\prime}}$ is assumed to take the form

(4.3) \begin{align} \begin{aligned} \boldsymbol{{U}^{\prime}}=\boldsymbol{\hat {U}}{{e}^{\lambda t}} \end{aligned} , \end{align}

here, $\boldsymbol{\hat {U}}$ and $\lambda$ represent the eigenfunction and eigenvalue, respectively. By substituting (4.3) into (4.2), we derive an eigenvalue problem that can be expressed in the following general form:

(4.4) \begin{align} \begin{aligned} \lambda \boldsymbol{\hat {U}}=\mathcal{J}\boldsymbol{\hat {U}} \end{aligned} . \end{align}

The eigenvalue problem is addressed using the implicitly restarted Arnoldi method, as implemented in ARPACK (Lehoucq, Sorensen & Yang Reference Lehoucq, Sorensen and Yang1998), and interfaced with the MUMPS library (Amestoy et al. Reference Amestoy, Buttari, Duff, Guermouche, L’Excellent and Uçar2011) for efficient linear algebra operations. The same discretisation method as in the flow solver is used except that the inviscid fluxes are calculated using a central scheme in smooth regions, while the modified Steger–Warming scheme is applied to compute inviscid fluxes at flow discontinuities detected by the modified Ducros sensor (Hendrickson, Kartha & Candler Reference Hendrickson, Kartha and Candler2018), thereby minimising numerical dissipation. Boundary conditions are consistent with those implemented in the numerical simulations. In the turbulent regions, the linearised source term for turbulence in the Spalart–Allmaras model is incorporated as described by Crouch, Garbaruk & Magidov (Reference Crouch, Garbaruk and Magidov2007), without resorting to simplifications such as the frozen eddy-viscosity approach (Carini et al. Reference Carini, Airiau, Debien, Léon and Pralits2017). Furthermore, since the computational grids from the numerical simulations are directly applied to GSA, an examination of grid convergence for GSA is also presented in Appendix B.

The present GSA solver has been effectively applied to both laminar and turbulent flow regimes (Hao et al. Reference Hao, Cao, Wen and Olivier2021, Reference Hao, Fan, Cao and Wen2022; Hao Reference Hao2023). For each case investigated in this study, we seek to compute 50 eigenvalues using the GSA solver. The eigenvalues consist of a real part $\lambda _{r}$ and an imaginary part $\lambda _{i}$ , which correspond to the growth rate and angular frequency of the perturbations, respectively. A flow is considered globally stable if $\lambda _{r}$ is less than zero for all computed eigenvalues. Furthermore, the eigenmodes can be classified based on the imaginary parts: those with $\lambda _{i}$ = 0 are termed stationary modes, indicating no oscillation in time, while those with $\lambda _{i}\neq 0$ are categorised as oscillatory modes.

The eigenvalue spectra for all three cases are presented in figure 14( $a$ , $c$ , $e$ ). To facilitate a direct comparison with experimental measurements, frequencies are expressed dimensionally ( $f$ ) and non-dimensionally ( $St_L$ ), rather than as angular frequency ( $\omega$ ). An analysis of the spectra reveals consistently negative real parts across all eigenvalues in the three cases, thereby confirming global stability of the flow systems. A closer examination of the spectra shows that the growth rate of the least stable mode progressively increases with increasing nose bluntness. Moreover, a notable change in the characteristics of the least stable mode is observed: for cases RN0 and RN4, the least stable mode remains stationary; however, for case RN8, the least stable mode exhibits an oscillatory nature. Interestingly, this unsteady dominant mode in case RN8 possesses a damping rate very close to zero, classifying it as a marginally stable mode. The corresponding frequency of this mode is approximately $St_L$ = 0.112 ( $f$ = 2030 Hz), which aligns remarkably well with the dominant frequency identified in the experimental results.

Figure 14.

Eigenvalues spectra and real parts of the density perturbation of the least stable mode for cases ( $a{,}b$ ) RN0, ( $c{,}d$ ) RN4 and ( $e{,}f$ ) RN8.

The eigenfunctions corresponding to the least stable modes in each case are illustrated in figure 14( $b$ , $d$ , $f$ ). For cases RN0 and RN4, where stationary modes dominate, the perturbation $\rho '$ is predominantly concentrated within the shock-dominated regions – specifically the separation shock wave, detached bow shock wave and Mach shock wave reflection. A significant presence of perturbation is concurrently observed in the shear layer enveloping the separation bubble, whereas the recirculation zone within the bubble exhibits minimal modal activity. Interestingly, although the least stable mode in case RN8 is oscillatory, the spatial distribution of the perturbation $\rho '$ shows some similarity to the stationary modes in cases RN0 and RN4, with peak amplitudes likewise localised within the shock wave structures. This consistent spatial localisation of least stable modes aligns closely with the shock-associated coherent structures identified via SPOD in prior experimental investigations, confirming the physical relevance of the computed eigenmodes.

4.2. Resolvent analysis

The GSA results demonstrate that all three base flows behave as a noise amplifier. To investigate the response of these base flows to the perturbations, a small-amplitude forcing term is added to (4.2),

(4.5) \begin{align} \begin{aligned} \frac {\partial \boldsymbol{{U}^{\prime}}}{\partial t}=\mathcal{J}\boldsymbol{{U}^{\prime}}+\mathcal{B} \kern-1pt \boldsymbol{{f}'} \end{aligned} , \end{align}

where the matrix $\mathcal{B}$ constrains the forcing to a localised region. We assume that the forcing and response are temporally harmonic,

(4.6) \begin{align} \begin{aligned} \boldsymbol{{f}'}=\boldsymbol{\hat {f}}{{e}^{i\omega t}},\quad \boldsymbol{{U}^{\prime}}=\boldsymbol{\hat {U}}{{e}^{i\omega t}} \end{aligned} , \end{align}

where $\omega$ is the angular frequency. Substituting (4.6) into (4.5) yields

(4.7) \begin{align} \begin{aligned} \boldsymbol{\hat {U}}=\mathcal{R}\mathcal{B}\boldsymbol{\hat {f}}, \mathcal{R}=\left ( i\omega \mathcal{I}-\mathcal{J} \right ) \end{aligned} , \end{align}

where, $\mathcal{R}$ and $\mathcal{I}$ are the resolvent and identity operators, respectively. The resolvent analysis aims to identify the forcing and response pair that exhibits the maximum energy amplification, or the optimal gain $\sigma$ , which is defined as

(4.8) \begin{align} \begin{aligned} {{\sigma }^{2}} ( \omega)=\underset {{\boldsymbol{\hat {f}}}}{\mathop {\max }}\,\frac {{{ \| {\boldsymbol{\hat {U}}} \|}_{E}}}{{{ \| \boldsymbol{B \kern-1.5pt \hat {f}} \|}_{E}}} \end{aligned} . \end{align}

The Euclidean norm and Chu energy norm are employed to calculate the energies of both the forcing and the response. This optimisation problem can be reformulated as an eigenvalue problem, as demonstrated by Bugeat et al. (Reference Bugeat, Chassaing, Robinet and Sagaut2019). To solve this eigenvalue problem, we utilise the power iteration method, which is well-suited for identifying the dominant eigenvalue associated with maximum energy amplification. The discretisation methods used in this process remain consistent with those implemented in GSA, ensuring coherence in the numerical framework.

Figure 15.

Optimal gain distributions and real parts of the density perturbations for the selected mode at specific frequency for cases ( $a{,}b$ ) RN0, ( $c{,}d$ ) RN4 and ( $e{,}f$ ) RN8.

The resolvent analysis is conducted based on the base flow over a wide frequency range to investigate its response to the environmental noise. Notably, guided by the scale decoupling assumption, the RANS simulations effectively capture low-frequency, large-scale structures – such as shock-wave motion, while high-frequency, small-scale vortical structures are modelled using turbulence models. Consequently, this preliminary resolvent analysis focuses exclusively on the low-frequency range in the current study. In the resolvent analysis, perturbations were applied throughout the entire computational domain using a methodology similar to that employed by Sartor et al. (Reference Sartor, Mettot, Bur and Sipp2015a ,Reference Sartor, Mettot and Sipp b ). Figure 15( $a$ , $c$ , $e$ ) illustrates the variation of the optimal gains with frequency. The results reveal that for cases RN0 and RN4, the maximum optimal gain occurs at the ultralow frequency, with a gradual decrease as frequency increases. This behaviour aligns with trends observed in the experimental results and has been reported in previous studies of laminar separation (Bugeat et al. Reference Bugeat, Robinet, Chassaing and Sagaut2022). Strikingly, case RN8 exhibits fundamentally different dynamics: after an initial decline, the optimal gain surges to a sharp peak at $St_L$ = 0.112. To elucidate the low-frequency flow characteristics observed in RN0 and RN4, figures 15 $(b)$ and 15 $(d)$ present the optimal response modes at $St_L$ = 0.001, corresponding to the red points marked in figures 15( $a$ ) and 15( $c$ ). The response modes in the two cases resemble the corresponding least stable modes predicted by GSA in figures 14( $b$ ) and 14( $d$ ). To quantitatively compare the modes in resolvent analysis and GSA, the projection of them is calculated, which is defined as

(4.9) \begin{align} \begin{aligned} \xi =\frac {\big | {{\big \| {\boldsymbol{\hat {U}}_{GSA}},{\boldsymbol{\hat {U}}_{IO}} \big \|}_{E}} \big |}{{{\big \| {\boldsymbol{\hat {U}}_{GSA}} \big \|}_{E}}\times {{\big \| {\boldsymbol{\hat {U}}_{IO}} \big \|}_{E}}} \end{aligned} , \end{align}

where $\boldsymbol{\hat {U}}_{GSA}$ and $\boldsymbol{\hat {U}}_{IO}$ denote the least stable GSA mode and optimal resolvent mode, respectively. The notation $\|\boldsymbol{\cdot }\|_{E}$ represents the Chu norm (Chu Reference Chu1965). Thus, the projection coefficient $\xi$ = 0 signifies the two modes are orthogonal, while $\xi$ = 1 indicates a perfect alignment between the two modes. The calculation shows that $\xi$ = 0.875 and 0.915 for cases RN0 and RN4. This result reveals a strong alignment between the least stable GSA mode and the optimal resolvent mode for cases RN0 and RN4. It indicates that their broadband, low frequency response is a modal amplifier phenomenon. The globally stable flow preferentially amplifies the low frequency content of the environmental noise, with maximum gain occurring at the lowest frequencies where the system’s dynamics are governed by its least damped stationary mode. This is consistent with established results for laminar separation (Bugeat et al. Reference Bugeat, Robinet, Chassaing and Sagaut2022). For the peak frequency identified in case RN8, the optimal response pattern in figure 15 $(f)$ exhibits nearly perfect agreement with the least stable global mode predicted by GSA in figure 14( $f$ ), with the projection coefficient $\xi$ = 0.991. This confirms that the sharp spectral peak in RN8 originates from a true modal resonance. Here, the flow possesses a marginally stable oscillatory global mode that selectively and powerfully amplifies disturbances at its natural frequency (∼2 kHz). This contrasts with the broadband response in RN0/RN4, which arises from the modal amplification of noise by a stationary global mode.