2.1. Definition of asymmetric shock reflection in an overexpanded jet

Figure 2 is a schematic illustration of asymmetric Mach reflection in an overexpanded jet. The shaded triangles represent nozzle exit geometry, where the exit velocity is aligned with the x-axis. The nozzle height $A_{1}A_{1}'$ is $H_{A}$ , nozzle exit Mach number is $M_{0}$ and nozzle pressure ratio (NPR) is denoted as the ratio of the stagnation chamber pressure to the ambient pressure. The asymmetry originates from the offset in x-coordinates between upper and lower nozzles. To quantitatively analyse this asymmetry, we define a dimensionless form of the extension length of the lower nozzle ( $A'_{1}A_{2}$ in figure 2) as

(2.1) \begin{align} L_w=\frac {L_{A'_{1}A_{2}}}{H_A}. \end{align}

When $L_w \neq 0$ , expansion waves reflected at the lower jet boundary interact with the slipline first, causing flow deflection downward and resulting in an asymmetric overexpanded jet.

Figure 2.

Schematic illustration of asymmetric Mach reflection in an overexpanded jet.

Moreover, it is evident that $L_w$ cannot increase indefinitely. As shown in figure 3(a), when the extension line of the incident shock wave at the upper lip position $A_{1}$ intersects precisely with $A'_{1}A_{2}$ , the asymmetric shock reflection does not exist. Denoting the corresponding $L_{w}$ as $L_{w, \textit{max}}$ ,

(2.2) \begin{align} L_{w, \textit{max}}=\frac {1}{\tan \beta _{i1}} ,\end{align}

where $\beta _{i1}$ is the shock angle of the upper incident shock wave, which can be calculated using the classical oblique shock relations.

Figure 3.

Flow structure and its predictive parameter under $L_{w, \textit{max}}$ . (a) Schematic illustration of flow structure corresponding to $L_{w, \textit{max}}$ . (b) Variation of the predictive $L_{w, \textit{max}}$ at different NPR when given $M_{0}$ .

Given NPR, ambient pressure and $M_{0}$ , $p_{0}$ can be obtained via the stagnation pressure–static pressure relation for isentropic flow. Then, $\beta _{i1}$ can be obtained using oblique shock relations. From (2.2), the variation of $L_{w, \textit{max}}$ with NPR under different Mach numbers is derived, as shown in figure 3(b). Here, $L_{w, \textit{max}}$ increases with the increase of NPR and decreases with the increase of Mach number.

2.2. Numerical simulation method

The boundary conditions are illustrated in figure 4, where AH, BC are set as the wall, AB is set as the pressure inlet, EF is set as the pressure outlet, and CD, DE, FG, GH are set as the pressure far field. To facilitate subsequent verification of the computational domain independence, the dimensions of the computational domain are defined as $L_{1}$ and $L_{2}$ .

Figure 4.

Boundary condition for CFD calculation.

We solve the full set of nonlinear Euler equations in gas dynamics using the well-known Roe scheme based on the finite difference approximation and second-order upwinding for the flux. The solution is considered converged when the size and shape of the Mach stem remain unchanged after sufficient iterations. A structured grid is adopted, with refinement applied to regions featuring complex flow structures. It is assumed that the gas flows uniformly from the nozzle in the form of an inviscid ideal gas and the parameters conform to the isentropic flow relationship.

We first conduct domain size independence. Table 1 presents the Mach stem height of different domain size for the condition $M_{0}=3$ , $ \textit{NPR}=7$ and $L_{w}=0$ . Group A analyses the effects of the computational domain dimensions perpendicular to the flow direction ( $L_{1}$ ). When replacing the computational domain size from $10H_A\times 30H_A$ to $10H_A\times 15H_A$ , the Mach stem height decreases by only 1.26 %. Group B analyses the effects of the computational domain dimensions along the flow direction ( $L_{2}$ ). When changing from $20H_A\times 15H_A$ to $10H_A\times 15H_A$ , the increase is merely 1.29 %. Therefore, based on a consideration of accuracy and computational efficiency, $L_{1}$ is set as $10H_A$ and $L_{2}$ is set as $15H_A$ for subsequent computations.

Table 1.

Computational domain size independence test results for $M_{0}=3,\, \mathrm{NPR}=7,\, L_{w}=0$ .

Then, we conduct the grid density independence test. Table 2 presents the different grid configurations and corresponding Mach stem height for the condition $M_{0}=3$ , $\textit{NPR}=9$ and $L_{w}=0$ . The number of grid nodes distributed along the Mach stem is N. When replacing the grid from Grid C with Grid B, the Mach stem height difference is only 1.85 %. Moreover, figure 5 comparatively displays the Mach number contours near the triple point. The dash–dot–dotted lines, solid lines and dashed lines represent the computational results for Grid A, Grid B and Grid C, respectively. The Mach number contours of Grid B and Grid C are essentially identical. Thus, we conclude that a minimum of 150 grid points must be distributed along the Mach stem to ensure sufficient accuracy of numerical results.

Table 2.

Grid configurations and computational domain size independence test results for $M_{0}=3, \,\mathrm{NPR}=9,\, L_{w}=0$ . N represents the grid number along the Mach stem.

Figure 5.

CFD simulation results: Mach number contours of overexpanded jet near the upper triple point under $M_{0}=3$ , NPR = 9, $L_{w}$ = 0.

In summary, the grid used in this study is characterised by a computational domain size of $10H_A\times 15H_A$ , with at least 150 grid points guaranteed in the height direction of the Mach stem and a uniformly graded transition adopted in other regions.

2.3. Analysis of flow structures of asymmetric shock reflection in overexpanded jet

Figure 6(a) shows the CFD results for an asymmetric shock reflection under $M_{0}=3$ , NPR = 11, $L_{w}=0.1$ . The NPR exceeds the von Neumann criterion (von Neumann Reference von Neumann1943), resulting in an asymmetric regular reflection. Two incident shocks $A_{1}T$ and $A_{2}T$ form at the lip. These two opposite-side incident shocks reflect as two reflected shocks $B_{1}T$ and $B_{2}T$ . The reflected shocks interact with the jet boundaries $A_{1}B_{1}$ and $A_{2}B_{2}$ , generating two Prandtl–Meyer expansion waves.

Figure 6.

CFD simulation results: Mach number contours (grey solid lines) and characteristic lines (black solid lines) of asymmetric overexpanded jet under $L_{w}=0.1$ . (a) $M_{0}=3$ , NPR = 11. (b) $M_{0}=3$ , NPR = 7.

Figure 6(b) shows the CFD results for an asymmetric shock reflection under $M_{0}=3$ , NPR = 7, $L_{w}=0.1$ . The NPR is in dual solution domain (von Neumann Reference von Neumann1943), resulting in an asymmetric Mach reflection. Two incident shocks $A_{1}T_{1}$ and $A_{2}T_{2}$ are connected by a Mach stem $T_{1}T_{2}$ . At the triple points $T_{1}$ and $T_{2}$ , two incident shocks reflect as two reflected shocks $B_{1}T_{1}$ and $B_{2}T_{2}$ . The region downstream of the Mach stem is initially subsonic and then gradually becomes supersonic, and can be divided into three regions based on whether there is interference from upper/lower reflected expansion waves. Initially, neither the upper nor lower expansion waves interfere with the region downstream of the Mach stem. Subsequently, the lower expansion wave first interferes with the downstream region of the Mach stem. Finally, both the upper and lower expansion waves interfere with the downstream region of the Mach stem. Moreover, the sonic locations $E_{1}$ and $E_{2}$ , and the points of flow horizontalisation are displayed and their positions do not coincide.

3.1. Theoretical model of asymmetric shock reflection in an overexpanded jet

Figure 7 is the schematic illustration of asymmetric Mach reflection in an overexpanded jet. The flow is divided into seven regions: region (0) is upstream of the Mach stem and two incident shock waves; region (1) is between the incident shock wave $i_{1}$ and the reflected shock wave $r_{1}$ ; region (2) is between the reflected shock wave $r_{1}$ and the expansion wave reflected at point $B_{1}$ ; region (3) is downstream of the Mach stem; region (4) is the region of expansion wave reflected at point $B_{1}$ ; region (5) is between the incident shock wave $i_{2}$ and the reflected shock wave $r_{2}$ ; region (6) is between the reflected shock wave $r_{2}$ and the expansion wave reflected at point $B_{2}$ ; and region (7) is the region of the expansion wave reflected at point $B_{2}$ .

Figure 7.

Schematic illustration of asymmetric Mach reflection in an overexpanded jet.

Slipline $s_{1}$ and slipline $s_{2}$ are the upper slipline and lower slipline, respectively. Here, $\delta _{s_1}$ and $\delta _{s_2}$ are respectively the upper slipline deflection angle and the lower slipline deflection angle, and counterclockwise rotation from the positive x-axis is defined as positive here. There is a noticeable misalignment between the upper and lower interference segments of the sliplines. Therefore, when solving sliplines, the problem should be categorised into three cases: (i) both segments being free; (ii) one segment being free, while the other is an interference segment; and (iii) both segments being interference segments. Here, $H_{m}$ is the given Mach stem height. The height of the quasi-one-dimensional channel $H_{s}$ is the height of the quasi-one-dimensional isentropic channel perpendicular to the average flow direction of the throat. Specific expressions are provided in § 3.2.3. Additionally, $M_{s}$ and $p_{s}$ are the Mach number and pressure in region (3); $M _{s_1}$ and $p_{s_1}$ are the Mach number and pressure near the upper slipline $s_{1}$ in region (2) and region (4); $M _{s_2}$ and $p_{s_2}$ are the Mach number and pressure near the lower slipline $s_{2}$ in region (6) and region (7).

3.2. Analytical methods pertinent to flow field structures

3.2.2. Governing equations for slipline flow parameters

Section 2.3 describes the flow field characteristics associated with the free segments of sliplines ( $T_{1}K_{1}$ and $T_{2}K_{2}$ in figure 7). Accompanied by the contraction of the one-dimensional isentropic flow region and the resulting pressure drop, secondary expansion waves are generated along the free segment of the slipline to maintain the pressure balance. These secondary waves further influence the slope of the slipline.

This mechanism was applied by Omelchenko et al. (Reference von Neumann2003) to Mach reflection in an overexpanded jet and by Gao & Wu (Reference Gao and Wu2010) to that in wedge flows, and was systematically discussed by Silnikov & Chernyshov (Reference Silnikov and Chernyshov2015), who further analysed the possible geometric shapes of the free segment and the interference segment of the slipline. Furthermore, Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022) also incorporated this mechanism in their modelling and analysis of shock reflection problems in symmetric jet flows. Furthermore, secondary waves have been accounted for by Bai & Wu (Reference Bai and Wu2017) and Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022).

Therefore, we can present the governing equations for flow parameters near slipline $s_{1}$ in the upper region (2):

(3.1) \begin{align} \left \{ \begin{aligned} &\nu \big(M_{s_1} \big)-\nu \big(M_{2}^{T_1}\big) =-\delta _{s_1}+\delta _{m}^{T_1},\\[6pt] &p_{s_1}=p_{2}^{T_1}\left(\frac {\vartheta \big(M_{2}^{T_1} \big)}{\vartheta \big(M_{s_1}\big)}\right)^{\frac {\gamma }{\gamma -1}},\\ \end{aligned} \right . \end{align}

where $M_{2}^{T_1}, p_{2}^{T_1}$ is the Mach number and pressure of the airflow in region (2) near the triple point $T_{1}$ ; $\nu (M)$ is the Prandtl–Meyer relation and

(3.2) \begin{align} \left \{ \begin{aligned} &\vartheta (M)=1+\frac {\gamma -1}{2}M^2,\\[5pt]&\nu (M)=\sqrt {\frac {\gamma +1}{\gamma -1}}\arctan \sqrt {\frac {\gamma -1}{\gamma +1}(M^2-1)}-\arctan \sqrt {M^2-1}.\\ \end{aligned} \right . \end{align}

Additionally, the governing equations for flow parameters near slipline $s_{2}$ in region (6) are

(3.3) \begin{align} \left \{ \begin{aligned} &\nu \big(M_{s_2}\big)-\nu \big(M_{6}^{T_2}\big) = \delta _{s_2}-\delta _{m}^{T_2},\\[5pt]&p_{s_2}=p_{6}^{T_2} \left(\frac {\vartheta \big(M_{6}^{T_2}\big)}{\vartheta \big(M_{s_2}\big)}\right)^{\frac {\gamma }{\gamma -1}},\\ \end{aligned} \right . \end{align}

where $M_{6}^{T_2}, p_{6}^{T_2}$ are the Mach number and pressure of the airflow in region (6) near the triple point $T_{2}$ .

The flow characteristics of the interference segments of the sliplines ( $K_{1}E_{1}$ and $K_{2}E_{2}$ in figure 7) are also described in § 2.3. The pressure drop resulting from the flow passing through the reflected expansion waves is the dominant factor driving the pressure change. To match this pressure drop, the quasi-one-dimensional isentropic flow region contracts, ultimately leading to a change in the curvature of the slipline.

In the solution for a symmetric jet, Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022) employed the straight-line assumption to describe the reflected expansion waves and solved for the interference segment of the slipline. In this study, we further assume the reflected expansion waves to be small-curvature curves, based on which the governing equations for flow parameters near the interference segments of the sliplines are established. This assumption aligns with the flow characteristics depicted in figure 6. The derivation of the governing equations for flow parameters in the slipline interference segments can refer to the methodology of Bai (Reference Bai2023), with the distinction that in this model, the expansion waves do not interact with the reflected shock waves.

The expressions for solving the small flow deflection increment $\mathrm{d}\theta _{B_1}$ at position $B_1$ in region (4) and the angular change $\mathrm{d}\delta _{s_1}$ of slipline $s_{1}$ are

(3.4) \begin{align} \left \{ \begin{aligned} &\mathrm{d}\theta _{B_1}=\frac {L_{21}\varPhi _{B_1s_1}\mathrm{d}p_{s_1}+L_{31}}{L_{11}+L_{21}\big(1+\varPhi _{B_1s_1}p_{B_1}\varPsi _{p}\big(M_{B_1}\big)\big)},\\[6pt]&\mathrm{d}\delta _{s_1}=\big(1+\varPhi _{B_1s_1}p_{B_1}\varPsi _{p}\big(M_{B_1}\big)\big)\mathrm{d}\theta _{B_1}-\varPhi _{B_1s_1}\mathrm{d}p_{s_1},\\ \end{aligned} \right . \end{align}

and the expressions for solving the small flow deflection increment $\mathrm{d}\theta _{B_2}$ at position $B_2$ in region (7) and the angular change $\mathrm{d}\delta _{s_2}$ of slipline $s_{2}$ are

(3.5) \begin{align} \left \{ \begin{aligned} &\mathrm{d}\theta _{B_2}=\frac {L_{22}\varPhi _{B_2s_2}\mathrm{d}p_{s_2}+L_{32}}{L_{12}+L_{22} \big(1+\varPhi _{B_2s_2}p_{B_2}\varPsi _{p}\big(M_{B_2}\big)\big)},\\[6pt] &\mathrm{d}\delta _{s_2}= \big(1+\varPhi _{B_2s_2}p_{B_2}\varPsi _{p}\big(M_{B_2}\big)\big)\mathrm{d}\theta _{B_2}-\varPhi _{B_2s_2}\mathrm{d}p_{s_2},\\ \end{aligned} \right . \end{align}

where the $\mathrm{d}x_{s_1}$ , $\mathrm{d}y_{s_1}$ , $\mathrm{d}x_{s_2}$ , $\mathrm{d}y_{s_2}$ are the coordinate increments along the slipline $s_{1}$ and $s_{2}$ , $\mathrm{d}y_{s_1}=-\tan \delta _{s_1}\,\mathrm{d}x_{s_1}$ , $\mathrm{d}y_{s_2}=\tan \delta _{s_2}\,\mathrm{d}x_{s_2}$ ; $M_{B_1}$ and $p_{B_1}$ are the Mach number and pressure in region (4) near $B_{1}$ ; $M_{B_2}$ and $p_{B_2}$ are the Mach number and pressure in region (7) near $B_{2}$ ; $\theta _{11}$ , $\theta _{21}$ are the angles between the expansion wave and the x-axis at $B_{1}$ and near the slipline $s_1$ ; $\theta _{12}$ , $\theta _{22}$ are the angles between the expansion wave and the x-axis at $B_{2}$ and near the slipline $s_2$ ; for the detailed definitions of $L_{11}$ , $L_{21}$ , $L_{31}$ , $\varPhi _{B_1s_1}$ , $L_{12}$ , $L_{22}$ , $L_{32}$ , $\varPhi _{B_2s_2}$ , $ \varPsi _{p}(M)$ , refer to (B7), (B5), (B7). The detailed mathematical derivations can be found in Appendix B.

Moreover, the Mach number increments near $s_{1}$ can be obtained by

(3.6) \begin{align} \begin{aligned} \mathrm{d}M_{s_1}=\frac {\varPsi _{M}\big(M_{s_1}\big)}{\varPsi _{p}\big(M_{s_1}\big)}\frac {\mathrm{d}p_{s_1}}{p_{s_1}} \end{aligned} \end{align}

and the Mach number increments near $s_{2}$ can be obtained by

(3.7) \begin{align} \begin{aligned} \mathrm{d}M_{s_2}=\frac {\varPsi _{M}\big(M_{s_2}\big)}{\varPsi _{p}\big(M_{s_2}\big)}\frac {\mathrm{d}p_{s_2}}{p_{s_2}} .\end{aligned} \end{align}

3.2.3. Iterative solution method and procedure for sliplines

For the iterative solution of the sliplines, Lin et al. (Reference Lin, Bai and Wu2019) assumed that the average flow deflection angle in region (3) $\delta _{ave}$ was equal to the average of the flow deflection angles at triple points $T_{1}$ and $T_{2}$ of region (3). However, Tao et al. (Reference Tao, Liu, Fan, Xiong, Yu and Sun2017) assumed that $\delta _{ave}$ continually updates throughout the iteration, equalling the average of the upper/lower slipline angles $\delta _{s_1}$ and $\delta _{s_2}$ .

Figure 8.

Schematic illustration of sliplines’ iterative solution.

The sliplines’ iterative solution presented in this study assumes that the value of $\delta _{ave}$ remains equal to the average of the upper and lower slipline angles, and continuously updates throughout the iteration process. By integrating the quasi-one-dimensional isentropic flow theory, the pressure equilibrium condition across the slipline and the governing equations for flow parameters near the sliplines, a complete solution for the sliplines is achieved.

Figure 8 is a schematic illustration of the geometric relationship between the upper and lower sliplines from step $i-1$ to step $i+1$ during the iterative solution process. During the advancement, this model assumes that the line segment $s_{1}s_{2}$ is always perpendicular to the mean flow direction.

Taking the case where $s_{2}$ first interacts with the expansion wave as an example, the detailed solution steps are as follows.

1. (i) Given $s_{1,i-1}$ and $s_{2,i-1}$ , and the upper/lower slipline angles $\delta _{s_1,i-1}$ and $\delta _{s_2,i-1}$ , compute $s _{ave,i}$ and the angle $\delta _{ave,i-1}$ :

   (3.8) \begin{align} \left \{ \begin{aligned} &x_{s_{ave}, i}=x_{s_{ave}, i-1}+\mathrm{d}x_{ave},\\[3pt]&y_{s_{ave}, i}=y_{s_{ave}, i-1}+\tan \big(\delta _{ave,i-1}\big)\,\mathrm{dx}_{ave},\\[3pt]&\delta _{ave, i-1}=\frac {1}{2}\big(\delta _{s_1, i-1}+\delta _{s_2, i-1}\big). \end{aligned} \right . \end{align}
   Additionally, at the start of the slipline advancement and solution process, $x_{s_{ave},1}={1}/{2}(x_{T_1}+x_{T_2})$ , $y_{s_{ave},1}={1}/{2}(y_{T_1}+y_{T_2})$ and $\delta _{ave}={1}/{2}(-\delta _m^{T_{1}}+\delta _m^{T_{2}})$ .

2. (ii) Solve the $s_{2, i}$ and $s_{1, i}$ . The equation of the straight $s_{1,i}s_{2,i}$ is

   (3.9) \begin{align} \begin{aligned} &y-y_{s_{ave}, i}=-\frac {1}{\tan \big(\delta _{ave, i-1}\big)}\big(x-x_{s_{ave}, i}\big)\\ \end{aligned} \end{align}
   and the equations of straight $s_{1,i-1}s_{1,i}$ , $s_{2,i-1}s_{2,i}$ are
   (3.10) \begin{align} \left \{ \begin{aligned} &y-y_{s_{1}, i-1}=-\tan \delta _{s_1,i-1}(x-x_{s_{1}, i-1}),\\[5pt]&y-y_{s_{2}, i-1}=\tan \delta _{s_2,i-1}(x-x_{s_{2}, i-1}).\\ \end{aligned} \right . \end{align}
   By solving (3.9) and (3.10) simultaneously, the coordinates of $s_{2, i}$ and $s_{1, i}$ can be determined and $\mathrm{d}x_{s_{2}, i}=x_{s_{2}, i}-x_{s_{2}, i-1}$ , $\mathrm{d}x_{s_{1}, i}=x_{s_{1}, i}-x_{s_{1}, i-1}$ .

3. (iii) Solve the $H_{s,i}$ and $M_{s,i}$ , $p_{s,i}$ . Here, $H_{s,i}$ is the height of flow channel and $H_{s,i}=(y_{s_1,i}-y_{s_2,i})/\cos \delta _{ave,i-1}$ ; $M_{s,i}$ and $p_{s,i}$ follow the quasi-one-dimensional isentropic flow theory relation

   (3.11) \begin{align} \left \{ \begin{aligned} &\frac {H_{s,i}}{H_{m}}=\frac {M_{m}}{M_{s,i}}\left(\frac {\vartheta \big(M_{s,i}\big)}{\vartheta (M_{m})}\right)^{( {\gamma +1}/{2(\gamma -1))}}, \\&p_{s,i}=p_{m}\left(\frac {\vartheta (M_{m})}{\vartheta \big(M_{s,i}\big)}\right)^{({\gamma }/{\gamma -1})} ,\\ \end{aligned} \right . \end{align}
   where $M_{m}$ is the average Mach number and can be expressed as the ratio of the average velocity to the average sonic speed (Li & Ben-Dor Reference Li, He, Zhang, Hao, Wu and Liu1997a ,Reference Li, He, Zhang, Hao, Wu and Liu b )
   (3.12) \begin{align} \begin{aligned} &M_m=\frac {2\big(\rho _3^{T_1}u_3^{T_1}\cos \delta _m^{T_1}+\rho _Gu_G\big)}{\big(\rho _3^{T_1}+\rho _G \big)\big(a_3^{T_1}+a_G \big)},\\ \end{aligned} \end{align}
   where $u_3^{T_1}=a_3^{T_1}M_3^{T_1}$ and $u_G=a_GM_G$ ; $\rho _3^{T_1}$ , $a_3^{T_1}$ , $M_3^{T_1}$ are the density, sonic speed and Mach number at region (3) near $T_1$ ; $\rho _G$ , $a_G$ $M_G$ can be obtained using the normal shock relations. Here, $p_{m}$ is the average pressure behind the Mach stem.

4. (iv) Solve $\delta _{s_1,i}$ and $\delta _{s_2,i}$ . Based on the pressure equilibrium relation across the slip lines, it follows that $p_{s_1,i}=p_{s,i}=p_{s_2,i}$ or $\mathrm{d}p_{s_1,i}=(p_{s,i}-p_{s,i-1})=\mathrm{d}p_{s_2,i}$ . The specific solution methods are as follows: when both $s_{1}$ and $s_{2}$ are in the free segment, $\delta _{s_1,i}$ and $\delta _{s_2,i}$ can be solved using (3.1) and (3.3); when $s_{2}$ is in the interference segment and $s_{1}$ is in the free segment, $\delta _{s_1,i}$ can be obtained from (3.1), $\mathrm{d}\delta _{s_2,i-1}$ is derived from (3.5) and $\delta _{s_2,i}=\delta _{s_2,i-1}+\mathrm{d}\delta _{s_2,i-1}$ ; when both $s_{1}$ and $s_{2}$ are in the interference segment, $\mathrm{d}\delta _{s_1,i-1}$ is derived from (3.4) and $\delta _{s_1,i}=\delta _{s_1,i-1}+\mathrm{d}\delta _{s_1,i-1}$ , $\mathrm{d}\delta _{s_2,i-1}$ is derived from (3.5) and $\delta _{s_2,i}=\delta _{s_2,i-1}+\mathrm{d}\delta _{s_2,i-1}$ .

5. (v) Solve $s_{ave, i+1}$ , the average flow angle $\delta _{ave, i}$ , then repeat steps (i)–(iv) to determine $s_{1, i+1}$ , $\delta _{s_1, i+1}$ , $s_{2, i+1}$ and $\delta _{s_2, i+1}$ .

3.3. Global solution based on sonic throat compatibility condition

Bai & Wu (Reference Bai and Wu2017) proposed the sonic throat compatibility condition that determines the Mach stem height for classical shock reflection. This condition states that, if the slipline angle is zero $\delta _{s}=0$ at the exact location where the sonic throat condition occurs $M_{s}=M_{s}^{+}=1$ , then this solution is compatible. Lin et al. (Reference Lin, Bai and Wu2019) subsequently refined this criterion by proposing a throat compatibility condition for asymmetric cases, indicating that the sonic throat ( $E_{1}E_{2}$ ) locates at the position where the angles of the upper and lower slip lines are equal.

As noted in § 2.3, for an asymmetric overexpanded jet, the sonic throats do not coincide with the points of flow horizontalisation. Therefore, in this work, we adopt the throat compatibility condition given by Lin et al. (Reference Lin, Bai and Wu2019):

(3.13) \begin{align} \left \{\begin{aligned} &\delta _{s_1}^{E_1} - \delta _{s_2}^{E_2} = 0, \\[3pt] & M_{s_1}^{E_1} = M_{s_2}^{E_2}=1, \\[3pt] &\frac {y_{E_2} - y_{E_1}}{x_{E_2} - x_{E_1}} = -\frac {1}{\tan \delta _{\text{ave}}}, \end{aligned} \right . \end{align}

which means that at the throat location $E_{1}$ and $E_{2}$ , the Mach number $M_{s_1}^{E_1}$ and $M_{s_2}^{E_2}$ equal 1, and the flow deflection angles of the two sliplines $\delta _{s_1}^{E_1}$ and $\delta _{s_2}^{E_2}$ are equal.

Figure 9.

Schematic illustration of global solution flowchart.

According to the sonic throat compatibility condition for the global model of flow field structures in asymmetric overexpanded jet, the iterative solution flowchart is presented as figure 9, with the specific steps as follows.

1. (i) Solve the flow parameters near the triple points $T_{1}\,\text{and}\,T_{2}$ , the triple-shock relation can simply estimate the average values of flow parameters. The pressure, Mach number, temperature near the triple points can be solved, and the shock angles and flow deflection angles of the incident shocks, reflected shocks and jet boundaries can also be obtained.

2. (ii) According to the nozzle geometry, a solution interval for the Mach stem height is preset as $(H_{\textit{min}},H_{\textit{max}})$ . At the beginning of the solution process, the upper and lower bounds of the solution interval are preset based on the nozzle geometry, as illustrated in figure 10. The initial value of $H_{\textit{min}}$ is 0 and the initial value of $H_{\textit{max}}$ is $H_{A}(1-L_{w}\tan \beta _{i_1})$ . During the iterative solution process, the Mach stem height solution interval is continuously updated, while the model-predicted Mach stem height remains 0.5 ( $H_{\textit{min}}+H_{\textit{max}}$ ).

3. (iii) Solve the coordinates of triple points $T_{1}$ and $T_{2}$ . Based on the geometric configuration among the incident shock wave, reflected shock wave and Mach stem, the triple point coordinates are determined. Treating the value of $H_{m}$ obtained from step (ii) as a known quantity, the system yields seven fundamental equations (i.e. (A1)–(A5) in Appendix A) with seven unknowns during the triple point coordinate calculation, thereby admitting a unique solution.

4. (iv) Solve the sliplines. The detailed iterative procedure for solving the sliplines has been provided in § 3.2.3. During the solution process, attention should be paid to identifying the positions of $K_{1}$ and $K_{2}$ to select the appropriate governing equations.

5. (v) Update the Mach stem solution interval using the sonic throat compatibility condition. When computing the Mach stem height $H_{m}$ in step (ii), the shape of the slipline will be uniquely given whether the solution satisfies the sonic throat compatibility condition or not. However, only solutions that meet the sonic throat compatibility conditions are physical. Therefore, during the solution process, the Mach stem solution interval should be continuously updated based on the following conditions until the throat compatibility condition is met. Additionally, the updated conditions are that if $\delta _{s_1}^{E_1}-\delta _{s_2}^{E_2}\neq 0,M_{s_1}=M_{s_2}=1$ , then we make $H_{\textit{min}}=H_{m}$ and if $\delta _{s_1}^{E_1}-\delta _{s_2}^{E_2}=0,M_{s_1}=M_{s_2}\neq 1$ , then we make $H_{\textit{max}}=H_{m}$ . In this way, steps (i)–(v) are repeated until the solution to the subsonic region corresponding to the accuracy requirement or the height of Mach stem just meets the sonic throat compatibility condition.

Figure 10.

Schematic illustration of initial interval $(0,H_{A}(1-L_{w}\tan \beta _{i_1}))$ of Mach stem solution.

4.1. Accuracy of the theoretical model for symmetric Mach reflection in an overexpanded jet

Figure 11 presents Mach stem height predictions of current model and CFD results, comparing model predictions and numerical results proposed by Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022) under different NPR conditions at $M_{0}=2.75$ . Significant discrepancies exist between the numerical results of this study and those reported by Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022). In their study, Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022) performed computations for overexpanded jet shock reflection under symmetric conditions using a 400 $\times$ 200 structured grid within a computational domain of $4H_{A}\times 2H_{A}$ . Notably, the uniform grid distribution employed proved insufficient for resolving small Mach stems. When the Mach stem height spans a few cells, it becomes easily overlooked. This may explain why their numerical results indicate the disappearance of the Mach stem before reaching the von Neumann point. In contrast, our numerical simulations were conducted on a significantly larger computational domain of $10H_{A}\times 15H_{A}$ , while maintaining adequate grid resolution along the Mach stem direction. Particularly for cases with higher NPR values (corresponding to smaller Mach stem height), we ensured a minimum of 150 grid points across the Mach stem height to enhance capture accuracy. This enables our numerical calculations to successfully capture the Mach stem at $M_{0}=2.75$ under NPR conditions of 7, 7.2 and 7.4 (prior to reaching the von Neumann point), thereby demonstrating the reliability of our computational approach. Consequently, subsequent comparisons will focus on comparing model predictions with the numerical results from this study. Though all theoretical models accurately predict the NPR corresponding to the von Neumann condition, only our model successfully captures the trend of an increasing slope with NPR and predicts significant Mach stem height even at very low NPR conditions. In general, the Mach stem height predictions from our model demonstrate the best agreement with our numerical results.

Figure 11.

Comparison among current CFD results, current theoretical model results and CFD results, and theoretical model results of Paramanantham et al. (Reference Paramanantham, Janakiram and Gopalapillai2022) under $M_{0}$ = 2.75.

Table 3.

Comparison of Mach stem height between current model and CFD for symmetric Mach reflection in an overexpanded jet.

Figure 12.

Comparison between current model and CFD under $M_{0}$ = 3, NPR = 7.

A comparison of the CFD results with the predictions of our model is presented in figure 12 under the conditions of $M_{0}=3$ and NPR = 7. The predicted Mach stem height deviation (the relative error of the theoretical model with respect to the CFD results) is 3.70 $\%$ . Mach number contours obtained by CFD are shown with the solid lines, while the dashed lines depict the model predictions. The predicted triple points and reflected shock waves exhibit a slight upstream shift compared with the CFD results. In general, the incident/reflected shock waves, jet boundaries and sliplines predicted by the current model show good agreement with the CFD results.

Table 3 presents further results of the Mach stem height obtained from the theoretical model and CFD, along with the corresponding deviations, under various nozzle exit conditions. The model-predicted Mach stem height is slightly larger than the CFD results. Maximum deviation of the model-predicted Mach stem height is 4.48 % for $M_{0}$ = 3. Maximum deviation of the model-predicted Mach stem height is 4.74 % for $M_{0}$ = 4. In general, the model demonstrates high accuracy in predicting the Mach stem height for symmetric shock reflection.

Figure 13.

Influence analysis of individual nozzle exit parameters. (a) Influence of NPR on Mach stem height under different $M_{0}$ . (b) Influence of $M_{0}$ on Mach stem height under different NPR.

Additionally, figure 13 provides further results on the influence of the Mach stem height with nozzle exit parameters as predicted by the model. Figure 13(a) shows the influence of Mach stem height with NPR under different $M_{0}$ conditions. The Mach stem height decreases as NPR increases and the rate of decrease accelerates with larger NPR. Furthermore, in logarithmic coordinates, as $M_{0}$ increases, the curve between NPR and Mach stem height becomes increasingly flat. Figure 13(b) shows the influence of Mach stem height with $M_{0}$ under different NPR conditions. The Mach stem height increases as $M_{0}$ increases and the rate of increase decelerates with larger NPR. Furthermore, under different NPR conditions, the curves between $M_{0}$ and Mach stem height are approximately parallel.

In summary, the numerical method employed in this paper provides more accurate capture of flow field structures in symmetric jet Mach reflection. Additionally, our model demonstrates an excellent predictive capability for the flow field structure. The maximum deviation between model-predicted Mach stem height and numerical results does not exceed 4.74 % for both the $M_{0}=3$ and $M_{0}=4$ cases. Additionally, the Mach stem height decreases as NPR increases and $M_{0}$ decreases.

4.2. Accuracy of the theoretical model for asymmetric Mach reflection in an overexpanded jet

Figure 14 presents a comparison between the current model and CFD under $M_{0}$ = 2.75, NPR = 6 and $L_{w}$ = 0.1. The streamline passing through y = 0 is depicted. Due to the geometric asymmetry of the nozzle ( $L_w$ = 0.1), the flow is deflected downward, resulting in streamline deviation towards the lower region. The predicted Mach stem height deviation is 2.76 $\%$ . The Mach number contour obtained by CFD is shown with the solid lines, while the dashed lines depict the model predictions. The predicted triple points exhibit good agreement compared with the CFD results. In general, the incident/reflected shock waves, jet boundaries and sliplines predicted by the model show good agreement with the CFD results.

Figure 14.

Comparison between current model and CFD under $M_{0}$ = 2.75, NPR = 6 and $L_{W}$ = 0.1.

Table 4 presents the height of the Mach stem obtained by the current model, CFD and deviations under different nozzle exit conditions. The model-predicted Mach stem height remains slightly larger than the CFD results. Compared with the prediction accuracy for symmetric Mach reflection, the model’s accuracy for asymmetric Mach reflection shows a slight decrease. Maximum deviation of the model-predicted Mach stem height does not exceed 3.82 % for $M_{0}=3$ and NPR = 7. The maximum deviation of the model-predicted Mach stem height does not exceed 4.07 % for $M_{0}=4$ and NPR = 16. In general, the model demonstrates high accuracy in predicting the Mach stem height for asymmetric jet flow. Additionally, both CFD and current model results consistently demonstrate that the Mach stem height decreases as $L_{w}$ increases.

Table 4.

Comparison of Mach stem height between current model results and CFD results for an asymmetric jet.

When the $M_{0}$ and $\mathrm{NPR}$ are given, the triangles $A_{1}B_{1}T_{1}$ and $A_{2}B_{2}T_{2}$ (figure 7) formed by the incident shocks, reflected shocks and jet boundaries above and below are similar. Changing the unilateral $L_{W}$ actually alters the geometric size of the corresponding triangle. When $L_{W}$ increases, the lower triangle $A_{2}B_{2}T_{2}$ continuously decreases in size. This change is transmitted to the slipline via expansion waves, altering the geometric shape of the boundary of the subsonic region downstream of the Mach stem (i.e. the slipline). With the nozzle exit flow parameters remaining unchanged, the shock angle does not change; instead, the Mach stem moves in the downstream direction, ultimately resulting in a reduction in its height.

4.3. Influence of nozzle exit conditions on Mach stem height

For asymmetric Mach reflection in an overexpanded jet, there is an additional parameter $L_w$ that affects the Mach stem height comparing with the symmetric Mach reflection. Referring to the theoretical formula for the parameters influencing the Mach stem height for symmetric Mach reflection (Paramanantham et al. Reference Paramanantham, Janakiram and Gopalapillai2022), we can provide the theoretical formula for asymmetric Mach reflection as

(4.1) \begin{align} \begin{aligned} \frac {H_{m}}{H_{A}}=f(M_{0}, \mathrm{NPR},L_{w}, \gamma ). \end{aligned} \end{align}

Figure 15.

Influence analysis of individual nozzle exit parameters. (a) Influence of NPR on Mach stem height with increasing $L_{w}$ under $M_{0}=3$ . (b) Influence of $M_{0}$ on Mach stem height with increasing $L_{w}$ under NPR = 7.

Figure 16.

Influence analysis of individual nozzle exit parameters. (a) Influence of NPR on the effect of $L_w$ reduction on Mach stem height under $M_{0}=3$ . (b) Influence of NPR on Mach stem height with different $M_{0}$ under $L_{w}=0.1$ .

Figure 17.

Influence analysis of individual nozzle exit parameters. (a) Influence of $M_{0}$ on the effect of $L_w$ reduction on Mach stem height under NPR = 7. (b) Influence of $M_{0}$ on Mach stem height with different NPR under $L_{w}=0.1$ .

Next, we will analyse the effects of several parameters on the Mach stem height. The model-predicted results under different nozzle exit conditions are presented in figures 15, 16 and 17.

Figure 15(a) presents the variation of Mach stem height with $L_{w}$ under different NPR and $M_{0}=3$ . Mach stem height slightly decreases as $L_{W}$ increases. The model predicts that the Mach stem height at $L_w = 0.2$ is 9.64 % lower than that at $L_w = 0$ for NPR = 7. In comparison, the model predicts that the Mach stem height at $L_w = 0.2$ is 10.23 % lower than that at $L_w = 0$ for NPR = 9. This indicates that an increase in NPR enhances the reducing effect of $L_{w}$ on the Mach stem height.

Figure 15(b) presents the variation of Mach stem height with $L_{w}$ under different $M_{0}$ and NPR = 7. The model predicts that the Mach stem height at $L_w = 0.2$ is 12.14 % lower than that at $L_w = 0$ for $M_{0}$ = 3.1. In comparison, the model predicts that the Mach stem height at $L_w = 0.2$ is 11.42 % lower than that at $L_w = 0$ for $M_{0}$ = 2.7. This indicates that an increase in $M_{0}$ enhances the reducing effect of $L_{w}$ on the Mach stem height.

Figure 16(a) presents the variation of Mach stem height with NPR under different $L_{w}$ and $M_{0}=3$ . The Mach stem height decreases as NPR increases. The model predicts that the Mach stem height at NPR = 9 is 50.08 % lower than that at NPR = 7 for $L_w = 0$ . In comparison, the model predicts that the Mach stem height at NPR = 9 is 48.69 % lower than that at NPR = 7 for $L_w = 0.3$ . This indicates that an increase in $L_{w}$ weakens the reducing effect of NPR on the Mach stem height.

Figure 16(b) presents the variation of Mach stem height with NPR under different $M_{0}$ and $L_{w}=0.1$ . When compared with the results of the symmetric case ( $L_{w}=0$ ) in figure 13(a), the curve forms remain largely consistent. Only a slight decrease in Mach stem height is observed, ultimately decreasing to zero at the NPR corresponding to the von Neumann condition.

Figure 17(a) presents the variation of Mach stem height with $M_{0}$ under different $L_{w}$ and NPR = 7. The Mach stem height increases as $M_{0}$ increases. The model predicts that the Mach stem height at $M_{0}=3$ is 205.55 % higher than that at $M_{0}=2.7$ for $L_w = 0$ . In comparison, the model predicts that the Mach stem height at $M_{0}=3$ is 201.80 % higher than that at $M_{0}=2.7$ for $L_w = 0.3$ . This indicates that an increase in $L_{w}$ weakens the enhancement effect of $M_{0}$ on the Mach stem height.

Figure 17(b) presents the variation of Mach stem height with $M_{0}$ under different NPR and $L_{w}=0.1$ . Similarly, when compared with the symmetric case ( $L_{w}=0$ ) results in figure 13(b), the curve forms remain largely consistent. Only a slight decrease in Mach stem height is observed and the Mach stem height increases with the Mach number.

In summary, an increase in $L_w$ weakens both the reducing effect of NPR and the enhancement effect of $M_{0}$ on the Mach stem height. Furthermore, the reducing effect of $L_{w}$ is slightly enhanced by the increase in NPR and the increase in $M_{0}$ .

4.4. Critical geometric condition for complete misalignment of the interference segment of the upper and lower sliplines

As $L_{w}$ gradually increases, the flow field structure will evolve from figure 14 into the form shown in figure 18(a), which presents the Mach number contour under $L_{w}$ = 0.3, $M_{0}$ = 4 and NPR = 18. It can be observed that the interference segment of upper and lower sliplines become completely misaligned, which means the last expansion wave from one side interacts with the slipline, while the first expansion wave from the opposite side does not interact with the slipline. The Mach number contour line of 1 is marked and the sonic throat location shifts.

A critical flow field structure can be defined, schematically illustrated in figure 18(b). This critical state occurs when the last expansion wave generated by the reflected shock on one side ( $B_2F_2$ in figure 18 b) precisely impinges on the slipline, while the first expansion wave from the opposite reflected shock ( $B_1K_1$ in figure 18 b) just begins to interact with the slipline. The theoretical model used in this study can be employed to predict $L_{w}$ for this condition, which is termed $L_{w}^{\textit{crit}}$ .

Figure 18.

(a) Mach number contours (grey solid lines) and characteristic lines (black solid lines) under $M_{0}$ = 4, NPR = 18, $L_{w}=0.3$ . (b) Schematic illustration of the flow field structures of the complete misalignment of the interference segment of upper/lower sliplines.

Figure 19.

Influence of NPR on $L_{w}^{\textit{crit}}$ at different $M_{0}$ .

Figure 19 is the prediction of corresponding $L_{w}^{\textit{crit}}$ at different NPR when given $M_{0}$ . For lower nozzle exit Mach numbers ( $M_{0}$ = 2.75 and $M_{0}$ = 3), the curves show a slight upward convexity. Before the inflection point, $L_{w}^{\textit{crit}}$ increases rapidly with NPR; after the inflection point, the growth rate slows down, which indicates that the enhancement effect of the NPR on $L_{w}^{\textit{crit}}$ gradually weakens under low-Mach-number conditions. In contrast, for higher Mach numbers ( $M_{0}$ = 4, 5, 6), the curves show a slight concavity downward. Before the inflection point, $L_{w}^{\textit{crit}}$ increases slowly with NPR; after the inflection point, the growth rate accelerates, which indicates that the enhancement effect of the NPR on $L_{w}^{\textit{crit}}$ gradually intensifies under high-Mach-number conditions. For a given NPR, $L_{w}^{\textit{crit}}$ decreases as $M_0$ increases. At $ \textit{NPR}=6.5$ , $L_{w}^{\textit{crit}}$ is 0.26 for $M_0=2.75$ and 0.18 for $M_0=3$ , representing a decrease of 30.77 % when $M_0$ increases from 2.75 to 3. For a given Mach number, $L_{w}^{\textit{crit}}$ increases with increasing NPR. At $M_0=4$ , $L_{w}^{\textit{crit}}$ is 0.14 for $\mathrm{NPR}=13$ and 0.41 for $\mathrm{NPR}=27.5$ , representing an increase of 192.86 % when NPR increases from 13 to 27.5.

The trend of how $L_{w}^{\textit{crit}}$ changes with NPR and $M_{0}$ can be easily understood by considering the degree of overexpansion. When other conditions remain constant, decreasing the Mach number or increasing the NPR both lead to a reduction in the intensity of the incident shock wave. In conditions of lower overexpansion (i.e. at lower NPR), the weakened incident shock wave intensity signifies a decrease in the shock wave angle. When the shock wave angle of the incident shock wave decreases, it geometrically results in an increase in the streamwise length of the interaction segments (interference segments) between expansion waves and sliplines. However, increasing $L_{w}$ actually reduces the length of the lower interference segment and shifts it upstream. When the upper interference segment lengthens, the lower interference segment has more space to move upstream, thereby achieving a larger $L_{w}^{\textit{crit}}$ .