5.1. Two-dimensional flow conditions

Table 1 summarises the flow conditions for a choked 2-D case with $J=1.17$ (slightly above the critical threshold). These conditions will later be compared with 3-D simulations to quantify the dimensional effects on choking. Despite their limitations, the 2-D results enhance the intuition about the underlying physics, such as SBLIs and blockage effects, which remain relevant in 3-D flows.

Table 1.

The flow conditions of the two-dimensional jet in supersonic crossflow simulation.

Table 1. Long description

The table has two rows of data and three columns: temperature in kelvin, static pressure in bar, and Mach number. The crossflow row gives a temperature of 1400 kelvin, a static pressure of 0.75 bar and a Mach number of 3.2. The jet row gives a temperature of 300 kelvin, a static pressure of 9 bar and a Mach number of 1.

The 2-D channel configuration, adapted from the 3-D computational framework (§ 3), exhibited choking at significantly lower momentum flux ratios (J) than its 3-D counterpart. Key modifications included: (i) extending the slip-wall condition to x = 7 cm to suppress boundary-layer development, (ii) positioning the injector at x = 8 cm for flow stabilisation and (iii) applying AMR prior to x = 72 cm with an upper slip wall to isolate inviscid shock dynamics. These adjustments eliminated 3-D effects – turbulent mixing, sidewall vortices and boundary-layer separation – thereby reducing momentum losses and advancing the choking onset. The 1152 cm domain length ensured outflow boundary independence, while the slip condition simplified the flow topology, yielding coherent Mach stems and unambiguous choking signatures.

The observed leftward shift in the choking threshold underscores the dominant role of dimensionality in compressible flow resistance. In two dimensions, the absence of spanwise dissipation pathways allows shock trains to stabilise with weaker adverse pressure gradients, explaining the reduced J requirement. The AMR-resolved flow fields confirmed quasi-steady shock structures, validating the numerical framework’s ability to capture inviscid choking criteria without 3-D artefacts. While this simplification clarifies fundamental mechanisms – such as the direct correlation between choking and bulk flow deceleration – it also highlights the necessity of reintroducing 3-D effects for practical systems, where turbulent mixing and boundary-layer interactions critically modulate the choking limit.

The flow evolution was systematically analysed through temporal snapshots at five characteristic instants, with pressure contours depicted in figure 4. The simulations reveal that Mach stem formation commenced immediately upon initialisation (figure 4 a), indicating rapid establishment of choking conditions. Subsequent temporal evolution demonstrated progressive downstream pressure accumulation (figure 4 a–c), which ultimately generated upstream-propagating disturbances culminating in unstart initiation. This process manifested through two distinct yet coupled phenomena: (i) the formation of a normal shockwave upstream of the injector and (ii) the temporal development and upstream migration of the downstream-originating shock structure.

Figure 4.

Temporal evolution of pressure contours (in barye) illustrating the choking process in the 2-D JISCF configuration. Snapshots at (a) t = 0.3 ms, (b) t = 0.7 ms, (c) t = 1.59 ms, (d) t = 2.389 ms and (e) t = 3.3 ms show the formation of the Mach stem, downstream pressure accumulation and the upstream propagation of disturbances leading to unstart. The locations of pressure probes used for spectral analysis are indicated.

The observed dynamics highlights a critical feedback mechanism wherein downstream pressure buildup drives upstream shock formation, which in turn amplifies Mach stem growth. This mutual reinforcement between downstream and upstream wave systems ultimately leads to complete flowpath unstart. The time-resolved visualisation provides clear evidence of the shock interactions that characterise this transitional process, with particular emphasis on the spatial progression of wavefronts and their role in establishing the final unstart configuration.

In short, the 2-D configuration serves as a controlled numerical experiment to isolate the inviscid core of the global mode. Slip conditions on the upper wall are applied to: (i) eliminate boundary-layer disturbances that could mask system-scale dynamics, (ii) suppress upstream acoustic propagation through the upper wall while maintaining an upstream communication pathway through the no-slip lower wall (a separate supersonic confinement effect) and (iii) verify that the mode originates from bulk shock–flow interaction, not from random wall-generated turbulence. This approach clarifies the fundamental wave feedback responsible for choking.

Table 2.

The flow conditions of the three-dimensional jet in supersonic crossflow simulations at varying momentum flux ratios (J).

Table 2. Long description

The table has two rows of data and three columns: temperature in kelvin, static pressure in bar, and Mach number. The crossflow row gives a temperature of 1400 kelvin, a static pressure of 0.75 bar and a Mach number of 2.7. The jet row gives a temperature of 300 kelvin and a Mach number of 1, with three static pressures of 7.6659, 15.33 and 30 bar. These correspond to momentum flux ratios of 1.5, 3.0 and 5.87 respectively.

5.2. Three-dimensional flow conditions

The 3-D refined computational domain extends from x = 0 to 40 cm. Table 2 summarizes the flow conditions for the 3-D cases. A brief slip-wall region ( $x\lt 2$ cm) is applied at the inlet to prevent spurious numerical reflections (Rising et al. Reference Rising, Ching and Johnson2026), after which no-slip boundary conditions are enforced on all walls to allow natural boundary-layer development. The injector, with a radius of 1.6 mm, is positioned at x = 8.4 cm, providing a development length of $\approx 6.3$ cm from the start of the no-slip region. Over this distance, the local theoretical Reynolds number based on streamwise distance from the effective origin reaches $Re_x\approx 4.7\times 10^5$ – a value that lies within the natural transitional range for a flat-plate boundary layer. No artificial tripping is employed; the inflow profile develops self-consistently from the governing equations under the specified boundary conditions. This approach avoids introducing arbitrary trip parameters and ensures the inflow is determined solely by the governing equations and specified boundary conditions. The flow was initialised to the inflow conditions.

To ensure accurate flame and shear-layer resolution, this work adopts the experimental observations of Gamba & Mungal (Reference Gamba and Mungal2015), who report an average OH layer thickness of 0.3 mm – a scale which is associated with the shear-layer thickness in the simulations. This follows the resolution criteria established by Attili et al. (Reference Attili, Bisetti, Mueller and Pitsch2014), which prescribe 3–10 grid points across the shear layer. This resolution guarantees proper capture of both the flame structure and mixing dynamics, while remaining computationally tractable for the domain size under consideration.

To assess the choking behaviour across different momentum flux ratios (J), pressure iso-surfaces were analysed at 1.5 ms to identify the onset of flow obstruction. Figure 5 illustrates the pressure distribution using three distinct iso-surfaces (1, 3 and 5 bar). For J = 1.5, the flow remains largely unperturbed, with pressure levels consistently near 1 bar throughout the domain (figure 5 a). In contrast, at J = 3, a moderate pressure rise is observed downstream, particularly along the lower wall, though the peak pressure remains below 3 bar. This suggests increased flow resistance but not yet sufficient to induce full choking.

Figure 5.

Instantaneous pressure iso-surfaces (1, 3 and 5 bar) in the 3-D JISCF at t = 1.5 ms for momentum flux ratios (a) J = 1.5, (b) J = 3.0 and (c) J = 5.87. The onset of significant flow blockage and high-pressure regions downstream signals the transition to choking at J = 5.87.

The most pronounced effects are seen at J = 5.87, where a clear pressure buildup emerges downstream, with the 3 bar iso-surface extending between x = 30 and 40 cm. This localised high-pressure region signals the initiation of choking, as the flow encounters significant blockage due to the jet’s momentum dominance. While the pressure does not reach 5 bar within the observed time frame, the progressive rise with increasing J underscores the role of momentum flux ratio in dictating choking onset. These findings align with theoretical expectations, where higher J values enhance jet penetration and flow restriction, ultimately leading to choking at sufficiently elevated ratios.

The heat release characteristics were quantified through spatial averaging and integration of the volumetric HRR field. Spanwise and wall-normal averaging yielded the local HRR distribution $\langle \dot {q} \rangle _{y,z}$ (units: $\mathrm{erg\,cm^{-3}\,s^{-1}}$ ), while streamwise integration of these profiles provided the total power output $\dot {Q}_{total}$ (units: $\mathrm{erg\,s^{-1}}$ ). This dual approach captures both the local combustion intensity and global energy release, with the integration performed across the entire reaction zone to account for all chemical heat release.

7.1. Spectral signature of choking

The slip boundary condition at the upper wall eliminated boundary-layer vortices, enabling unambiguous identification of pressure fluctuations behind the Mach stem and secondary reflected shock in the 2-D case (figure 9). Power spectra revealed distinct frequency signatures tied to these features: a dominant 305 Hz peak behind the Mach stem, and weaker modes downstream. The low-frequency mode will be studied in the future discussions, while the secondary frequency emerges from some large-scale convective structures or a nonlinear coupling, but it is out of the scope of this study.

Figure 9.

Power spectral density (PSD) of pressure perturbations at four streamwise locations (x = [9, 20, 30, 50] cm, and y = 3.9 cm) in the 2-D simulation. The dominant 305 Hz peak is a global feature of the choked flow.

Time-resolved pressure signals were extracted from selected locations along the wall. Due to non-uniform time-step outputs inherent to AMR and Courant–Friedrichs–Lewy-based time stepping, the non-uniform fast Fourier transform (NUFFT) algorithm was employed to compute the spectral content. Each signal consisted of at least 1000 temporally resolved data points spanning a physical window of approximately 3 to 5 ms for the unchoked and the supercritical case.

The NUFFT implementation allows for accurate estimation of spectral components without resampling or interpolation, avoiding spectral distortion or aliasing artefacts associated with uneven sampling. This approach enabled robust identification of dominant low-frequency instabilities associated with shock-induced oscillations and global flow unsteadiness.

The 305 Hz frequency is that of an emergent mode that the entire system locks into. The Mach stem oscillates and contributes energy, but it is a participant in this system-wide resonance, not the sole source of its specific frequency. The elevated energy in the region near the Mach stem confirms its role as a key driver of the instability. The mode itself arises from a balance between large-scale wave propagation and confinement effects that coordinate the entire flow field.

The secondary mode emerges downstream, produced by convective structures as the flow interacts with the walls. This frequency’s appearance coincides with a significant reduction ( ${\gt} 35\,\%$ ) in the primary mode’s amplitude, illustrating an energy cascade from large-scale shock oscillations to finer-scale turbulence. The damping of the 305 Hz signal reflects convective losses, damping of nonlinear waves and vorticity production by Kelvin–Helmholtz instabilities, which disrupt coherent oscillations while enhancing mixing.

The persistence of the primary frequency downstream confirms its connection to global choking dynamics. This dual-scale behaviour aligns with compressible flow theory, where shock-driven instabilities (305 Hz) seed shear-layer vortices, redistributing energy across scales. The spectra collectively quantify how choking-induced oscillations evolve into dissipative structures. The persistence of the dominant mode’s frequency at downstream locations suggests it is a global feature of the flow, rather than a local, spatially dependent instability. In turn, the Lagrangian point can be chosen to be located anywhere downstream of the Mach stem. However, this claim requires further investigations using statistical analysis.

In figure 10(d), the phase angles between four probe locations were analysed to elucidate wave propagation at the dominant mode. The analysis revealed a complex, bidirectional wave structure (table 3). At Probe 3, a phase lead of $64^\circ$ relative to Probe 1 was observed, corresponding to a downstream phase speed of 358.2 m s⁻¹ (figure 10 e). This subsonic velocity, characteristic of large-scale convective structures, indicates a mechanism for downstream energy transport from the primary shock system.

Figure 10.

Statistical analysis of the 2-D choked flow. (a) Correlation matrix of pressure signals from four probes. (b) Magnitude-squared coherence spectrum. (c) Streamwise velocity phase difference at dominant low-frequency mode. (d) Pressure phase difference at dominant low-frequency mode. (e) Phase speed between probe pairs, revealing bidirectional (downstream convective and upstream acoustic) wave propagation.

Table 3.

Phase angle and phase speed between the four sensors.

Table 3. Long description

The table has six rows, one per probe pair, and four columns: phase speed in metres per second, pressure phase angle, streamwise velocity phase angle, and the difference between the pressure and velocity phase angles. All angles are in degrees. Pair 1 to 2: 148.3, 81.5, 61.7, 19.8. Pair 1 to 3: 358.2, 64.4, 54.0, 10.4. Pair 1 to 4: 1116.6, 40.3, 36.1, 4.2. Pair 2 to 3: negative 642.5, negative 17.1, negative 7.7, negative 9.4. Pair 2 to 4: negative 800.8, negative 41.1, negative 25.6, negative 15.5. Pair 3 to 4: negative 913.2, negative 24.0, negative 17.9, negative 6.2. Positive phase speeds correspond to downstream convective propagation and negative values to upstream acoustic propagation.

At Probe 2, a phase lead of $81.5^\circ$ – midstream probes – and nearly half a cycle – was observed relative to Probe 1. The corresponding phase speed of 148.3 m s⁻¹ (figure 10 e) indicates convective communication mediated by large-scale vortical structures. This downstream propagation mechanism is consistent with that identified between Probes 1 and 3.

Conversely, the measured phase lag of $-24^\circ$ between Probes 4 and 3 yields a phase speed of −913 m s⁻¹ (figure 10 e), a value that matches the local speed of sound and unambiguously identifies an upstream-travelling acoustic wave. This wave is interpreted as the key component of a feedback loop: it is the acoustic precursor generated by downstream shock oscillations, which travels upstream to trigger new instabilities, completing a cycle that sustains the emergent global mode.

The phase relationships for other probe pairs further reinforced the presence of upstream acoustic feedback. Relative to Probe 2, Probe 3 lagged by $-17^\circ$ and Probe 4 by $-41.1^\circ$ , yielding upstream phase speeds of −642 and −800.8 m s⁻¹ (figure 10 e), respectively. These negative velocities are a clear signature of slow acoustic communication propagating upstream from the choking zone.

The net communication across the entire domain is captured by the relationship between Probe 4 and Probe 1, where a $40.3^\circ$ phase lead and a high phase speed of 1116.4 m s⁻¹ were measured. This value, which aligns with the local speed of sound in the high-temperature sections of the flow, confirms the presence of system-spanning acoustic communication with an average phase speed of −122.2 m s⁻¹ (figure 10 e). This global acoustic wave is responsible for phase-locking the entire flow field, from the upstream Mach stem to the downstream choking zone, on the same dominant low-frequency mode.

To definitively characterise the nature of the observed waves, the phase relationship between the streamwise velocity and pressure fluctuations ( $\phi _p - \phi _u$ ) was examined. This analysis is critical to distinguish between active energy transport and passive kinematic phase speeds; a phase difference near $0^\circ$ confirms progressive acoustic waves capable of transporting energy, whereas a value near $\pm 90^\circ$ indicates standing waves with no energy transport. The measured $\phi _p - \phi _u$ phase differences across the four sensor locations provide definitive evidence that moves beyond a purely kinematic interpretation, revealing a dynamic and spatially evolving acoustic environment where each value corresponds to a distinct wave mechanism.

In table 3, the phase speed of the pressure signal, the phase angle of the pressure and the streamwise velocity, and the difference between the phase angle of the velocity and the pressure are given. The significant $+19.8^\circ$ and $10.4^\circ$ ( ${\lt} 30^\circ$ ) phase differences, measured in the region of the 148.3 and 358.2 m s⁻¹ downstream phase speed, are the hallmark of a progressive, travelling convective–acoustic wave for Probes 2 and 3 relative to Probe 1. This value, close to the theoretical ideal of $0^\circ$ , confirms that the subsonic wave propagation from the primary shock system is an efficient, energy-transporting mechanism, where pressure and velocity fluctuations are nearly in phase as they are carried by the mean flow.

Conversely, the distinct $-6.2^\circ$ , $-15.5^\circ$ and $-9.4^\circ$ phase differences are the definitive signature of the upstream-propagating acoustic feedback wave, corresponding to the −913 , −800.8 and −642.5 m s⁻¹ phase speeds (figure 10 e). This slight negative value arises from the wave fighting against the mean flow, where fluid inertia causes a slight lag in the velocity response compared with the pressure forcing. This measurement provides conclusive proof that the high-speed upstream propagation is an active, energy-carrying wave communicating the downstream blockage back towards the shock system.

The intermediate $+4.2^\circ$ phase difference corresponds to the net 1116 m s⁻¹ acoustic communication measured across the entire domain (table 3). This value represents the integrated effect of the wave field – a superposition of the downstream convective structures and the upstream acoustic feedback. Its positive, near-zero value confirms that the system-scale phase locking is achieved through genuine progressive acoustic waves, capable of spanning the combustor at speeds matching the local sound speed in high-temperature regions.

In summary, the $p'{-}u'$ phase analysis of the 2-D simulation elevates the observed spectral dynamics from a set of frequencies to a coherent physical mechanism. The measured angles in table 3 collectively prove the existence of a bidirectional acoustic energy circuit in the simplified geometry. This circuit is the essence of a dual-wave mechanism that coordinates the global instability: convective transport distributes energy downstream, while intense acoustic feedback propagates upstream, creating a self-sustaining loop. This framework provides the communication channel that phase-locks the shock system and confinement boundaries into a globally resonant mode, underscoring a mechanism absent from classical 1-D theory. The diagnostic power of this $p'{-}u'$ phase analysis establishes the methodology that will be applied to the 3-D flow to characterise the wave dynamics of its dominant low-frequency mode.

The statistical cross-correlation of the pressure signals was calculated. Figure 10(a) presents the correlation matrix, revealing that midstream Probes 2 and 3 and downstream Probe 4 are highly correlated with each other ( $Corr. \gt 0.7$ ). In contrast, Probe 1 shows only moderate correlation with the midstream probes (0.405 and 0.512) but stronger correlation with the downstream probe ( $Corr. = 0.745$ ). For the flow evolution from the highly nonlinear chaotic region near Probe 1 – characterised by Mach stem oscillations, bow shock dynamics and vortex shedding – the incoherent nonlinear waves decay spatially. Local structures would show monotonic correlation decay with distance. This selective dissipation preserves the coherent global mode, explaining both the low correlation between Probe 1 and midstream probes and the recovery of higher correlation further downstream.

The magnitude-squared coherence spectrum showed dominant peaks in the low-frequency range (0–1 kHz), with coherence decaying at higher frequencies – a characteristic signature of a global mode where system-scale dynamics dominates over local turbulent fluctuations. This inverse relationship between frequency and coherence further confirms the global nature of the 305 Hz oscillation.

To further investigate the nonlinear dynamics, magnitude-squared bicoherence analysis was performed on the pressure signals from all four probes. Figure 11(a) shows the magnitude-squared bicoherence contour plot at Probe 1, revealing three significant peaks indicative of nonlinear triad interactions. Two of these peaks involve synchronised coupling between low-frequency and moderate-frequency modes (approximately 1700 and 305 Hz), suggesting energy transfer across scales. The magnitude-squared bicoherence was calculated using the following equation:

(7.1) \begin{equation} {\textit{bic}}^2(f_1,f_2)=\frac {|\sum _n X_n(f_1)X_n(f_2)X_n^*(f_1+f_2)|^2}{\sum _n|X_n(f_1)X_n(f_2)|^2\sum _n|X_n^*(f_1+f_2)|^2}, \end{equation}

where $X_n(f)$ is the Fourier transform at this specific frequency of segment n, and the asterisk indicates the conjugate of a function (Maccarone & Schnittman Reference Maccarone and Schnittman2005). The overlapping ratio which was used to generate the following results was 0.5. The magnitude-squared coherence was calculated using the following equation:

(7.2) \begin{equation} \text{MSC} = \frac {|P_{xy}|^2}{|P_{xx}||P_{yy}|}, \end{equation}

where $P_{xx}$ is the auto-power spectral density (PSD) of signal x(t), $P_{yy}$ is the auto-PSD of signal y(t) and $P_{xy}$ is the cross-PSD of the signals x(t) and y(t) (González et al. Reference González, Rodríguez, Sagartzazu, Schumacher and Isasa2010).

Figure 11.

Magnitude-squared bicoherence analysis of pressure signals at Probe 1 for the 2-D simulation. (a) Contour plot showing nonlinear triad interactions. (b) Cross-section at f/2. (c) Statistical distribution of bicoherence values, indicating the prevalence of nonlinear coupling upstream (Probe 1).

To further clarify the role of magnitude-squared bicoherence and magnitude-squared coherence in this study, this work defines each metric in physical terms. The magnitude-squared bicoherence, ${\textit{bic}}^2(f_1,f_2)$ , quantifies the fraction of energy at the sum frequency $f_1+f_2$ that arises from quadratic phase coupling between $f_1$ and $f_2$ . A value of ${\textit{bic}}^2=0.8$ implies that the bicoherence $bic=\sqrt {0.8}\approx 0.89$ , indicating approximately $90\,\%$ probability that the mode at $f_1+f_2$ results from a nonlinear interaction between $f_1$ and $f_2$ . In this context, a dominant frequency X may be generated via self-interaction – where two waves of frequency $ {X}/{2}$ and similar phase produce X, or through triad interaction, where two distinct frequencies $f_1$ and $f_2$ combine such that $f_1+f_2=X$ .

In the present study, however, the observed global mode is unlikely to originate from such harmonic interactions. This is supported by its bidirectional convective–acoustic nature, spatial coherence and linear independence. Moreover, the mode emerged as a dominant spectral peak only under choked conditions in both 2-D and 3-D simulations, indicating a direct connection to the choking process.

As for the magnitude-squared coherence, a value of 0.81 indicates that $81\,\%$ of the power at a given frequency in one signal can be predicted by the other signal, reflecting a high degree of linear correlation and shared dynamics. The remaining $19\,\%$ is attributed to noise or unrelated fluctuations.

Figure 11(b) demonstrates that magnitude-squared bicoherence values range between 0.6 and 0.635 for low-frequency interactions (0–305 Hz), confirming substantial nonlinear coupling in this frequency band. This supports two possible interpretations: either the 305 Hz global mode emerges from nonlinear interactions between other modes and hence not being a global mode, or alternatively, nonlinear coupling actively feeds and sustains the global mode.

The latter interpretation is strengthened by the presence of significant magnitude-squared bicoherence peaks within the $305 \times 305$ Hz region (figure 11 a), where bicoherence values exceed 0.6. Statistical analysis reveals that while only $4.4\,\%$ of frequency pairs exhibit strong coupling ( $bic. \gt 0.6$ ), a substantial $41.9\,\%$ show moderate nonlinear interaction ( $bic. \gt 0.5$ ). This distribution indicates that the global mode is sustained by selective nonlinear energy transfer mechanisms.

The magnitude-squared bicoherence analysis of midstream probes (figures 12 a and 13 a) reveals remarkably low values ( $bic. \lt 0.15$ ), providing clear evidence of the global mode’s predominantly linear propagation characteristics in these regions. The absence of significant nonlinear mode interactions indicates that once established, the low-frequency oscillation propagates downstream with minimal energy transfer to other frequencies.

Figure 12.

Magnitude-squared bicoherence analysis of pressure signals at Probe 2 for the 2-D simulation. (a) Contour plot showing nonlinear triad interactions. (b) Cross-section at f/2. (c) Statistical distribution of bicoherence values, indicating the decay of nonlinear coupling in midstream regions (Probe 2).

Figure 13.

Magnitude-squared bicoherence analysis of pressure signals at Probe 3 for the 2-D simulation. (a) Contour plot showing nonlinear triad interactions. (b) Cross-section at f/2. (c) Statistical distribution of bicoherence values, indicating the decay of nonlinear coupling in midstream regions (Probe 3).

This linear behaviour is further emphasised by cross-sectional analysis at half the global mode frequency (figures 12 b and 13 b), which shows almost no nonlinear coupling with other oscillating modes. Statistical quantification confirms this observation: only $1\,\%$ of mode interactions at Probe 2 exhibit moderate coupling ( ${\textit{bic}}^2 \gt 0.5$ ), while Probe 3 shows complete absence of significant nonlinear interactions ( $0\,\%$ above this threshold).

To complete the nonlinear analysis, magnitude-squared bicoherence was computed for the downstream probe (Probe 4). The contour plot (figure 14 a) reveals re-excited nonlinear interactions, with distinct peaks between low-frequency (180 Hz) and moderate-frequency (1400 Hz) modes. This demonstrates that the global mode remains dynamically active far downstream, capable of reinitiating local nonlinearities in the flow recovery region, as a result for the presence of more than one peak on the 305 Hz line (figure 14 a).

Figure 14.

Magnitude-squared bicoherence analysis of pressure signals at Probe 4 for the 2-D simulation. (a) Contour plot showing nonlinear triad interactions. (b) Cross-section at f/2. (c) Statistical distribution of bicoherence values, indicating the generation of nonlinear waves due to the nonlinear interaction between the dominant mode and other low- and moderate-frequency modes in downstream regions (Probe 4).

Critically, the magnitude-squared bicoherence analysis shows no significant nonlinear interactions directly driving the dominant mode at this location (figure 14 b), indicating that the global mode now self-sustains without external nonlinear forcing – preserving its coherence through linear wave dynamics.

Statistical quantification confirms this interpretation: while $65.8\,\%$ of frequency pairs exhibit moderate nonlinear coupling ( ${\textit{bic}}^2 \gt 0.5$ ) and $34.7\,\%$ show strong coupling ( ${\textit{bic}}^2 \gt 0.6$ ), these interactions occur around rather than with the global mode. This suggests that the mode’s oscillation has transitioned from being nonlinearly driven to becoming a stable, self-sustaining wave that can passively modulate other flow instabilities.

To validate the Fourier-based identification of the global mode against the non-stationary transient, a time–frequency analysis was applied to the pressure signal at Probe 2. This midstream location is representative, as established by the linear propagation characteristics of the mode shown in the bicoherence analysis. The resultant spectrogram (figure 15 a) reveals a persistent concentration of spectral energy at a frequency of approximately 305 Hz throughout the duration of the choking transient. The temporal evolution of this mode’s normalised energy at each probe (figure 15 b) shows a sustained activity. This independent analysis confirms that the 305 Hz spectral peak identified via fast Fourier transform represents a coherent, long-lived global mode intrinsic to the choked flow state, rather than a transient artefact.

Figure 15.

Time–frequency analysis of the global mode of the 2-D case. (a) Time–frequency at midstream Probe 2 showing persistent energy concentration at 305 Hz. (b) Temporal evolution of normalised mode energy at all four probes, confirming sustained activity with transient suppression linked to shockwave passage.

The transient suppression of the global mode, sequenced from upstream to downstream and back (figure 15 b), supports the aforementioned acoustic–convective feedback loop claim in the time–frequency domain. The observed propagation (P1 $\rightarrow$ P4 $\rightarrow$ P3 $\rightarrow$ P2) is consistent with a disturbance originating in the upstream shock system, propagating downstream through the convective path and subsequently reflecting or generating an upstream-travelling wave through the acoustic path. This event temporarily disrupts the phase locking of the global mode at each local measurement point as the disturbance passes. The mode’s robust recovery following this major perturbation underscores its stability as the system’s dominant attractor state. This analysis reveals that the 305 Hz global mode is not a passive spectral feature but an active, dynamic entity that interacts with – and can be momentarily disrupted by – finite-amplitude waves within the very feedback circuit that sustains it. The magnitude-squared bicoherence reveals quadratic interactions not captured in linear cross-spectra (Hunt & Gamba Reference Hunt and Gamba2019).

To complete the discussion, the extension to the 3-D problem is performed. The three distinct cases were selected to show that there is no low-frequency dominant mode appearing in the unchoked cases. The PSD analysis of the 3-D flow reveals a fundamental reorganisation of energy as the momentum flux ratio ( $J$ ) increases, culminating in an order-of-magnitude amplification of fluctuation energy at $J = 5.87$ . This transition is marked by the emergence of a dominant 600 Hz peak ( $\textit{St}\approx 0.014$ ) near the upper wall, as shown in figure 16. This low Strouhal number signifies a global instability mode, which is spectrally distinct from the broadband turbulent signature observed at lower $J$ values. The energy for this global mode is supplied by 3-D mechanisms – including large-scale Mach stem oscillations, bow shock dynamics and vortical structures – which become synchronised to feed the instability. Secondary spectral features confirm that boundary-layer turbulence persists in a subordinate role, establishing a hierarchical energy distribution where the global mode dominates the flow dynamics under choked conditions.

Figure 16.

The PSD of pressure perturbations at (x = 14 cm, y = 4.2 cm) for the 3-D simulations. The transition from broadband turbulence (J = 1.5, 3.0) to a dominant, low-frequency global mode ( $\textit{St}\approx 0.0115{-}0.016$ ) under choked conditions (J = 5.87) is evident. The lower panel (log scale) highlights both the dominant peak and broadband background; upper panels (linear scale) show mode energy and temporal evolution.

The characteristic Strouhal scaling ( $\textit{St}\approx 0.014$ ) observed in the 3-D simulation suggests the global mode frequency is governed by a system-scale time scale, consistent with a feedback mechanism between the shock system and downstream flow separation. However, the specific convective and acoustic speeds – and therefore the precise delay time scale – are expected to differ significantly from those of the 2-D case due to the presence of spanwise vortical structures, 3-D shock topology and additional dissipation pathways unique to the confined 3-D geometry.

Comparative analysis across J values reveals a critical transition: for the unchoked cases, the PSD exhibits broadband turbulence without dominant peaks at the low-frequency spectrum, whereas choked flows develop a sharp spectral signature. The persistence of this mode across the choked flow field confirms its role as the primary energy carrier, overwhelming finer-scale turbulent structures that persist as background noise.

The Strouhal scaling provides a quantitative criterion for predicting shock-dominated regimes in high-speed flows, with direct relevance to scramjet unstart prediction and control. The spectral dominance of shock oscillations at high J suggests that traditional turbulence models may require augmentation with shock motion physics in choked flow simulations. These findings establish PSD analysis as a diagnostic tool for identifying flow regime transitions in compressible systems, linking spectral peaks to specific instability mechanisms.

Magnitude-squared bicoherence analysis of the 3-D choked case ( $J=5.87$ ) reveals a significantly intensified nonlinear sustenance mechanism compared with the 2-D case. For a probe located at (14, 4.2, 0) cm, the 600 Hz global mode emerges within a region of strong nonlinear coupling, characterised by two distinct peaks at approximately 360 and 600 Hz (figure 17 a). This was confirmed by the nonlinear coupling region with highly strong nonlinear interactions shown in figure 17(b). The magnitude-squared bicoherence statistical analysis (figure 17 c) demonstrates that more than $15\,\%$ of frequency pairs exhibit strong nonlinear coupling ( ${\textit{bic}}^2 \gt 0.6$ ), substantially higher than the $4.4\,\%$ observed in the 2-D simulation. This elevated nonlinear activity indicates that 3-D confinement enhances and broadens the energy transfer pathways available to sustain the global mode, with the 360 Hz component serving as a primary feeder to the 600 Hz oscillation through intensified triad interactions.

Figure 17.

Magnitude-squared bicoherence analysis at a probe in the 3-D choked flow (J = 5.87). (a) Contour plot reveals intensified nonlinear interactions. (b) Cross-section at f/2. (c) Statistical distribution shows a higher percentage of strong nonlinear coupling compared with the 2-D case, indicating enhanced sustenance mechanisms for the global mode.

It is worth noting that the $\phi _p - \phi _u$ analysis and the DMD spatial phase measure fundamentally different aspects of the wave field – the former diagnoses the local energy-transport character of the waves, while the latter reveals the spatial coherence of the dominant mode. The near-zero $\phi _p - \phi _u$ values and the near-uniform DMD spatial phase together suggest a wave structure that is neither purely travelling nor purely standing, but rather a superposition of bidirectional waves phase-locked by the global instability.

7.2. Dynamic mode decomposition

Dynamic mode decomposition is a data-driven modal decomposition technique that extracts spatiotemporal coherent structures directly from time-resolved flow snapshots, each associated with a single oscillation frequency, and spatial structure (Schmid Reference Schmid2010). Unlike proper orthogonal decomposition, which ranks modes by energy content, DMD provides spectral information analogous to a global stability analysis without requiring knowledge of the governing equations.

Let $\boldsymbol{x}_n \in \mathbb{R}^{M}$ denote the vectorised schlieren of the $n$ th snapshot within the analysis window at the central plane, as schlieren cleanly exposes the major flow structure oscillations. The snapshot matrices

(7.3) \begin{equation} \boldsymbol{X} = \bigl [\boldsymbol{x}_1,\, \boldsymbol{x}_2,\, \ldots ,\, \boldsymbol{x}_{N-1}\bigr ], \quad \boldsymbol{X}' = \bigl [\boldsymbol{x}_2,\, \boldsymbol{x}_3,\, \ldots ,\, \boldsymbol{x}_{N}\bigr ], \end{equation}

are constructed from the mean-subtracted fluctuation field $\boldsymbol{x}_n' = \boldsymbol{x}_n - \bar {\boldsymbol{x}}$ , where $\bar {\boldsymbol{x}}$ is the time-averaged dataset. The DMD method assumes a linear map $\boldsymbol{x}_{n+1} = \boldsymbol{A}\,\boldsymbol{x}_n$ and seeks the best-fit operator $\boldsymbol{A}$ in the least-squares sense: $\boldsymbol{A} = \boldsymbol{X}'\boldsymbol{X}^{+}$ .

Following the exact DMD algorithm of Tu et al. (Reference Tu, Rowley, Luchtenburg, Brunton and Kutz2014), the high-dimensional operator $\boldsymbol{A}$ is never formed explicitly. Instead, the economy SVD of the data matrix,

(7.4) \begin{equation} \boldsymbol{X} = \boldsymbol{U}\boldsymbol{\varSigma }\boldsymbol{V}^{*}, \end{equation}

is truncated to rank $r$ to suppress noise, yielding the reduced operator

(7.5) \begin{equation} \tilde {\boldsymbol{A}} = \boldsymbol{U}_{\!r}^{*}\boldsymbol{X}'\boldsymbol{V}_{\!r}\boldsymbol{\varSigma }_r^{-1} \in \mathbb{R}^{r \times r}. \end{equation}

The eigendecomposition $\tilde {\boldsymbol{A}}\boldsymbol{W} = \boldsymbol{W}\boldsymbol{\varLambda }$ yields discrete eigenvalues $\lambda _k$ and exact DMD modes

(7.6) \begin{equation} \boldsymbol{\phi }_k = \frac {\boldsymbol{X}'\boldsymbol{V}_{\!r}\boldsymbol{\varSigma }_r^{-1}\boldsymbol{w}_k} {\lambda _k}, \end{equation}

where $\boldsymbol{w}_k$ is the $k$ th column of $\boldsymbol{W}$ . Each eigenvalue maps to a continuous-time frequency and growth rate via

(7.7) \begin{equation} f_k = \frac {\mathrm{Im}(\ln \lambda _k)}{2\pi \,\Delta t}, \quad \sigma _k = \frac {Re}(\ln \lambda _k){\Delta t}, \end{equation}

where $\Delta t = 3.69$ ${\unicode[Arial]{x03BC}}\mathrm{s}$ is the inter-snapshot interval.

Standard exact DMD eigenvalue estimates are biased towards the interior of the unit circle in the presence of measurement noise (Dawson et al. Reference Dawson, Hemati, Williams and Rowley2016). To mitigate this, the forward–backward DMD algorithm is employed: the forward operator $\boldsymbol{A}_{{fwd}}$ and the inverse of the backward operator $\boldsymbol{A}_{{bwd}}$ are computed independently, and the debiased eigenvalues are taken as the element-wise geometric mean:

(7.8) \begin{equation} \lambda _k^{fb} = \sqrt {\lambda _k^{{fwd}} \boldsymbol{\cdot }\left (\lambda _k^{{bwd}}\right )^{-1}}, \end{equation}

which pulls estimates towards the unit circle and reduces the systematic underestimation of $|\lambda _k|$ that would otherwise cause all modes to appear artificially decaying.

The complex amplitude of each mode is obtained by projecting the initial snapshot onto the DMD basis: $\boldsymbol{b} = \boldsymbol{\varPhi }^{+}\boldsymbol{x}_1'$ , where $\boldsymbol{\varPhi } = [\boldsymbol{\phi }_1, \ldots , \boldsymbol{\phi }_r]$ . Modes are ranked by their amplitude-based energy fraction (Jovanović et al. Reference Jovanović, Schmid and Nichols2014):

(7.9) \begin{equation} E_k = \frac {|b_k|^2}{\displaystyle \sum _{j:\,f_{\!j}\gt 0} |b_{\!j}|^2}, \end{equation}

computed exclusively over positive-frequency modes to avoid double-counting complex-conjugate pairs. This metric is preferred over the singular-value spectrum of $\boldsymbol{X}$ , which measures snapshot variance and carries no direct correspondence to individual DMD mode contributions.

Figure 18.

Dynamic mode decomposition of numerical schlieren snapshots on the central plane for the choked case (J = 5.87): (a) amplitude spectrum showing dominant mode at 630 Hz of the first 20 modes, (b) spatial magnitude of the dominant mode and (c) spatial phase showing coherent domain-wide oscillation with $0^\circ$ and $\pm 180^\circ$ phase regions.

The analysis is applied to $N = 1354$ schlieren snapshots spanning a physical duration of 5 ms at a frame-equivalent sampling frequency of 270.8 kHz. The spatial domain is restricted to the window which isolates the refined region while excluding the far field. The rank truncation is set to $r = 200$ , retaining sufficient degrees of freedom to represent the broadband schlieren content while suppressing low-energy noise modes.

While the initial DMD amplitudes ( $|b_k|$ ) of the 630 and 840 Hz modes are comparable (figure 18 a), the time-integrated energy metric $|b_k|^2/(1-|\lambda _k|^2)$ , which accounts for both amplitude and temporal persistence, reveals that the 630 Hz mode contributes $56.9\,\%$ of the total dynamical energy compared to $35.5\,\%$ for the 840 Hz mode. This confirms the 630 Hz mode as the dominant dynamical structure, consistent with its identification as the global mode through independent probe-based spectral analysis. This clearly identifies the 630 Hz mode as the dominant and most persistent coherent structure, in excellent agreement with the spectral peak at 600 Hz obtained from independent probe-based PSD analysis (figure 16).

Figure 18(b) displays the spatial magnitude of the dominant 630 Hz mode. The magnitude is not uniform; it is concentrated in major flow structures: the shock train, the core of the CVP, the separation shock, the barrel shock and downstream the leeward region of the injector. High magnitude near the shock train indicates that the mode strongly modulates the oscillatory motion of the Mach stem and the reflected shocks. Equally high magnitude inside the CVP core shows that the vortical structures are not passive – they actively participate in the global oscillation, acting as a distributed energy source which will be discussed in future work. The lower-wall separation zone, where the flame anchors, also exhibits significant modal activity, linking the instability to the combustion region. The spatial extent of the mode spans the entire channel height and a streamwise distance, confirming its system-scale, global character. No single local feature can account for this mode; it is a genuine coupled instability of the whole confined supersonic reacting flow.

The spatial phase of the same mode is shown in figure 18(c). The phase distribution is not random but exhibits a coherent pattern: large upstream and downstream regions oscillate nearly in phase (red), while a narrow band crossing the shock train and the CVP core shows a rapid transition to the opposite phase (blue). This phase shift of approximately $180^\circ$ indicates that the shock train and the CVP act as a source or a sink of the oscillation. In physical terms, when the pressure in the shock-train region increases, the pressure in the CVP core decreases in a nearly anti-correlated manner – a hallmark of a global acoustic–convective feedback loop. The DMD phase map therefore provides a direct visual confirmation of the bidirectional wave mechanism that orchestrates choking, complementing the quantitative phase-speed measurements. Together, the frequency ranking, the spatial magnitude and the phase pattern leave no doubt that the 630 Hz mode is the shock-driven global instability that governs the thermal choking transition.

7.3. Shock–boundary layer interaction

The onset of choking and subsequent transition is further elucidated through the evolution of the shock train, which significantly modulates the effective flow area. As depicted in the centreline plane (Z = 0; figure 19), the shock train manifests as a series of oblique and normal (Mach stem) shocks, progressively compressing the core flow and thickening the boundary layer downstream. This compression reduces the effective cross-sectional area, exacerbating flow blockage and sustaining the choked state. The Mach number contours in figure 19 reveal strong SBLIs along the centreline, while figure 20 highlights their 3-D nature, with lateral separation bubbles and transverse pressure gradients further distorting the flow field. The coupling between the shock train and boundary-layer displacement results in a feedback loop: the effective area reduction amplifies local Mach numbers upstream, which in turn strengthens the shock train. This self-sustaining mechanism explains the hysteresis observed during choking onset and the abrupt transitions between flow regimes.

Figure 19.

Mach number contours on the central plane (Z = 0) at t = 1.5 ms for (a) J = 1.5, (b) J = 3.0 and (c) J = 5.87. The panels show the shock-train structure and the intensification of SBLIs leading to flow separation and choking at high J.

Figure 20.

Mach number contours on a lateral plane (y = 1.5 cm) at t = 1.5 ms for (a) J = 1.5, (b) J = 3.0 and (c) J = 5.87. The 3-D nature of the flow separation, vortical structures and the interaction with the jet plume is highlighted.

The spatial and temporal evolution of the Mach stem foot is critically governed by the momentum flux ratio (J), revealing distinct regimes of shock stabilisation. For subcritical conditions (J = 1.5 and 3), the Mach stem forms at progressively upstream locations (x = 15 and 13 cm, respectively), consistent with enhanced stabilisation driven by increased momentum flux. In contrast, under choked flow conditions (J = 5.87), the system exhibits nonlinear dynamics: the Mach stem initially stabilises farther upstream (x = 11–12 cm) but subsequently undergoes temporal oscillations due to downstream pressure accumulation. This unsteadiness is mediated by a propagating normal shockwave, which perturbs the local flow topology and couples with the Mach stem’s inherent instability. The results highlight a dual role of momentum flux. This competition between stabilisation and accumulation mechanisms yields a time-dependent shock configuration, underscoring the sensitivity of supersonic interactions to flow confinement and compressibility effects.

The boundary-layer response to increasing J exhibits distinct anisotropy between the lateral and transverse directions, as visualised in figures 19 and 20. Upstream of the interaction, the boundary layer thickens progressively with J, but this growth is abruptly modulated by the incident bow shock. While the subcritical cases (J = 1.5 and 3) show comparable boundary-layer thicknesses in the spanwise direction, the higher-J case triggers earlier flow separation upstream, reflecting stronger adverse pressure gradients. This trend intensifies under choked conditions (J = 5.87), where the boundary layer undergoes significant thickening, indicative of enhanced separation dynamics. The asymmetry between lateral and transverse growth further underscores the 3-D nature of shock-induced separation, with the transverse direction exhibiting greater sensitivity to momentum flux variations.

Recent investigations of flow choking mechanisms in supersonic combustors (Jin et al. Reference Jin, Cai, Hong, Zhang and Liang2024) reveal that the spatial distribution of heat release fundamentally alters the shock structure. Figure 20 demonstrates how the developing shear layer modulates the shock-train dynamics. The transverse fuel injection creates a jet plume that behaves as a flow obstruction, generating a characteristic bow shock structure. Downstream of this interaction, exothermic reactions increase both static pressure and temperature, thereby reducing the local Mach number through the thermodynamic relationship $M \propto T^{-0.5}$ .

The flow-field complexity escalates with higher momentum flux ratios, which amplify vortical structures including the horseshoe vortex system. These vortices interact strongly with incident shockwaves, progressively constricting the effective flow area. Concurrently, the interaction of shear layer and boundary layer promotes the growth of separation bubbles (figures 5 c and 7), with the ensuing flow instability introducing additional oscillation modes and chaotic structures (figure 19). This cascade of events ultimately triggers the transition to large-scale flow separation (figures 19 and 20).