## 1. Introduction

Detonation waves are supersonic combustion waves that can sustain themselves in a premixed combustible mixture (Zhang, Ng & Lee Reference Zhang, Ng and Lee2012; Zhang *et al.* Reference Zhang, Mehrjoo, Ng, Lee and Bai2014; Shi, Uy & Wen Reference Shi, Uy and Wen2020). Due to their self-ignition, fast energy release and high thermal cycle efficiency, detonation phenomena as a combustion mode have attracted increasing attention in the development of hypersonic propulsion systems (Kailasanath Reference Kailasanath2000; Wolański Reference Wolański2013; Lu & Braun Reference Lu and Braun2014). Among the different kinds of detonation-based engines (Roy *et al.* Reference Roy, Frolov, Borisov and Netzer2004; Higgins Reference Higgins2006; Ma *et al.* Reference Ma, Luan, Xia, Wang, Zhang, Yao and Wang2020), the oblique detonation engine (ODE) is expected to have more applications in hypersonic air-breathing propulsion at high flight Mach numbers because of the additional stabilization and reduced compression losses of oblique detonation waves (ODWs) (Qin & Zhang Reference Qin and Zhang2018; Ren *et al.* Reference Ren, Wang, Xiang and Zheng2019; Jiang *et al.* Reference Jiang, Zhang, Liu, Wang and Luo2021). Hence, a comprehensive understanding of ODWs is fundamental for their successful utility in ODEs, and substantial progress in understanding ODWs has been achieved in recent decades. For example, Li, Kailasanath & Oran (Reference Li, Kailasanath and Oran1994) conducted numerical simulations in which wedge-induced ODWs were found to be initiated from a non-reactive oblique shock wave (OSW), which was confirmed by subsequent experimental observations (Viguier *et al.* Reference Viguier, da Silva, Desbordes and Deshaies1996; Lin, Zhang & Zhou Reference Lin, Zhang and Zhou2007; Han Reference Han2013). The detailed initiation structures of ODWs were further classified into two patterns, namely smooth and abrupt OSW–ODW transition patterns (Teng & Jiang Reference Teng and Jiang2012; Miao *et al.* Reference Miao, Zhou, Liu and Cai2018; Teng *et al.* Reference Teng, Tian, Zhang, Zhou and Ng2021). With regard to the initiation mechanisms of ODWs, their initiation can be controlled by chemical kinetics or wave dynamics depending on different free-stream parameters and mixture reactivities (Teng, Ng & Jiang Reference Teng, Ng and Jiang2017; Zhang *et al.* Reference Zhang, Fang, Ng and Teng2019). Moreover, similar to normal detonation waves (NDWs), ODW surfaces are always inherently unstable and composed of multiscale cell-like structures and moving transverse waves (Choi *et al.* Reference Choi, Kim, Jeung, Ma and Yang2007; Han, Wang & Law Reference Han, Wang and Law2019; Yang *et al.* Reference Yang, Teng, Ng and Jiang2019*b*).

Among the various key techniques required for the development of an air-breathing detonation-based engine (Schwartzentruber, Sislian & Parent Reference Schwartzentruber, Sislian and Parent2005; Alexander, Sislian & Parent Reference Alexander, Sislian and Parent2006; Wolański Reference Wolański2013; Zhang *et al.* Reference Zhang, Ma, Zhang, Han, Liu and Jiang2020*b*, Reference Zhang, Wen, Yuan, Liu, Han, Wang and Jiang2022), detonation-stabilization control is imperative to the engine's stable operations; nevertheless, achieving such control remains challenging because of the fast response of detonation waves to the incoming flow conditions and the extremely high relative propagation speed (thousands of metres per second). To stabilize an NDW in an engine combustor, the key point is to match the detonation propagation speed (always overdriven) with the incoming flow speed. If the propagation speed of the NDW is greater than the flow speed in the combustor, the NDW propagates upstream; conversely, the NDW attenuates downstream if its propagation speed is lower than that of the incoming flow. To dynamically stabilize the detonation in the engine combustor, Cai *et al.* (Reference Cai, Liang, Deiterding and Lin2016, Reference Cai, Deiterding, Liang, Sun and Dong2019, Reference Cai, Deiterding, Liang, Mahmoudi and Sun2021) employed a cavity as a flame holder to support an NDW in the supersonic incoming flow and used an expansion wall to attenuate the detonation wave downstream; the opposing effects of the cavity and the expansion wall on the propagation of the detonation wave compete and ultimately equilibrate with each other, which facilitates the stabilization of the NDW in the combustor. Moreover, Cai *et al.* (Reference Cai, Deiterding, Liang, Sun and Dong2020) further attempted to control the propagation of an NDW in a straight channel by setting a series of suction slots behind the detonation front; their numerical results suggested that the dynamic stabilization of the detonation wave can be achieved.

Many experiments (Sterling *et al.* Reference Sterling, Cummings, Ghorbanian, Pratt, Sobota, Brock, Brown, Segall and Debarber1998; Verreault & Higgins Reference Verreault and Higgins2011; Maeda, Kasahara & Matsuo Reference Maeda, Kasahara and Matsuo2012; Maeda *et al.* Reference Maeda, Sumiya, Kasahara and Matsuo2013) and simulations (Teng *et al.* Reference Teng, Ng and Jiang2017; Han *et al.* Reference Han, Wang and Law2019; Yang, Ng & Teng Reference Yang, Ng and Teng2019*a*) have previously been conducted on the stabilization of ODWs, suggesting that it is easy to stabilize an ODW over isolated wedges or cones, which is one of the major advantages of ODEs. However, the stabilization of an ODW induced by an isolated wedge or cone is only a typical external flow problem, whereas the stabilization of an ODW in an ODE combustor is an internal flow problem because of the geometric constraints of the combustor's internal walls. In particular, the additional factors involved in internal flows, such as the reflection of OSWs/ODWs (Wang *et al.* Reference Wang, Teng, Yang and Ng2020*a*,Reference Wang, Zhang, Yang and Teng*b*; Wang, Yang & Teng Reference Wang, Yang and Teng2021), shock/detonation–shock/detonation interactions (Xiang *et al.* Reference Xiang, Li, Zhang, Xie and Zhang2021*a*,Reference Xiang, Zhang, Gao, Li and Huang*b*), shock/detonation–boundary layer interactions (Cai *et al.* Reference Cai, Deiterding, Liang and Mahmoudi2017, Reference Cai, Deiterding, Liang, Sun and Mahmoudi2018*a*,Reference Cai, Liang, Deiterding, Mahmoudi and Sun*b*) and shock/detonation-induced boundary layer separation (Ess, Sislian & Allen Reference Ess, Sislian and Allen2005; Miao *et al.* Reference Miao, Xu, Song and Yu2020), increase the flow complexity in a space-confined combustor and cannot be overlooked in the stabilization of ODWs. These fundamental differences in detonation stabilization between spatially confined and unconfined flow regimes were also pointed out by Higgins (Reference Higgins1997) through experiments, where hypervelocity blunt projectiles were fired into detonation chambers with significantly different chamber diameters and various detonation propagation phenomena were observed. In early experimental investigations of ram accelerators (e.g. Knowlen *et al.* Reference Knowlen, Higgins, Bruckner and Hertzberg1995; Higgins, Knowlen & Bruckner Reference Higgins, Knowlen and Bruckner1998; Yatsufusa & Taki Reference Yatsufusa and Taki2002), which share geometries similar to those of ODE combustors, it was demonstrated that boundary layer separation induced by shock waves and/or combustion is always responsible for the upstream motion of shock waves over the projectiles, resulting in the unstart phenomena and, consequently, the failure of acceleration and operation. Moreover, the effects of multiple shock reflections on the upstream propagation of ODWs over projectiles in superdetonative ram accelerators and even unstart of the systems have also been noted in experiments (Seiler *et al.* Reference Seiler, Patz, Smeets and Srulijes1998, Reference Seiler, Patz, Smeets and Srulijes2000). Therefore, despite the inherent stability of ODWs in an unconfined space, their stabilization control is always necessary in an ODE combustor; to this end, it is valuable to comprehensively understand the effects of the aforementioned factors on the stabilization of ODWs and the relevant mechanisms.

The effects of boundary layer separation on the failures of ramjets and scramjets have also been well documented in the literature (Curran & Murthy Reference Curran and Murthy2001; Im & Do Reference Im and Do2018; Seleznev, Surzhikov & Shang Reference Seleznev, Surzhikov and Shang2019). Accordingly, various techniques, such as boundary layer blowing/bleeding and mass flow spillage, have been proposed and successfully applied (Chang *et al.* Reference Chang, Li, Xu, Bao and Yu2017; Urzay Reference Urzay2018). Recently, these techniques have also been employed to overcome separation-induced destabilization of ODWs in ODE combustors. For example, Alexander & Sislian (Reference Alexander and Sislian2008) and Wang & Sislian (Reference Wang and Sislian2010) employed an air jet to blow down the boundary layer upstream of the ODW reflection locus to suppress the detonation-induced separation of the boundary layer. Zhang *et al.* (Reference Zhang, Ma, Zhang, Han, Liu and Jiang2020*b*, Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*, Reference Zhang, Wen, Yuan, Liu, Han, Wang and Jiang2022) adopted a floor bleed structure at the combustor's entrance to spill out the incoming boundary layer from the inlet; as a consequence, the continuous upstream motion of the boundary layer separation bubble, also induced by the reflection of the ODW, ceases at the entrance because of the interruption of the boundary layer. With boundary layer separation being addressed properly, destabilization of the flow field in the ODE combustor may still occur under the influences of the aforementioned inviscid factors (i.e. shock/detonation reflections, shock/detonation–shock/detonation interactions). In other words, the inviscid stabilization of ODWs, which is related to the geometric structure of the combustor, appears as a precondition of overall stabilization and must be considered prior to adopting other stabilization-control measures in the design of an ODE combustor.

Wave reflection is a major feature in the inviscid ODW flow field in a space-confined ODE combustor (Lu, Fan & Wilson Reference Lu, Fan and Wilson2006; Wang *et al.* Reference Wang, Yang and Teng2021; Xiang *et al.* Reference Xiang, Zhang, Gao, Li and Huang2021*b*), and its effects on ODW stabilization and the relevant mechanisms are important for the initial geometric design of the engine combustor. Figure 1 schematically shows two typical configurations of ODEs presented in previous literature (Schwartzentruber *et al.* Reference Schwartzentruber, Sislian and Parent2005; Alexander *et al.* Reference Alexander, Sislian and Parent2006; Wolański Reference Wolański2013), namely the external-compression ODE and the mixed-compression ODE. These two types of ODEs have an obviously similar flow pattern of ODW reflection in the combustor; that is, the ODW reflects over an internal wall with an expansion corner, and the ideal reflection locus is precisely at the expansion corner. However, due to the general unsteadiness of the high-speed fuel–air mixture flowing into the combustor and/or the potential off-design operation of the engine (Yang *et al.* Reference Yang, Ng and Teng2019*a*, Reference Yang, Ng and Teng2021; Ren, Wang & Zheng Reference Ren, Wang and Zheng2021), it is impossible to ensure that the ODW always reflects exactly at the expansion corner, implying that the ODW always reflects before or behind the expansion corner. Therefore, it is imperative to investigate the effects of the position of the reflection point relative to the expansion corner and other geometric parameters on the ODW flow structure in the combustor and to analyse the relevant stabilization characteristics. Because of the release of heat, a Mach reflection pattern can appear easily when the ODW reflects off a wall, and a Mach stem (MS; which is, in fact, an overdriven NDW) forms as a result. Due to the existence of a subsonic zone behind the MS, downstream disturbances may propagate upstream to the detonation front and consequently destabilize the ODW reflection pattern, eventually causing the ODE to fail. Accordingly, this paper mainly investigates the relevant flow structures, stabilization characteristics and inherent mechanisms of a Mach reflection occurring more easily when the ODW reflects off the wall before the expansion corner.

The remainder of this paper is organized as follows. The numerical details, including the geometry of the simulated ODE combustor, the numerical methods and the cases considered herein, are briefly introduced in § 2. Next, the flow fields and the relevant stabilization characteristics of ODWs that reflect precisely at the expansion corner at different flight Mach numbers are discussed in § 3.1. Then, the flow structures of the stabilized Mach reflection of an ODW reflecting before but close to the expansion corner and the evolutions of the relevant MS are presented in § 3.2. Thereafter, the evolutions of the destabilized Mach reflection of an ODW reflecting far upstream of the expansion corner and the relevant destabilization mechanism are analysed in § 3.3. In § 3.4, the inherent destabilization phenomenon of the Mach reflection of an ODW at a low flight Mach number is investigated, followed by a discussion of its relevant mechanism in § 3.5 with a formulation providing the stabilized locations of the MS. Finally, the conclusions are presented in § 4.

## 2. Physical and mathematical models

Since the flow patterns of ODW reflection are similar in the combustors of both configurations of ODEs (figure 1), the combustor of an external-compression ODE, whose detailed geometry is shown in figure 2, is taken as an example in this study. As shown in figure 2, the combustor is contained within points $\overline {\textrm{DCBOE}} $, where point C is the leading tip of the cowl and point O is the expansion corner with an expansion angle of 25°. $\overline {\textrm{CD}} $ or $\overline {\textrm{CG}} $ represents the engine's cowl wall. $\overline {\textrm{BO}} $ is one part of the combustor's lower wall before the expansion corner with an inclination angle of 25°, while $\overline {\textrm{OE}} $ (which is parallel to the cowl wall) is the other part behind the expansion corner. $\overline {\textrm{EF}} $ is a part of the nozzle's wall with an expansion angle of 25°, while $\overline {\textrm{AB}} $ is a part of the inlet's wall sharing the same inclination angle as $\overline {\textrm{BO}} $. The height $(\overline {\textrm{DE}} )$ of the combustor's straight channel is *H* = 6 cm, while its length $(\overline {\textrm{OE}} )$ is *L* ^{′} = 20 cm. The total length of the combustor *L* depends on the designed reflection locus of the ODW (i.e. point R), which varies among the different cases considered in this study. A similar ODE combustor geometry has been adopted for previous studies by the authors (Zhang *et al.* Reference Zhang, Ma, Zhang, Han, Liu and Jiang2020*b*, Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*).

The effects of ODW reflection on the stabilization/destabilization characteristics of ODWs as well as the inviscid destabilization mechanisms are emphasized in this study; hence, viscous effects are neglected. However, the effects of viscosity on the transition between regular and Mach reflections should be noted. It has been demonstrated that shock-induced boundary layer separation has significant effects on the shifting of the dual-solution regime of inert shock-wave reflections (Tao, Fan & Zhao Reference Tao, Fan and Zhao2014; Matheis & Hickel Reference Matheis and Hickel2015; Grossman & Bruce Reference Grossman and Bruce2018; Xue, Wang & Cheng Reference Xue, Wang and Cheng2021). That is, it delays the transition from regular to Mach reflection, and its effect on the transition from Mach to regular reflection is inverse (i.e. promotion). For ODW reflection transitions, although similar viscous effects may be expected, the specific transition criteria are not the main focus of this study. Therefore, the complex effects of boundary layer separation on ODW reflection transitions are not considered but, instead, will be further studied in the future. Without considering viscous effects, two-dimensional time-dependent multispecies reactive Euler equations (Fang *et al.* Reference Fang, Zhang, Hu and Deng2019; Teng, Liu & Zhang Reference Teng, Liu and Zhang2020) are solved. Jachimowski's H_{2} combustion mechanism (Wilson & Maccormack Reference Wilson and Maccormack1992), which involves 19 reversible elementary reactions among 9 species (H_{2}, H, O_{2}, O, OH, HO_{2}, H_{2}O_{2}, H_{2}O and N_{2}), is adopted to model the combustion rates. This detailed chemical reaction mechanism has been demonstrated to be capable of providing experimentally consistent laminar flame speeds, adiabatic flame temperatures and ignition delay times (Jachimowski Reference Jachimowski1988) and has been widely applied in simulations of scramjets (Fureby *et al.* Reference Fureby, Chapuis, Fedina and Karl2011; Chapuis *et al.* Reference Chapuis, Fedina, Fureby, Hannemann, Karl and Schramm2013), shock-induced combustion (Choi, Jeung & Yoon Reference Choi, Jeung and Yoon2000; Teng *et al.* Reference Teng, Liu and Zhang2020) and detonation (Choi, Shin & Jeung Reference Choi, Shin and Jeung2009; Zhang *et al.* Reference Zhang, Ma, Zhang, Han, Liu and Jiang2020*b*, Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*).

It should be noted that the temperature behind the detonation wave of a H_{2}–air mixture can reach 3000 K or higher; hence, vibration excitations of some involved polyatomic species (i.e. O_{2}, N_{2}, H_{2}O, HO_{2} and H_{2}O_{2}) with relatively low vibrational characteristic temperatures become noticeable. To include the excited vibrational energies and other energy modes in the internal energy of each species, in this study, the species thermodynamic database of NASA (McBride, Gordon & Reno Reference Mcbride, Gordon and Reno1993) is employed to evaluate the equilibrium thermodynamic properties of each species by a piecewise fourth-order polynomial of temperature. In recent studies (Shi *et al.* Reference Shi, Shen, Zhang, Zhang and Wen2017, Reference Shi, Uy and Wen2020; Uy *et al.* Reference Uy, Shi, Hao and Wen2020), thermal non-equilibrium has been demonstrated to have non-negligible effects on detonation instabilities and the resulting detonation cellular structures in a multidimensional space. However, the aforementioned detonation dynamic features are not the main focus of this study. Moreover, a comparison of the flow fields in a similar ODE combustor under thermal equilibrium and non-equilibrium conditions (conducted in the authors’ previous study: Zhang *et al.* (Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*)) showed that the overall flow structures or their stabilization characteristics are not significantly affected by thermal non-equilibrium effects. Therefore, thermal non-equilibrium effects are neglected in the present modelling; that is, all simulation results discussed in this study are under the thermal equilibrium assumption.

A quadrilateral grid-based finite-volume method (Chakravarthy Reference Chakravarthy1999) is employed to numerically solve the governing equations. A second-order upwind total variation diminishing scheme is used to spatially discretize the equations, and the nonlinear Harten–Lax–van Leer contact approximate Riemann solver (Toro Reference Toro2013) is adopted to evaluate the interface fluxes to satisfy the entropy and positivity conditions simultaneously. Additionally, to suppress spurious oscillations near flow discontinuities, a minmod-type limiter is used. In the time direction, the equations are explicitly integrated by the fourth-order Runge–Kutta method, and an operator-splitting method (Ropp & Shadid Reference Ropp and Shadid2009) is employed to treat the chemical source terms, which always cause stiffness problems originating from the mismatch of the time scales between chemical reactions and flows.

In the present study, numerical simulations are carried out with different flight Mach numbers and different designed ODW reflection locations before or precisely at the expansion corner, as summarized in table 1. Notably, in classic shock reflection configurations in which an OSW generated by a finite wedge impinges on a flat surface, the height of the MS depends on the height of the flow path (i.e. the vertical distance from the wedge tail to the flat surface) (Li & Ben-Dor Reference Li and Ben-Dor1997; Ben-Dor Reference Ben-Dor2007; Mouton & Hornung Reference Mouton and Hornung2007). According to shock reflection theory (Ben-Dor Reference Ben-Dor2007), the formation of a sonic throat under the interaction between the expansion fan and the slip line (SL) is important to the stabilized location and, consequently, the stabilized height of the MS. In such a type of shock reflection configuration, the expansion fan emits from the tail of the finite wedge on the opposite side of the reflection locus; hence, the interaction between the expansion fan and the SL depends on the vertical distance from the wedge tail to the flat surface. As a result, the height of the MS in such configurations depends on the height of the flow path. Comparatively, in the ODW reflection configuration depicted in figure 2, the expansion fan emits from an expansion corner on the same side of the reflection locus. As a result, the stabilized height of the MS discussed in this study (if not influenced by downstream wave structures) does not depend on the height of the combustor's straight channel (i.e. *H*) but on the designed reflection location of the ODW (i.e. $\overline {\textrm{RO}} $). However, destabilization of ODW Mach reflection may still be caused by downstream wave structures that result from the multiple reflections of shock waves between the combustor's walls. Obviously, with the designed reflection location of ODW unchanged, the relative locations of these downstream wave structures to the expansion corner are determined by the channel's height. In other words, the stabilization/destabilization characteristics of ODWs in the combustor considered in this study depend on not only the designed reflection location but also the geometric scale of the combustor (i.e. the height of the combustor's straight channel, *H*). Hence, to present and discuss the results in a more general way, the non-dimensional parameter *ζ* is introduced in this study to describe the relative location of the designed ODW reflection point:

A value of *ζ* = 0 means that the ODW reflects exactly at the expansion corner (the dashed red line in figure 2), while *ζ* > 0 implies that the ODW reflects off the wall before the expansion corner, as shown by the solid red line in figure 2. Notably, without considering the chemical reaction characteristic length scales (i.e. in cold or chemical equilibrium flows), the non-dimensionalization of (2.1) is justifiable. The flow field in the combustor should depend on the non-dimensional parameter *ζ*, rather than the dimensional reflection location or channel height individually. In the present study, because chemical non-equilibrium effects are not dominant (i.e. near-equilibrium overall, to be discussed later in § 3.1), (2.1) still approximately holds, and the overall flow fields with the same *ζ* value but different channel heights should be similar.

Notably, the exact reflection locations of ODWs on the wall (i.e. $\overline {\textrm{RO}} $) and the subsequent *ζ* values are difficult to determine when Mach reflection occurs or the flow structures in the combustor become destabilized because an MS forms before the ODW reflecting off the wall, or the ODW is even non-stationary. To overcome this difficulty, a preinvestigation simulation of a pure ODW induced only by the cowl wall (i.e. $\overline {\textrm{CD}} $) is conducted for each flight Mach number without considering the combustor's lower wall (i.e. $\overline {\textrm{BOE}} $), and the corresponding reflection locus of the ODW (i.e. point R) is determined by its intersection point with the location at which $\overline {\textrm{BOE}} $ should be. Note again that the *ζ* values discussed in this study are generally smaller than one; hence, *ζ* values of different cases are expressed in fractional terms with a constant denominator of 60 (i.e. the channel height of 60 mm) in table 1 and the remainder of this paper, for intuitional comparison.

Four flight Mach numbers, namely 9, 9.5, 10 and 10.5, are taken into account, and the flight altitude is fixed at 30 km for the ODE. Since this study focuses on the effects of reflection on the stabilization characteristics of ODWs in a space-confined ODE combustor and the relevant inherent mechanisms of stabilization/destabilization, the complex processes of fuel (H_{2}) injection and its subsequent mixing with air in the ODE inlet are not considered here. Rather, following the simplification of the inflow of the ODE combustor in previous studies (Ren *et al.* Reference Ren, Wang, Xiang and Zheng2018; Fang *et al.* Reference Fang, Zhang, Hu and Deng2019; Xiang *et al.* Reference Xiang, Li, Sun and Chen2019), the free stream of the high-altitude atmosphere is assumed to be precompressed twice by two weak OSWs generated by the two 12.5° ramps in the inlet (see figure 1*a*), and the injection of H_{2} and its subsequent mixing with air are assumed to be completed downstream of these two OSWs and before entering the ODE combustor. Consequently, a uniformly premixed stoichiometric H_{2}–air inflow is assumed to enter the combustor parallel to the end of the inlet and the wall before the expansion corner (i.e. 25° with respect to the *x* direction), and this serves as the inflow condition for the simulations in this study, as shown in figure 2. The corresponding inflow parameters for the different flight Mach numbers are summarized in table 1. Note again that the vibrational relaxation time scales of H_{2}–air mixture in the inlet are significantly smaller than the corresponding flow residence time scale; hence, it is reasonable to set the inflow of the combustor (i.e. the outflow of the inlet) to be under thermal equilibrium.

For the other boundary conditions, the slip wall condition is set at all wall boundaries due to the assumed inviscid flow, and the supersonic outflow condition is set at all outflow boundaries. Quadrilateral grid cells are used in this study, and the characteristic grid size in the main combustion zone is set as 0.1 mm, which has been carefully validated by resolution tests to exclude the grid dependence of the simulation results. One of these grid resolution studies is presented in § 3.2 as an example.

## 3. Results and discussion

### 3.1. Stabilized reflection of ODWs with *ζ* = 0

The reflection patterns of ODWs in the combustor at different flight Mach numbers with the reflection point located precisely at the expansion corner (i.e. *ζ* = 0) are analysed first. The time-dependent simulation results suggest that the flow structures with *ζ* = 0 all remain stabilized in the combustor at the investigated Mach numbers (9, 9.5, 10 and 10.5). Taking *Ma* = 9 and 10 as examples, the stabilized flow structures are depicted in figure 3 by shadowgraphs, temperature contours, H_{2} mass fraction contours and pressure isolines. When the high-speed premixed H_{2}–air mixture flows into the combustor, an OSW first forms over the leading tip of the cowl, inducing combustion of the mixture downstream. Away from the cowl wall, the combustion zone approaches the OSW front and eventually interacts with the shock-wave front over a short distance. Consequently, the flame front is coupled tightly with the OSW front, and the shock-wave angle increases, implying the formation of an ODW. Across this ODW, the flow temperature abruptly increases to approximately 3000 K, and the H_{2} mass fraction drops to a rather low level; these phenomena are the other two typical features of a detonation wave. Then, the ODW reflects at the expansion corner within the combustor, forming a reflected shock wave (RfSW), which is weak due to the parallel geometry of the flow channel behind the expansion corner. Downstream in the combustor, the RfSW reflects repeatedly between the internal walls and further weakens gradually before it meets another expansion fan at the outlet of the combustor. Finally, the high-temperature combustion products exit the combustor and expand and accelerate in the nozzle.

A comparison of the flow structures between these two Mach numbers in figure 3 reveals that the initiation length under *Ma* = 10 is shorter than that under *Ma* = 9. A transverse wave (TW) connecting to the OSW–ODW transition point exists in the case of *Ma* = 9, implying an abrupt OSW–ODW transition pattern at this relatively low Mach number, while no obvious TW can be observed within the initiation zone in the case of *Ma* = 10, meaning the OSW–ODW transition pattern at this relatively high Mach number is smooth. These findings regarding the ODW initiation structure coincide with those found in previous numerical studies (Teng & Jiang Reference Teng and Jiang2012; Teng *et al.* Reference Teng, Ng and Jiang2017, Reference Teng, Tian, Zhang, Zhou and Ng2021). Notably, as compared to the geometric scales of the combustor, the OSW–ODW transition zones are rather small, implying that the overall flow fields in the combustor are in near-equilibrium states, although some non-equilibrium effects still occur locally. Further, a higher Mach number leads to a smaller oblique detonation angle, and thus the leading tip of the cowl occurs farther upstream from the expansion corner for *ζ* = 0 (with the same height of the combustor). Additionally, the RfSW angles and its further reflected shocks are smaller for higher Mach numbers. Moreover, the combustion temperature in the combustor at *Ma* = 10 is slightly higher than that at *Ma* = 9, which is attributed to the higher total temperature but the same deflection angle of the inflow entering the combustor.

### 3.2. Stabilized Mach reflection of ODWs with small *ζ* > 0

When the ideal reflection point of the ODW is located on the wall before the expansion corner (i.e. for *ζ* > 0), the Mach reflection pattern of the ODW coincides with the formation of an MS. For small *ζ* > 0, the flow field under *Ma* = 10 and *ζ* = 10/60 is taken as an example. Figure 4 depicts the evolution of the flow field in this case by shadowgraphs overlaid with sonic lines. To quantitatively exhibit the evolution of the formed MS and the stabilization characteristics of the Mach reflection of the ODW, the temporal evolutions of the location and corresponding motion speed of the pressure-jump point along the wall before the expansion corner are plotted in figure 5, where the location of the pressure-jump point is evaluated according to its distance from the expansion corner. Obviously, this pressure-jump point represents the location of the ODW reflection point for a regular reflection pattern or the location of the MS for a Mach reflection pattern. Therefore, figure 5 also presents the temporal evolutions of the formed MS when a Mach reflection pattern appears in the late period of the evolution of the flow field. In addition, points at different time instants corresponding to the shadowgraphs shown in figure 4 are also marked on the motion speed curve in figure 5 for convenience.

As indicated in figure 4(*a*), after the ODW has formed over the cowl's internal wall but before it becomes fully established, its reflection point is located on the wall behind the expansion corner, and no pressure jump occurs along the wall before the expansion corner. At this moment, the location of the pressure-jump point and its corresponding motion speed are simply set to zero (point ‘a’ in figure 5). Hence, when the ODW reflects precisely at the expansion corner later in its evolution (figure 4*b*), the motion speed of the pressure-jump point abruptly rises (point ‘b’ in figure 5). Next, a regular reflection pattern occurs over the wall before the expansion corner (figure 4*c*) because of the relatively small oblique detonation angle of the tail of the ODW during its evolution, and the reflection point moves upstream along the wall continuously with a gradually decreasing motion speed (point ‘c’ in figure 5). Then, the regular reflection pattern of the ODW transits to a Mach reflection pattern (figure 4*d*), and the motion speed of the pressure-jump point reaches another obvious peak (point ‘d’ in figure 5) during the transition process. As depicted in figures 4(*e*) and 4( *f*), an MS forms and grows thereafter while continuously moving upstream along the wall, which is again accompanied by a decreasing motion speed. In addition to the formation of the MS and RfSW connected to the triple point, an SL is also emitted from the triple point in the Mach reflection pattern of the ODW, thereby separating the flows across the ODW and MS. When the RfSW of the ODW reflects off the cowl wall, a secondary reflected shock wave (RfSW2) forms and reflects later over the SL, forming another reflected shock wave and, notably, a transmitted MS (tMS) between the SL and the lower wall of the combustor. As indicated in figure 5, after the formation of the Mach reflection pattern, the motion speed of the formed MS rapidly drops and then gradually decreases towards zero. After a sufficiently long period, the motion speed of the MS finally drops to zero, with the MS finally reaching its most upstream location. As a result, all the flow structures remain stabilized in the combustor, and the flow field reaches its steady state, as shown in figure 4(*g*).

The temperature and H_{2} mass fraction contours along with the pressure isolines of the final steady flow field in this case are exhibited in figure 6. As indicated in figures 4(*g*) and 6, in addition to the ODW, a high temperature reaching approximately 3500 K and a rather low H_{2} mass fraction are also featured just behind the MS, implying the tight coupling between the flame front and the MS. In other words, the MS is an overdriven NDW. According to Zhang *et al.* (Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*), two stabilized detonation combustion modes, namely the oblique detonation combustion mode and the normal detonation combustion mode, simultaneously exist in the ODE combustor in this case. Moreover, a subsonic zone capable of propagating pressure waves upstream exists behind the MS and extends to the expansion fan emitting from the expansion corner (figure 4*g*). However, no obvious downstream pressure wave acting on or propagating in the subsonic zone is observed, and hence the MS is not disturbed and remains stabilized at a certain position before the expansion corner.

To exclude the dependence of the simulation results on the grid size, grid resolution tests are carefully carried out by halving the sizes of the quadrilateral grid cells in both directions, leading to characteristic grid sizes of 0.05 mm in the main combustion zone for the fine grids and 0.1 mm for the coarse grids. As a result, the total number of fine grid cells is almost four times that of coarse grid cells. Taking the case of *Ma* = 10 and *ζ* = 10/60 as an example again, comparisons of the flow fields obtained with different grid sizes are shown in figure 6, while comparisons of the flow parameter distributions are depicted in figure 7. Differences in the flow structures are difficult to distinguish by using different grid sizes. Furthermore, the distribution curves of the various flow parameters obtained with different grid sizes almost overlap with each other, especially the locations of the jumps or drops of the flow parameters, implying that the simulations accurately capture the locations of key flow structures such as detonation fronts and shock waves. Therefore, a coarse grid with a characteristic grid size of 0.1 mm is sufficient to resolve the key flow features and is adopted in this study.

It should be noted that, although the overall flow structures concerned in this paper are shown to be independent of grid size approximately when a characteristic grid size of 0.1 mm is used, this grid size is still not fine enough to precisely predict the induction length of an ODW or capture the potential cellular structures on the detonation surface. The grid-size requirement for predicting the aforementioned detonation features generally needs to ensure that there are at least 20 grid points within the half reaction zone of the corresponding Chapman–Jouguet detonation wave (Teng, Jiang & Ng Reference Teng, Jiang and Ng2014; Shen & Parsani Reference Shen and Parsani2017; Shi *et al.* Reference Shi, Uy and Wen2020). For example, for the inflow conditions of *Ma* = 10 listed in table 1, the length of the half reaction zone is approximately 0.4 mm, implying a minimum grid-size requirement of 0.02 mm, which is significantly smaller than that used in the present study (0.1 mm). Fortunately, the precise prediction of the induction length of an ODW does not affect the overall flow structures in the combustor, and the cellular characteristics of the detonation surfaces are not the main focus of this study; hence, a much finer grid is not employed to save computational resources.

For the other cases with small positive *ζ* values (i.e. *ζ* = 2/60 and 5/60) at *Ma* = 10, the flow fields of the Mach reflection of the ODW are also stabilized, and similar flow structures can be observed, as shown in figure 8. In other words, an MS, an SL, an RfSW and a subsonic zone form near the reflection point of the ODW. Downstream of the ODW, the RfSW reflects off the cowl wall to form an RfSW2, and the RfSW2 further reflects over the SL, forming another reflected shock wave and a tMS between the SL and the lower wall of the combustor. As *ζ* increases, i.e. as the designed reflection point of the ODW moves upstream, the final stabilized location of the MS also moves upstream; additionally, the length of the MS increases, as does the size of the subsonic zone behind the MS. Further, the locations of the RfSW2 and tMS downstream of the subsonic zone move farther upstream and approach the subsonic zone with increasing *ζ*. For *Ma* = 9.5 or 10.5, the Mach reflection patterns of the ODW similarly become stabilized at small positive *ζ* values, as summarized in figure 9, and stabilized flow structures similar to those discussed in this section can be obtained.

### 3.3. Destabilized Mach reflection of ODWs with large *ζ* > 0

Again, *Ma* = 10 is taken as an example for discussion in this section. Different from the stable cases with small *ζ*, when *ζ* further increases to a larger value, for example *ζ* = 20/60 or 30/60, the flow field of the Mach reflection of the ODW becomes destabilized in the ODE combustor, as depicted in figure 9. Similar to figure 5, the temporal evolutions of the location of the pressure-jump point along the wall before the expansion corner and its corresponding motion speed in the cases of *ζ* = 20/60 and 30/60 at *Ma* = 10 are shown in figure 10. According to the above discussion in § 3.2, these curves represent the evolutions of the formed MS, where the MS forms at the time instant of the second peak in the motion speed curve. In the early period, the motion speed of the MS, which continuously moves upstream, is greater than zero but decreases, rapidly at first and then gradually. Accordingly, the upstream motion of the MS approaches zero and seemingly stops at a certain position, and the MS appears to be stabilized. However, as the flow field continues to evolve, at a certain time instant (approximately *t* = 0.9 ms for *ζ* = 20/60 or *t* = 0.7 ms for *ζ* = 30/60), the motion speed of the MS suddenly begins to increase, after which the MS accelerates while moving upstream. Finally, the MS moves upstream far from the expansion corner, indicating that the MS has moved out of the combustor, and the MS continues moving upstream with a very large speed. In other words, the Mach reflection of the ODW becomes destabilized in the combustor at a large value of *ζ*, which may lead to the engine becoming unstart. Therefore, to help design the geometry of the combustor, it is imperative to determine the mechanisms responsible for the destabilized Mach reflection of ODWs as well as the relevant influencing factors.

Taking *ζ* = 30/60 as an example, the evolution of the flow field of the destabilized Mach reflection of an ODW at *Ma* = 10 is illustrated in figure 11 by shadowgraphs along with sonic lines. The time instants corresponding to the exhibited flow fields are also marked along the motion speed curve of the pressure-jump point in figure 10 for convenience. After the Mach reflection pattern of the ODW is established (figure 11*a*), flow structures similar to the stabilized structures in § 3.2, including the MS, SL, RfSW, RfSW2, tMS and subsonic zone, also appear in the flow field of this destabilized case. In the early period of evolution of this destabilized Mach reflection of the ODW, the tMS is located relatively far downstream from the subsonic zone behind the MS (figure 11*a*). As time proceeds, the upstream motion of the MS drives the following flow structures (e.g. the RfSW, RfSW2 and tMS) upstream; consequently, the location of the tMS approaches the subsonic zone behind the MS (figure 11*b*). As indicated in figure 10(*a*), the flow structures tend to stabilize since the motion speed decreases to zero during this period. However, at approximately *t* = 0.585 ms (figure 11*c*), the tMS acts on the subsonic zone before the flow structures are fully stabilized. Thereafter, the pressure rise caused by the tMS rapidly propagates upstream through this subsonic zone, as shown in figures 11(*d*) and 11(*e*). Note that the symbol ‘tMS’ in these two figure panels is enclosed in parentheses, implying that the pressure rise caused by the tMS is no longer as sharp as that of a shock wave due to the diffusion effects of propagating through a subsonic zone.

The upstream propagation of the pressure rise caused by the tMS is also clearly visualized in figure 12, showing the evolution of the pressure distributions along the wall before the expansion corner (after the tMS acts on the subsonic zone). As shown in figure 12, the pressure rise caused by the tMS propagates upstream through the subsonic zone and finally arrives at the detonation front of the MS at approximately *t* = 0.713 ms (see also figure 11*f*), resulting in an increase in the pressure behind the MS and consequently the pressure ratio across the MS. Notably, the pressure ratio across an overdriven detonation wave must match its propagation speed relative to the incoming flow, and vice versa. Therefore, this increase in the pressure ratio across the MS results in an increase in its relative propagation speed; that is, the MS accelerates while moving upstream, as indicated by the sudden increase in the motion speed of the MS (point ‘f’) in figure 10(*a*). Then, the upstream motion of the MS enters a so-called ‘positive feedback’ period: the more upstream the MS travels, the closer the distance between the MS and RfSW2 (figure 11*g*), the higher the pressure ratio across the MS induced by the RfSW2 pressure rise and, consequently, the larger the upstream motion speed of the MS. As depicted in figure 10(*a*), the motion speed of the MS increases continuously and reaches a large value after the pressure rise produced by the tMS arrives at the detonation front. At approximately *t* = 1.345 ms (figure 11*h*), the MS reaches the leading tip of the cowl, and the ODW vanishes, resulting in only an NDW existing at the entrance of the combustor. Subsequently, the detonation front moves out of the combustor; instead, a bow detonation wave (BDW) forms and continues moving upstream, as shown in figure 11(*i*). In this destabilized case of the Mach reflection of an ODW (*Ma* = 10, *ζ* = 30/60), because the detonation front that serves as the main combustion organization structure in the ODE moves out of the combustor and thereafter continues to move upstream with a large motion speed, the combustor (and thus the engine) cannot work normally and may even fail.

For *ζ* = 20/60 at the same *Ma*, analyses of the instantaneous flow fields in combination with the evolutions of the MS location and its motion speed suggest that the sudden upstream acceleration of the MS during its apparent stabilization (figure 10*b*) can also be attributed to the action of the tMS on the subsonic zone behind the MS (as shown in figure 13), which causes a pressure rise that propagates upstream to the MS and subsequently results in the incompatibility of the pressure ratio across the MS and its relative propagation speed. Then, the upstream motion of the MS enters a positive feedback period, and the detonation front moves out of the combustor. At other flight Mach numbers (i.e. *Ma* = 9.5 and 10.5), destabilized Mach reflection patterns of ODWs also occur at large *ζ* values, as shown in figure 9, and similar analyses suggest that the relevant destabilization mechanisms are also similar. Since this kind of destabilization phenomenon is induced by the action of a shock wave (the tMS in this section), this process is referred to as wave-induced destabilization in this paper.

### 3.4. Inherently destabilized Mach reflection of ODWs at low *Ma*

As depicted in figure 9, for a relatively large Mach number (*Ma* = 9.5, 10 or 10.5), when the *ζ* (>0) value is small (i.e. the designed ODW reflection point is located close to the expansion corner), the flow structures of the Mach reflection of an ODW remain stabilized in the combustor, whereas the detonation front moves out of the combustor upstream for large *ζ* values (i.e. the designed ODW reflection point is located far upstream from the expansion corner), implying the destabilization of the Mach reflection pattern of the ODW. However, a special phenomenon of stabilization/destabilization occurs for a relatively low Mach number, for example *Ma* = 9, as shown in figure 14. For *Ma* = 9 and *ζ* = 12/60, the Mach reflection of an ODW is destabilized in the combustor, and the evolution of the MS is similar to the evolutionary characteristics in the destabilized cases of *Ma* = 10. After the Mach reflection pattern of an ODW forms (corresponding to the first peak on the speed curve in figure 14*a* for this case), the upstream motion speed of the MS initially drops rapidly and then decreases gradually. During the stabilization of the MS, the motion speed of the MS suddenly rises, and the MS accelerates upstream. As a result, the MS moves out of the combustor to a location far upstream from the expansion corner. According to the previous analyses, one may expect to stabilize the Mach reflection of the ODW in the combustor by decreasing the *ζ* value, i.e. by setting the designed ODW reflection point closer to the expansion corner. However, while this guideline prevails for *Ma* = 9.5, 10 and 10.5, this solution fails for *Ma* = 9. By continuously decreasing the *ζ* value from 12/60 to 6/60, 3/60 or even 1/60, the motion speed of the MS still exhibits a sudden increase, eventually leading to the destabilization of the flow field. One of the key differences among the different *ζ* values is the time instant at which the motion speed of the MS rises; in other words, the time instant of the turning point occurs later for smaller *ζ* values. Notably, the designed reflection point is only 1 mm upstream of the expansion corner for a 6 cm high combustor with *ζ* = 1/60, but the destabilization of the Mach reflection of an ODW still occurs for this extremely small value of *ζ*. Remarkably, the flow structures at this low *Ma* are stable with *ζ* = 0, where the ODW reflects precisely at the expansion corner and the Mach reflection pattern of an ODW does not form. The above analyses suggest that destabilization is inherent in the Mach reflection of an ODW at this low *Ma*, when the ODW reflects off the wall before the expansion corner.

To determine the inherent mechanism causing the sudden increase in the motion speed of the MS in these low-*Ma* cases that directly destabilizes the Mach reflection of ODWs, the instantaneous flow fields and pressure distributions along the wall before the expansion corner are analysed near the turning time instant, as shown in figures 15 and 16, respectively, taking *Ma* = 9 and *ζ* = 12/60 as an example. Figure 15(*a*) demonstrates that the RfSW2 crosses the SL directly and then reflects off the lower wall at this low *Ma* rather than reflecting over the SL and forming a tMS at a relatively high *Ma* (e.g. figure 11*a*). Although slightly different flow structures are formed at this low *Ma*, the RfSW2 acts on the subsonic zone behind the MS at approximately *t* = 0.446 ms (figure 15*b*) as the MS moves upstream continuously. Thereafter, the pressure rise caused by the RfSW2 propagates upstream through this subsonic zone towards the detonation front of the MS, as shown in figures 15(*c*–*e*) and 16. After this pressure wave arrives at the detonation front of the MS at approximately *t* = 0.560 ms, the pressure behind the MS and the pressure ratio across the MS both increase; the resulting incompatibility between the pressure ratio across the MS and its relative propagation speed leads to a sudden increase in the motion speed of the MS (point ‘e’ in figure 14*a*), ultimately causing the destabilization of the detonation front in the combustor (figures 15*f* and 15*g*).

It is suggested that the Mach reflection of an ODW before an expansion corner (*ζ* > 0) can be directly destabilized by the action of the RfSW2 or its tMS on the subsonic zone behind the MS regardless of the flight Mach number. Indeed, this is easy to understand with the above analyses. As the MS moves upstream, the RfSW2 and/or its tMS moves upstream and approaches the MS, and when the MS moves sufficiently far upstream from the expansion corner, the RfSW2 and/or its tMS acts on the subsonic zone and eventually causes destabilization. For a relatively high *Ma* (e.g. *Ma* = 10), the stabilized distance of the MS from the expansion corner increases as *ζ* increases. Consequently, a small *ζ* value results in a short distance of the MS from the expansion corner, which prevents the tMS of the RfSW2 from acting on the subsonic zone. On the other hand, a potentially long distance of the MS from the expansion corner is implied with a large *ζ* value, leading to destabilization. However, destabilization always occurs at *Ma* = 9, even with an extremely small *ζ* value of 1/60, implying that a stabilized MS location may not exist for this low Mach number and consequently that the MS may continuously move upstream once the Mach reflection pattern of the ODW forms. In other words, destabilization may be inherent in the Mach reflection of an ODW at a low *Ma*, and the action of the RfSW2 or its tMS on the subsonic zone behind the MS may only accelerate the upstream motion of the MS. Hence, this kind of destabilization phenomenon at a low *Ma* is referred to as inherent destabilization in this paper, and an in-depth theoretical analysis is further conducted to understand the relevant mechanisms.

### 3.5. Formulating the stabilized location of the MS

In the following formulation, the effects of the RfSW2 and its tMS are neglected since they are found to only accelerate the inherent upstream motion of the MS at a low *Ma*, as discussed in § 3.4. According to the above analyses of the instantaneous flow fields, the SL is inclined downward relative to the incoming flow, as is schematically shown in figure 17. Hence, the subsonic flow crossing the MS accelerates towards becoming sonic in the contracted flow path bounded by the SL and the lower wall, and the geometric throat is positioned at the expansion corner (point O). Assume that the MS is a straight line perpendicular to the lower wall and that the SL is also a straight line with an inclination angle of *α* relative to the incoming flow. According to the geometric relationship exhibited in figure 17, the length of the MS can be calculated by

where *β* is the oblique detonation angle of the ODW relative to the incoming flow. The distance from the expansion corner to the SL, namely $\overline {\textrm{OL}} $, which is considered as the width of the geometric throat, can be expressed as

where $\overline {\textrm{ON}} $ is a segment parallel to MS. Therefore, the contraction ratio of the flow path from the MS to the geometric throat $\overline {\textrm{OL}} $ appears as

Equation (3.3) reveals that as the S point (the triple point) of the MS moves upstream along the ODW (i.e. the triple point is always located on the ODW), the ratio of $\overline {\textrm{RO}} \textrm{/}\overline {\textrm{MO}} $ (<1) decreases, and the contraction ratio of the flow path, $\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} $, decreases accordingly. If the MS is located close to the expansion corner, the contraction ratio of the flow path behind the MS is large enough that a sonic state may be achieved before the geometric throat $\overline {\textrm{OL}} $, as schematically shown in figure 17(*a*) (the location of the carmine dashed line corresponds to the contraction ratio of the flow path of the sonic state). Notably, the streamlines behind the MS are generally curved to adjust the flow direction to conform to the flow path between the SL and the lower wall. Further contraction of the flow path behind the sonic line would directly choke the flow. As a result, the MS would move upstream to adjust the flow behind the MS to ensure that the sonic state occurs at the geometric throat $\overline {\textrm{OL}} $. Obviously, the geometric converging effects of the flow path before the expansion corner contribute to sustaining the MS. In contrast, if the MS is located far upstream from the expansion corner (figure 17*b*), the contraction ratio becomes so small that the local Mach number at the geometric throat $\overline {\textrm{OL}} $ is still smaller than one. In other words, under this circumstance, the flow behind the MS may be subsonic throughout the flow path before and across the expansion corner. Then, the geometric diverging effects of the flow path behind the expansion corner could propagate upstream through left-running characteristic lines across the expansion corner and arrive at the MS, which weakens the geometric converging effects of the flow path before the expansion corner. As a result, the MS would be blown away by the upstream supersonic flow. It appears that the stabilized location of the MS (i.e. $\overline {\textrm{MO}} $) corresponds to a critical condition in which the local Mach number at the geometric throat $\overline {\textrm{OL}} $ is equal to one, which is similar to the Kantrowitz limit in the unstart problems of scramjet inlets (Curran & Murthy Reference Curran and Murthy2001). In such a problem, the Kantrowitz-limit contraction ratio is determined by the one-dimensional isentropic area ratio that accelerates a subsonic flow (behind a normal shock wave) to a sonic flow at the throat. In the present problem regarding the stabilization of the Mach reflection of an ODW, the similar Kantrowitz-limit contraction ratio, *CR _{Kantrowitz}*, can be evaluated for a specific NDW, and it can be employed to determine the theoretical stabilization location of the MS.

From the above analysis, the critical condition for the stabilization of the MS appears as

or

Substituting (3.3) into (3.5) yields the formula for calculating the stabilized location of the MS, expressed as

Obviously, with the parameters *α*, *β* and *CR _{Kantrowitz}* being known, the stabilized location of MS (i.e. $\overline {\textrm{MO}} $) corresponding to the designed reflection location of ODW (i.e. $\overline {\textrm{RO}} $) can be evaluated directly from (3.6).

To evaluate the oblique detonation angle of the ODW, *β*, the inclination angle of the SL, *α*, and the Kantrowitz-limit contraction ratio behind the MS, *CR _{Kantrowitz}*, for different cases, the chemical equilibrium hypothesis is employed in the present study to account for changes in the chemical compositions and gas thermodynamic properties and to more accurately predict the heats of chemical reactions; that is, the chemical reaction rates are assumed to be infinitely fast, and chemical equilibrium states are assumed to be achieved everywhere in the flow field behind the leading shock/detonation fronts (i.e. behind the ODW and MS). Further, a two-step cyclic iterative solution method proposed by the authors (Zhang Reference Zhang2020

*a*; Zhang

*et al.*Reference Zhang, Wen, Zhang, Liu and Jiang2021

*a*) is adopted to solve the chemical equilibrium solutions of different processes with high accuracy, including the chemical equilibrium normal shock relationship, the chemical equilibrium strong and weak oblique shock relationships (where ‘strong’ and ‘weak’ refer to the strong and weak branches of the oblique shock solutions, respectively) and the chemical equilibrium isentropic relationship. For clarity, the solution processes of the above chemical equilibrium relationships are briefly presented in Appendices A–E; for the detailed derivations and evaluations of the solution accuracy, the reader may refer to Zhang (Reference Zhang2020

*a*) and Zhang

*et al.*(Reference Zhang, Wen, Zhang, Liu and Jiang2021

*a*). Specifically, for different cases, the oblique detonation angle of the ODW,

*β*, can be directly evaluated by the chemical equilibrium weak oblique shock relationship; the predicted values of

*β*in different cases agree well with those observed in the numerical simulations in the present study. For the evaluation of the Kantrowitz-limit contraction ratio behind the MS,

*CR*, the flow state behind the MS is first obtained through the chemical equilibrium normal shock relationship; then, to evaluate the MS–throat area ratio through the chemical equilibrium isentropic relationship, the subsonic flow state needs to accelerate to the sonic state in the contracted flow path between the SL and the lower wall, during which the chemical compositions change due to shifts in chemical equilibrium. Moreover, to determine the inclination angle of the SL,

_{Kantrowitz}*α*, the three-shock theory (Ben-Dor Reference Ben-Dor2007) is applied to the neighbourhood of the triple point, as schematically shown in figure 18. The MS is slightly inclined with respect to the perpendicular direction of the incoming flow in the neighbourhood of the triple point S, and hence the flow state behind this small segment of the MS needs to be evaluated through the chemical equilibrium strong oblique shock relationship, while the change in the flow state across the RfSW (from region 2 to region 3, also based on the assumption of chemical equilibrium) is calculated by the chemical equilibrium weak oblique shock relationship. Thereafter, to determine the inclination angle of the SL,

*α*, the pressure balance between both sides of the SL, namely between region 3 behind the RfSW that deflects the flow by an angle of $\theta - \alpha$ (

*θ*= 25° in this paper) and region 4 behind the inclined MS that deflects the flow by an angle of

*α*, i.e.

*p*

_{3}=

*p*

_{4}(figure 18), needs to be satisfied. With the chemical equilibrium strong and weak shock relationships, this pressure-balanced state can be solved by utilizing a simple dichotomizing search algorithm.

From the above analyses, the parameters *α*, *β* and *CR _{Kantrowitz}* can be uniquely determined (based on the assumption of chemical equilibrium) according to the inflow conditions before the combustor, which depend uniquely on

*Ma*, as shown in figure 19. The specific values of these parameters for the stabilized flight Mach numbers (i.e.

*Ma*= 9.5, 10 and 10.5) are summarized in table 2. As

*Ma*decreases, both the chemical equilibrium inclination angle of the SL,

*α*, and the chemical equilibrium oblique detonation angle of the ODW,

*β*, increase, whereas the chemical equilibrium Kantrowitz-limit contraction ratio behind the MS,

*CR*, decreases. Because the flow Mach number behind the chemical equilibrium MS is closer to one at a lower

_{Kantrowitz}*Ma*, it is easier for this subsonic flow to accelerate to the sonic state, leading to a smaller limiting contraction ratio. Substituting these parameters into (3.6), the predicted stabilized locations of the MS (i.e. $\overline {\textrm{MO}} $) at

*Ma*= 9.5, 10 and 10.5 as a function of the designed ODW reflection location (i.e. $\overline {\textrm{RO}} $) are plotted as the dashed lines in figure 20. For comparison, the real values of $\overline {\textrm{MO}} $ observed in the numerical simulations are also presented in figure 20, revealing that (3.6) underestimates the stabilized locations of the MS to a large extent compared with the numerical simulations.

These discrepancies in the stabilized locations of the MS between the theoretical predictions by (3.6) and the real observations via numerical simulations result from the following factors. First, the real SL between the triple point and the geometric throat is not strictly straight, as shown in figure 21(*a*) (taking the case of *Ma* = 10 and *ζ* = 5/60 as an example). Due to the limited space in the neighbourhood of the triple point compared with the chemical reaction characteristic length, the chemical equilibrium state is not fully achieved near the triple point, resulting in the local inclination angle of the SL being larger than the corresponding chemical equilibrium value, as depicted in figure 21(*b*). As the fluid flows downstream, the flow state varies to the chemical equilibrium state, and hence the local inclination angle of the SL decreases gradually to the corresponding chemical equilibrium value. However, the segment of the SL close to the geometric throat $\overline {\textrm{OL}} $ becomes increasingly inclined downward under the effects of the expansion fan emitted from the expansion corner (figure 21*a*), leading to an increase in the local inclination angle of the SL and another deviation from its chemical equilibrium value (figure 21*b*). Influenced by these two effects, the average inclination angle of the SL from the triple point to the geometric throat is undoubtedly larger than that predicted by the chemical equilibrium theory. Second, due to the local chemical non-equilibrium effects of the flow behind the MS, the change in the flow from the MS to the throat does not exactly follow the chemical equilibrium isentropic process. Furthermore, due to the special one-sided expansion geometry of the expansion corner and the finite-length contracted flow path in the present problem, the one-dimensional assumption in the traditional Kantrowitz limit is no longer satisfied (Curran & Murthy Reference Curran and Murthy2001), and hence the flow at the geometric throat $\overline {\textrm{OL}} $ is not exactly at the sonic state; in fact, the Mach number of this flow is smaller than one and varies along $\overline {\textrm{OL}} $, as shown, for example, in figures 4(*g*) and 8. As a result, the real contraction ratio from the MS to the geometric throat $\overline {\textrm{OL}} $, i.e. $\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} $, deviates slightly from the chemical equilibrium Kantrowitz-limit contraction ratio, *CR _{Kantrowitz}*.

To consider the above factors causing the discrepancies in the theoretical predictions of the stabilized locations of the MS compared to the numerical observations, a corrected parameter, *κ*, is introduced to the present theoretical formulation, and the critical condition for the stabilization of the MS, i.e. (3.5), is rewritten as

Then, a corrected formula for calculating the stabilized location of the MS is obtained by substituting (3.3) into (3.7), expressed as

To evaluate the corrected parameter, *κ*, (3.8) is rearranged into

which is substituted with the numerical values of $\overline {\textrm{MO}} $ in the stabilized cases at *Ma* = 9.5, 10 and 10.5 summarized in table 2. The value of *κ* ranges from approximately 1.19 to 1.33 for these stabilized Mach reflection cases. Although the corrected parameter, *κ*, does not change much between different cases, it appears to be a function of the flight Mach number, *Ma*, and the dimensionless reflection location, *ζ*, i.e. *κ* = *κ*(*Ma*, *ζ*), as depicted in figure 22.

One purpose of this section is to provide theoretical predictions of the stabilized locations of the MS at *Ma* = 9 and small *ζ* values, especially its asymptotic behaviour as *ζ* goes to zero at this relatively low *Ma*. Since the values of *Ma* and *ζ* do not change much compared to those in the stabilized cases summarized in table 2, a second-order Taylor expansion approximation of the function of $f(Ma,\zeta ) \equiv {(\kappa \text{--}1)^{ - 1}}$ can be applied to fit the corrected parameter, *κ*, as follows:

After rearranging the above equation, it can be rewritten as

where *a* _{1} = 187.3, *a* _{2} = −35.23, *a* _{3} = 67.63, *a* _{4} = 1.680, *a* _{5} = −5.364 and *a* _{6} = −5.937, as determined by the least squares method using the data of *κ* given in table 2. The fitting curves of *κ* = *κ*(*Ma*, *ζ*) by (3.11) are presented in figure 22. The fitting curves and the numerical data of the stabilized cases exhibit good agreement. Further, the corrected theoretical predictions of the stabilized locations of the MS (i.e. $\overline {\textrm{MO}} $) using (3.8) and (3.11) are also presented in figure 20 as the solid lines, demonstrating good consistency between the corrected theoretical predictions and the numerical observations of the stabilized cases.

For *Ma* = 9, the chemical equilibrium oblique detonation angle of the ODW and the inclination angle of the SL are *β* = 55.34° and *α* = 13.16°, respectively, while the corresponding chemical equilibrium Kantrowitz-limit contraction ratio of the subsonic flow behind the MS is *CR _{Kantrowitz}* = 1.3790. For

*ζ*= 1/60 (i.e. $\overline {\textrm{RO}} = 1\;\textrm{mm}$), the predicted value of the corrected parameter using (3.11) is

*κ*= 1.151. Substituting all the above parameters into (3.8), the theoretical prediction of the ratio of $\overline {\textrm{MO}} \textrm{/}\overline {\textrm{RO}} $ is equal to −7.454, which is smaller than zero. This result obviously contradicts the geometric relationship depicted in figure 17, i.e. $\overline {\textrm{MO}} \textrm{/}\overline {\textrm{RO}} > 1 > 0$, implying that a stabilized location of the MS (i.e. $\overline {\textrm{MO}} $) cannot be achieved in this low-

*Ma*case. Figure 23 shows the changes in the value of

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ as the MS moves upstream at

*ζ*= 1/60 (i.e. $\overline {\textrm{RO}} = 1\;\textrm{mm}$) for different

*Ma*, revealing that the value of

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ is greater than zero when $\overline {\textrm{MO}} $ is small for all

*Ma*. According to the previous analysis of figure 17 but with the corrected critical condition of (3.7), when

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}} > \textrm{ }0$, the flow becomes choked in the contracted flow path bounded by the SL and the lower wall, resulting in the MS moving upstream. In the cases of

*Ma*= 9.5, 10 and 10.5,

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ approaches its zero point as $\overline {\textrm{MO}} $ increases, and, finally, the MS remains stabilized at the corresponding location. In contrast, if the change in

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ passes through its zero point to a negative value in the cases of low

*Ma*(

*Ma*= 9.5, 10 and 10.5), the MS becomes attenuated downstream under the geometric diverging effects of the expansion corner; as a result,

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ returns to its zero point, and the MS remains stabilized. Therefore, the existence of zero points of

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ for these

*Ma*values implies the potential stabilization of the MS at a specific location. However, the value of

*κ*$\overline {\textrm{MS}} \textrm{/}\overline {\textrm{OL}} \text{--}C{R_{Kantrowitz}}$ is always greater than zero, and there is no zero point in the evolution of $\overline {\textrm{MO}} $ in the case of

*Ma*= 9. Therefore, the MS grows and moves upstream continuously after it forms, leading to the destabilization of the Mach reflection of the ODW in this low-

*Ma*case. Further, as depicted in figure 20, the predicted stabilized location of the MS at

*Ma*= 9 is always negative for

*ζ*> 0, implying that the Mach reflection of ODWs is inherently destabilized at this low

*Ma*due to flow choking, when the ODW reflects before an expansion corner, i.e.

*ζ*> 0, as discussed in § 3.4.

According to the previous analyses, the stabilization of the Mach reflection of an ODW before an expansion corner is feasible at small *ζ* because destabilization would be directly triggered by the action of the RfSW2 or its tMS on the subsonic zone behind the MS at large *ζ*. Due to the limited spatial resolution of numerical simulations, it is impossible to infinitely decrease *ζ* to check whether the Mach reflection of ODWs before an expansion corner at a specific *Ma* is inherently destabilized. However, the inherent stabilization or destabilization characteristics at a specific *Ma* can be confirmed by the asymptotic behaviour of $\overline {\textrm{MO}} \textrm{/}\overline {\textrm{RO}} $ when $\zeta \to 0 +$, as shown in figure 24. The limit of $\overline {\textrm{MO}} \textrm{/}\overline {\textrm{RO}} $ is greater than zero when *Ma* is greater than approximately 9.3, implying that when *ζ* is small enough that the RfSW2 or its tMS does not act on the subsonic zone behind the MS, the Mach reflection of an ODW before the expansion corner can be stabilized. For *Ma* smaller than approximately 9.3, the limit of $\overline {\textrm{MO}} \textrm{/}\overline {\textrm{RO}} $ is smaller than zero, implying that the Mach reflection of an ODW before an expansion corner is inherently destabilized regardless of how small *ζ* (>0) is. In other words, the critical Mach number of inherently stabilized/destabilized Mach reflection of ODWs before an expansion corner, *Ma _{cri}*, predicted by the present theory, is approximately 9.3. The present theoretical prediction is consistent with the numerical observations summarized in figure 9, in which the Mach reflections of ODWs can be stabilized with a small

*ζ*at

*Ma*≥ 9.5 but exhibit destabilization with all

*ζ*> 0 at

*Ma*= 9. The numerical prediction of

*Ma*appears between 9 and 9.5.

_{cri} Further, figure 19 shows that as *Ma* decreases, the inclination angle of the SL, *α*, increases, resulting in an increase in the real contraction extent of the flow path bounded by the SL and the lower wall, whereas the corresponding Kantrowitz-limit contraction ratio of the subsonic flow behind the MS decreases. Hence, the subsonic flow behind the MS is easily choked before the expansion corner at a low *Ma*. In other words, this subsonic flow may accelerate to the sonic state before the expansion corner and be choked by further area contraction, resulting in the MS continuously moving upstream and inherent destabilization, as shown in figure 17(*a*). In summary, the theoretical results presented in this section demonstrate that the inherent destabilization of the Mach reflection of an ODW before an expansion corner at a low *Ma* is attributed to the geometric flow choking phenomenon caused by the large flow path contraction behind the MS.

## 4. Conclusions

In this paper, the inviscid Mach reflection patterns of ODWs before an expansion corner in a space-confined ODE combustor are numerically studied by solving the two-dimensional time-dependent multispecies reactive Euler equations in combination with a detailed hydrogen combustion mechanism. The effects of different flight Mach numbers (*Ma*) and different dimensionless reflection locations (*ζ* ≥ 0) on the stabilization/destabilization characteristics of the flow structures in the combustor and the relevant mechanisms are emphasized. The results suggest that a stabilized regular reflection pattern occurs at all *Ma* values when the ODW reflects precisely at the expansion corner, i.e. *ζ* = 0. However, if the ODW reflects upstream of the expansion corner, i.e. *ζ* > 0, a Mach reflection pattern emerges, and complex flow structures, such as the MS, SL, RfSW, RfSW2 and subsonic zone behind the MS, form in the combustor. The stabilization/destabilization characteristics of these formed flow structures in the Mach reflection pattern of the ODW depend on the values of *Ma* and *ζ*.

When *Ma* is not very low, e.g. *Ma* = 9.5, 10 and 10.5, and the *ζ* value is small (>0), the flow structures of the Mach reflection of the ODW remain stabilized in the combustor. Under this circumstance, the stabilized location of the MS and the downstream flow structures move upstream as *ζ* increases. In particular, the distance between the RfSW2 or its tMS and the subsonic zone behind the MS closes with increasing *ζ*. Further, when *ζ* increases beyond a specific value, wave-induced destabilization emerges. The RfSW2 or its tMS acts on the subsonic zone behind the MS, and the relevant pressure rise propagates through this subsonic zone towards the detonation front of the MS, leading to sudden increases in the pressure behind the MS and the pressure ratio across the MS. Although the upstream motion of the MS is seemingly stabilized initially, the incompatibility of the pressure ratio across the MS and its propagation speed relative to the incoming flow causes the MS to suddenly accelerate upstream. Then, the upstream motion of the MS enters a positive feedback stage, and its upstream motion speed continuously increases. Finally, the detonation front of the ODW vanishes, and the MS moves out of the combustor, resulting in the destabilization of the Mach reflection of the ODW at this large *ζ* value.

For a low *Ma*, for example *Ma* = 9, the Mach reflection patterns of ODWs are destabilized at any *ζ* > 0 (even as small as *ζ* = 1/60), implying inherent destabilization, although the reflection flow field of an ODW is stabilized at *ζ* = 0 at this low *Ma* or at small *ζ* > 0 at a higher *Ma*. Analyses of the evolution of the flow fields and pressure distributions reveal that the inherent destabilization of the Mach reflection of an ODW at a low *Ma* is also directly induced by the action of the RfSW2 or its tMS on the subsonic zone behind the MS during the inherently continuous upstream movement of the MS. Further analyses conducted by theoretically formulating the stabilized location of the MS suggest that its inherently continuous upstream motion and consequently the inherent destabilization of the flow field are attributed to choking of the flow in the contracted flow path between the SL and the lower wall at low *Ma* values, where the inclination angle of the SL relative to the incoming flow is large and hence the geometric contraction of the flow path is large.

The above results imply that inviscid destabilization of the flow field occurs easily when the ODW reflects off the wall before the expansion corner of the combustor, causing a Mach reflection pattern to form, which is unfavourable to the operation of an ODE combustor because failure of the engine may occur. Therefore, this study clearly demonstrates that, to avoid the potential destabilization of the Mach reflection patterns of ODWs, the reflection point of the ODW should not be designed to be located before the expansion corner. Instead, the reflection point of the ODW should be designed downstream of the expansion corner, as suggested by the authors’ previous research (Zhang *et al.* Reference Zhang, Ma, Zhang, Han, Liu and Jiang2020*b*, Reference Zhang, Wen, Zhang, Liu and Jiang2021*b*), where the stabilized viscous flow fields of ODWs or ODW reflections were obtained, even for a low *Ma* (i.e. *Ma* = 9). Moreover, future studies can focus on (1) the effects of viscous factors, especially boundary layers and shock/detonation-induced boundary layer separations, on the stabilization characteristics of an ODW when the ODW is designed to reflect downstream of the expansion corner, (2) the relevant flow mechanisms and (3) the stabilization-control strategies in the engine combustor.

## Funding

This work is supported by the National Natural Science Foundation of China (grant nos 11772284, 11672312 and 11532014).

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Two-step cyclic iterative solution method for different chemical equilibrium relationships

Problem description: the pressure *p* _{1}, temperature *T* _{1}, velocity *u* _{1} and chemical compositions *χ _{i}*

_{1}(in molar fraction,

*i*= 1, …,

*nSp*, where

*nSp*is the number of species considered in the mixture) are known before a reactive normal shock wave or OSW or before starting an isentropic process; additionally, the deflection angle

*θ*is known for the reactive oblique shock problem, or the target pressure

*p*

_{2}is known for the isentropic process. By assuming infinite chemical reaction rates (i.e. chemical equilibrium), the pressure

*p*

_{2}(except for the isentropic process), temperature

*T*

_{2}, velocity

*u*

_{2}, chemical compositions

*χ*

_{i}_{2}and oblique shock angle

*β*(only for an oblique shock) behind the reactive normal shock wave or OSW or after the isentropic process need to be solved.

The method adopted to solve the above chemical equilibrium relationships, namely the chemical equilibrium normal shock relationship, the chemical equilibrium strong or weak normal shock relationship and the chemical equilibrium isentropic relationship, is the two-step cyclic iterative solution method, which was previously proposed by the authors (Zhang Reference Zhang2020*a*; Zhang *et al.* Reference Zhang, Wen, Zhang, Liu and Jiang2021*a*). A flow diagram illustrating the solution process is presented in figure 25, where *λ* = 0.4 is a relaxation factor to ensure stability and convergence in the iteration process and *k* is an index of the iteration step. Moreover, the convergence criteria in figure 25 can be chosen as follows:

For the calculation of the chemical equilibrium normal shock relationship, Step 1 in figure 25 is specified as the calculation of the multispecies normal shock relationship by giving *p* _{1}, *T* _{1}, *u* _{1}, *χ _{i}*

_{1}and

*χ*

_{i}_{2}, and the detailed calculation process is given in Appendix B.

For the calculation of the chemical equilibrium strong or weak normal shock relationship, Step 1 in figure 25 is specified as the calculation of the multispecies strong or weak oblique shock relationship by giving *p* _{1}, *T* _{1}, *u* _{1}, *θ*, *χ _{i}*

_{1}and

*χ*

_{i}_{2}, and the detailed calculation process is given in Appendix C.

For the calculation of the chemical equilibrium isentropic relationship, Step 1 in figure 25 is specified as the calculation of the multispecies isentropic relationship by giving *p* _{1}, *T* _{1}, *u* _{1}, *χ _{i}*

_{1},

*χ*

_{i}_{2}and

*p*

_{2}, and the detailed calculation process is given in Appendix D.

For the calculations of all the above chemical equilibrium relationships, the determination of the chemical equilibrium compositions, *χ _{i}*

_{2}, at the given

*p*

_{2}and

*T*

_{2}, i.e. Step 2 in figure 25, is given in Appendix E in detail.

## Appendix B. Calculation of the normal shock relationship with varying compositions

Problem description: the pressure *p* _{1}, temperature *T* _{1}, velocity *u* _{1} and chemical compositions *χ _{i}*

_{1}before a normal shock wave are known. Additionally, the chemical compositions

*χ*

_{i}_{2}behind the normal shock are also known, but they do not need to be the same as

*χ*

_{i}_{1}to account for chemical reactions. Then, the post-shock pressure

*p*

_{2}, temperature

*T*

_{2}and velocity

*u*

_{2}are solved. A single-variable Newton iteration method derived from the governing equations of a normal shock wave is adopted here to solve the post-shock velocity

*u*

_{2}(the iteration variable). The iteration function

*f*(

_{NSW}*u*

_{2}) is expressed as

where *R* _{1} and *R* _{2} are the gas constants of the mixtures before and behind the normal shock wave, respectively, and they can be easily calculated from the relevant chemical compositions of mixtures. Moreover, *h* _{1} and *h* _{2} are the specific enthalpies of the mixtures before and behind the normal shock wave, respectively, and they can be evaluated at the given temperature from a species thermodynamic database with known mixture compositions, for example, the piecewise fourth-order temperature polynomial fitting of NASA (McBride, Gordon & Reno Reference Mcbride, Gordon and Reno1993). The Newton iteration formula appears as

where the initial iteration value can be set as *u* _{2,0} = 0. In (B2), the derivative of the iteration function, ${f^{\prime}_{NSW}}({u_2})$, is expressed as

where *c _{p}*

_{2}is the specific heat at constant pressure of the mixture behind the normal shock wave. After calculating

*u*

_{2}, other unknown post-shock parameters can be obtained directly by substituting

*u*

_{2}into the normal shock governing equations.

## Appendix C. Calculation of the strong or weak oblique shock relationship with varying compositions

Problem description: the pressure *p* _{1}, temperature *T* _{1}, velocity *u* _{1} and chemical compositions *χ _{i}*

_{1}before an OSW and the flow deflection angle

*θ*across the oblique shock are known. Additionally, the chemical compositions

*χ*

_{i}_{2}behind the oblique shock are also known, but they do not need to be the same as

*χ*

_{i}_{1}to account for chemical reactions. Then, the post-shock pressure

*p*

_{2}, temperature

*T*

_{2}and velocity

*u*

_{2}and the oblique shock angle

*β*are solved. A similar single-variable Newton iteration method derived from the governing equations of an oblique shock wave is also adopted. Here, the iteration variable is the oblique shock angle

*β*, and the iteration function

*f*(

_{OSW}*β*) appears as

Then, the Newton iteration formula is

with the derivative of the iteration function, ${f^{\prime}_{OSW}}(\beta )$, expressed as

In general, there are two branches of oblique shock solutions, namely the strong and the weak oblique shock solutions. To ensure that the iteration process converges to the different solution branches of OSWs exactly, different iteration initial values are used as follows:

where *Ma* _{1} is the Mach number before the oblique shock. After *β* is known, other unknown post-shock parameters can be obtained directly by substituting *β* into the oblique shock governing equations.

## Appendix D. Calculation of the isentropic relationship with varying compositions

Problem description: the pressure *p* _{1}, temperature *T* _{1}, velocity *u* _{1} and chemical compositions *x _{i}*

_{1}(moles of species per unit mass of the mixture) before starting an isentropic process are known. Notably, the mixture must be in chemical equilibrium before starting an isentropic process with varying compositions. Additionally, the chemical compositions

*x*

_{i}_{2}(they also do not need to be the same as

*x*

_{i}_{1}to account for shifts in chemical equilibrium during the isentropic process) and the pressure

*p*

_{2}after the isentropic process are known. Then, the temperature

*T*

_{2}and velocity

*u*

_{2}after the isentropic process are solved. A single-variable Newton iteration method is employed again, which is directly derived from the isentropic condition of the mixture. The iteration variable is the temperature

*T*

_{2}, and the iteration function

*f*(

_{isen}*T*

_{2}) is expressed as

where $S_i^0(T)$ is the species standard molar entropy at temperature *T* and can also be obtained from a species thermodynamic database; *R* _{u} = 8.314 J mol^{−1} K^{−1} is the universal gas constant; and *p* _{atm} = 1.01325 × 10^{5} Pa is the pressure value defined under standard conditions. The Newton iteration formula appears as

with the derivative of the iteration function, ${f^{\prime}_{isen}}({T_2})$, given by

In (D3), $C_{pi}^0(T)$ is the species constant-pressure molar heat capacity at temperature *T*. After *T* _{2} is known, *u* _{2} can be obtained directly through the energy conservation equation. Furthermore, if an area ratio of the isentropic process is needed, it can be obtained easily by substituting *p* _{2}, *T* _{2} and *u* _{2} into the equation of state of the mixture and the mass conservation equation.

## Appendix E. Calculation of chemical equilibrium compositions at a specific temperature and pressure

Problem description: the initial compositions of the mixture are known, i.e. the total moles of atoms of element *j*, *b _{j}* (

*j*= 1, …,

*nEl*, where

*nEl*is the number of elements involved in the mixture), and the number of atoms of element

*j*in species

*i*,

*a*, are known. Then, the equilibrium compositions at a specific temperature

_{ij}*T*and pressure

*p*are solved. The minimization of free energy method of NASA (Gordon & Mcbride Reference Gordon and Mcbride1994) is used. Let

*y*(moles of species per unit mass of the mixture) be an estimate of the equilibrium compositions; then, the new estimate

_{i}*x*can be evaluated by

_{i}
where *d _{i}* is defined as

In (E2), $G_i^0$ is the species standard molar Gibbs free energy. Moreover, the unknowns *∑x _{i}*/

*∑y*and

_{i}*λ*in (E1) are the solutions of the following linear equations:

_{j}
After *x _{i}* has been solved, it serves as the new

*y*in an iterative process until convergence is achieved.

_{i}