An Approximate Analytical Model of a Jet Flow with Mach Reflection and Pulsed Energy Supply at the Main Shock

Abstract

The supersonic flow of a reactive gas mixture with Mach reflection of oblique shocks and pulsed energy supply at the Mach stem is considered within the framework of the Chapman–Jouguet theory. An approximate analytical model is proposed that quickly determines the shape and size of the shock-wave structure as well as the flow parameters in various flow regions. As an example of the application of the proposed analytical model, the “first barrel” of a highly overexpanded jet flow of an air-methane mixture with a high supersonic velocity, is studied. Flows of hydrogen–air and hydrogen–oxygen mixtures were also considered for comparison with preceding numerical results. The height of the triple point of the Mach reflection is determined in the presence of a change in the chemical composition of the mixture and an isobaric pulsed energy supply at the main shock.

1. Introduction

The analysis of flows with an irregular (Mach) reflection of shock waves is important in the gas-dynamic design and optimization of supersonic air intakes, nozzles, jet apparatus, launch complexes, as well as in the development of blast technologies and means of suppressing the damaging effect of an explosion [1]. According to [1,2,3], various flow parameters (stagnation pressures, velocities, dynamic pressures, etc.) in flow regions separated by a slipstream issuing from the triple point of Mach reflection can differ significantly. It can determine the design of gas-dynamic devices. For example, the idea of a new combined ramjet engine was proposed and theoretically substantiated in [4,5,6]. According to this concept, the gas flow behind the reflected shock $j_2$ (in region 2 in Figure 1a), which has a significantly higher total pressure than the flow behind the main (Mach) shock $j_3$ (in region 3), can be used in the thermodynamic cycle of a classic ramjet engine. At the same time, the flow behind the main shock has a significantly higher temperature, especially at high supersonic speeds, which can initiate the detonation of the gas mixture, thus it can be used in the thermodynamic cycle of a detonation engine. For the successful separation of two streams behind the triple point T of the Mach reflection, it is important to determine the height triple point $y_T$ and shape $y\left(x\right)$ of the slipstream $\tau$ that emanates from it, as well as the shape and size of all other discontinuities in the flowfield.

One of the first approximate analytical models of a planar supersonic flow with Mach reflection was formulated in [7]. Within the framework of this model, it was reasonably assumed [8] that the flow in region three forms a so-called “virtual nozzle” with subsonic flow acceleration behind the shock $j_3$ up to the critical sound speed ($M_3=1$) in the narrowest section (CC’ in Figure 1b). The acceleration to the critical speed must coincide with the turn of the slipstream $\tau$ in the horizontal direction (the angle of flow $\theta =0$ in point C) under the influence of a rarefaction wave ${\psi }_4$ that falls from the trailing edge of the wedge (Figure 1a), or from the jet boundary (Figure 1b), or formed in some other way. However, it was unreasonably assumed in [7] that the slipstream $\tau$ is rectilinear, and the critical section of the “virtual nozzle” three corresponds to the point of incidence of the first characteristic BD of the rarefaction wave ${\psi }_4$. It led to large (up to 50–90%) errors in determining the height of the triple point and other flow parameters. The application of the Grib–Ryabinin model [9] to region two of the supersonic flow, carried out in [10], led to a slight refinement of the results with a simultaneous complication of the mathematical model.

The analytical model proposed in [11] and detailed in [12,13] for flows in supersonic jets and narrowing channels considers flow region two as a simple Prandtl–Meyer rarefaction wave conjugated with region three by the condition of equality of pressures on the sides of an upwardly convex slipstream $\tau$ [14]. Taking into account the variation in the angle $\theta$ of the inclination of the slipstream significantly improves the accuracy of calculations of the supersonic part of the flow (as shown in [11,12,13], the error in determining the height of the triple point is about 0.5–2% compared to the results of calculations of the supersonic part of the flow by the method of characteristics). The opposite turn of the slipstream is carried out on its final segment DC. It can be based on the solution of the problem of the incidence of the rarefaction wave ${\psi }_4$ on the slipstream.

Interest in methods for rapid estimation of the parameters of the shock-wave structure of planar flows with Mach reflection, including asymmetric ones [15,16,17], has noticeably increased recently [15,16,17,18]. This fact can be associated with the development of aviation and rocket technology, flying at high supersonic speeds in the atmosphere.

The presence of strong shock waves in a supersonic flow of a reactive gas (fuel–air gaseous mixture) can initiate chemical reactions and detonation effects. In this regard, it is necessary to generalize the approximate analytical models [11,12,13], which have proven themselves well for the flows of a perfect non-reacting gas, to the case of a change in the chemical composition and pulsed energy release on emerging shocks. In [19,20,21], the classical shock-wave relations are generalized to the case of a change in the chemical composition and the impulse energy efflux within the framework of Chapman–Jouguet stationary detonation theory. In particular, “detonation” polars have been constructed instead of classical shock polars [19,20,21]; generalized criteria for regular or Mach shock reflection transition were derived, and their shift compared to the flow without energy supply and change in chemical composition was analyzed.

The numerical and theoretical results achieved in [19,20,21,22,23] are not always directly applicable in practice. In particular, the important analytical relations obtained in [19] are based on the assumption of the presence of detonation effects on both the incident ($j_1$) and main ($j_3$) shocks. At the same time, as it was shown in [5,6], it is a noticeably higher temperature of the gas mixture downstream from the main shock that initiates detonation (see also Figure 2). At the same time, neither a change in the chemical composition nor the energy release is observed on the incident shock or the reflected one. It is taken into account in study [20], where, however, to describe the flow on the main shock, the average (between the initial mixture and its combustion products) adiabatic index was introduced, but the mechanism of its calculation was not explained. The calculation of the chemical kinetics of initiated reactions [22,23] also does not assist the development of an approximate analytical technique for fast calculations.

This study generalizes the approximate analytical flow model with Mach reflection [11,12,13] to the case of a change in the chemical composition of the gas mixture and a pulsed energy supply at the main shock. On the example of a highly overexpanded jet of a methane–air mixture, we present an algorithm for calculating the parameters of the shock wave structure. The primary results of its application are shown, and a comparison is made with similar data for a highly overexpanded jet of non-reacting gas.

2. Model and Methods

2.1. Calculation of Parameters in Vicinity of the Triple Point

When a sufficiently strong oblique shock wave $j_1$ reflects in a flow with the Mach number $M$, a triple configuration of the second type (TC-2) according to the classification [24,25] forms at the triple point T. The parameters of the shocks $j_1-j_3$ are related by the following conditions of equality of static pressures and co-directionality of flows across a slipstream $\tau$:

$$J_1{J}_2=J_3,$$

(1)

$${\theta }_1+{\theta }_2={\theta }_3.$$

(2)

Here, $J_1=p_1/p$, $J_2={p_2/p}_1$, and $J_3=p_3/p$ are the strengths of incident shock ($j_1$), the reflected one ($j_2$), and the main one (Mach stem $j_3$), i.e., the ratios of static pressures on their sides. Additionally, $p$ is the static pressure in the incoming stream, $p_1-p_3$ are the pressures behind the shocks $j_1-j_3$, and ${\theta }_1-{\theta }_3$ are the flow deflection angles on those shocks in the vicinity of the triple point related to corresponding shock strengths:

$$\left|{\theta }_1\right|=\arctan \left[\sqrt{\frac{\left(1+\epsilon \right)M^2-J_1-\epsilon }{J_1+\epsilon }}\cdot \frac{\left(1-\epsilon \right)\left(J_1-1\right)}{\left(1+\epsilon \right)M^2-\left(1-\epsilon \right)\left(J_1-1\right)}\right]$$

(3)

$$\left|{\theta }_2\right|=\arctan \left[\sqrt{\frac{\left(1+\epsilon \right)M_1^2-J_2-\epsilon }{J_2+\epsilon }}\cdot \frac{\left(1-\epsilon \right)\left(J_2-1\right)}{\left(1+\epsilon \right)M_1^2-\left(1-\epsilon \right)\left(J_2-1\right)}\right]$$

(4)

on the incident shock and the reflected one [25], and by the relation [19]

$$\left|{\theta }_3\right|=\arctan \left[\frac{\left(J_3-1\right)\sqrt{F-1}}{\gamma M^2-\left(J_3-1\right)}\right],F=\frac{2\gamma M^2\left[\left(\gamma -{\gamma }_3\right)+\left(\gamma -1\right)\left(\left(J_3-1\right)-\left({\gamma }_3-1\right)\phi \right)\right]}{\left(\gamma -1\right)\left(J_3-1\right)\left[\left({\gamma }_3+1\right)\left(J_3-1\right)+2{\gamma }_3\right]}$$

(5)

for the Mach stem. Here, $\epsilon =\left(\gamma -1\right)/\left(\gamma +1\right)$, $\gamma$ is the adiabatic index of the gas mixture in the incoming stream, and ${\gamma }_3$ is the adiabatic index of the combustion products behind the Mach shock. Dimensionless quantity

$$\phi =\frac{\lambda }{\left(p/\rho \right)}=\frac{\gamma \lambda }{\left(\gamma -1\right)c_{p}T}$$

(6)

characterizes the pulse energy supply; the values of density $\rho$, temperature $T$, isobaric specific heat $c_p$, and released heat $\lambda$ (specific heat of combustion) are the parameters of the incoming stream.

The Mach numbers behind the incident and reflected shocks obey the relations

$$M_1=\sqrt{\frac{\left(J_1+\epsilon \right)M^2-\left(1-\epsilon \right)\left(J_1^2-1\right)}{J_1\left(1+\epsilon J_1\right)}},$$

(7)

$$M_{2T}=\sqrt{\frac{\left(J_2+\epsilon \right)M_1^2-\left(1-\epsilon \right)\left(J_2^2-1\right)}{J_2\left(1+\epsilon J_2\right)}}.$$

(8)

Classical relations based on the Chapman–Jouguet theory of stationary detonation wave are,

$$\rho u_n={\rho }_3{u}_{3n}, p+\rho u_n^2=p_3+{\rho }_3{u}_{3n}^2, u_{\tau }=u_{3\tau }, \frac{u_n^2+u_{\tau }^2}{2}+\frac{\gamma }{\gamma -1}\frac{p}{\rho }+\lambda =\frac{u_{3n}^2+u_{3\tau }^2}{2}+\frac{{\gamma }_3}{{\gamma }_3-1}\frac{p_3}{{\rho }_3},$$

where the indices “$n$” and “$\tau$” mean the velocity components normal and tangential to shock surface, and we determine the Mach number $M_{3T}$ in the vicinity of the triple point:

$$M_{3T}=\sqrt{\frac{\gamma M^2\left(E_3^2{\sin }^2{\omega }_3+{\cos }^2{\omega }_3\right)}{{\gamma }_3{E}_3{J}_3}}.$$

(9)

Here

$${\omega }_3=\arcsin \sqrt{\frac{J_3-1}{\gamma \left(1-E_3\right)M^2}}$$

is shock slope angle to flow direction upstream it, and

$$E_3=1-\frac{2\left[J_3-\left({\gamma }_3-1\right)/\left(\gamma -1\right)-\left({\gamma }_3-1\right)\phi \right]}{\left({\gamma }_3-1\right)+\left({\gamma }_3+1\right)J_3}$$

is the inverse ratio of the gas densities on the sides of the shock. At $\phi =0$ and ${\gamma }_3=\gamma$, it reduces to the ordinary Rankine–Hugoniot adiabat.

As shown in [20], a sufficiently large, pulsed energy supply $\phi$ displaces the “detonation” polar III (Figure 3), which obeys the relation (5), inside the shock polar I, corresponding to Equation (3). As it is seen from Figure 3, the strength of the steady Mach stem must belong to the interval

$$J_{3\min }\le J_3\le J_{3\max }.$$

(10)

Here, the value $J_3=J_{3\max }$ corresponds to a direct shock wave with pulsed energy supply, and the triple configuration with $J_3=J_{3\max }$ and $J_1=J_N$ (Figure 3) corresponds to an analogue of the von Neumann criterion [25,26] of shock reflection transition. Shock strengths $J_{3\min }$ and $J_{3\max }$ are obey the relations [23]

$$J_{3\min ,\max }=\frac{\gamma +1}{{\gamma }_3+1}\cdot \frac{J_m\left(M\right)+1}{2}\mp \frac{\sqrt{{\gamma }_{}^2{M}^4+{\gamma }_3^2-2\gamma M^2\left[\left({\gamma }_3^2-1\right)\phi +\left({\gamma }_3^2-\gamma \right)/\left(\gamma -1\right)\right]}}{{\gamma }_3+1},$$

(11)

$$J_m\left(M\right)=\left(1+\epsilon \right)M^2-\epsilon .$$

The only solution of the system (1-2) in the considered range of sufficiently large Mach numbers exists in the range of incident shock strengths

$$J_{1N}\le J_1\le J_{1P}.$$

(12)

Here, the value $J_1=J_{1P}$ corresponds to the limiting position of polar IIb of the reflected shock (its contact with polar III, see Figure 3). The flow behind the reflected shock $j_2$ is always supersonic ($M_{2T}>1$), while the flow behind the main shock $j_3$ is subsonic ($M_{3T}<1$).

As a rule [7,10,11,12,13,15,16,17,18] for an approximate analytical description of the flow in region 3 (Figure 1a,b), a model of a quasi-one-dimensional flow with some initial Mach number $M_{30}$ directly downstream from the main shock is used. The value $M_{30}$ can be determined by Formula (9), then $M_{30}=M_{3T}$, or a similar relation for $J_3=J_{3\max }$, which corresponds to the flow at point N (Figure 1) behind the direct shock (then $M_{30}=M_{3N}$), or the half-sum of these values. Then

$$M_{30}=\left(M_{3T}+M_{3N}\right)/2,$$

(13)

Corresponding to the approach adopted in [13,18]. The value of $M_{30}$ also can be based on more complex methods for flow parameters averaging [7,27]. Furthermore, we use approximation (13) to determine the initial Mach number of the flow along “virtual nozzle” 3.

Thus, the application of the proposed model makes it possible to establish the initial values of the flow angle (${\theta }_0={\theta }_3$) and Mach numbers ($M_{20}=M_{2T}$ and $M_{30}$) in regions 2 and 3 behind the reflected and main shocks.

It is necessary to add that energy supply should not be taken into account immediately at surface of the Mach stem, following the more advanced classical Zel’dovich–von Neumann–Döring (ZND) theory. Thus, solving the problem (1–2) at the triple point, we should use a relation similar to (3, 4) instead of (5) for the flow deflection angle at the shock $j_3$ and a relation similar to (7, 8) instead of (9) for Mach number after it. In this approach, “detonation” polar III in Figure 3 coincides with ordinary shock polar I, and the problem in the triple point is to be solved as an ordinary problem of Mach reflection in the non-reacting gas.

After it, according to ZND theory, a thin layer of chemical reactions initiated by high temperature after the Mach stem $j_3$ follows just after that shock. The static pressure drop after the Mach stem accompanies that chemical reaction, forming a so-called “von Neumann peak” or a “himpik” (as it is called in Soviet or post-Soviet literature). Vertical distance between “strong” upper branches of the polars I and III in Figure 3 characterizes that pressure drop from the “Hugoniot” parameters that corresponds to Mach reflection in non-reacting gas to the parameters that approximately correspond to Chapman–Jouguet theory applied in this study. Modernizing the proposed approximate analytical model, it is possible to take into account finite-rate kinetics downstream from the Mach stem. However, as the chemical reaction zone is very thin (about 10⁻⁴ m in typical considering conditions [28]), comparing with typical lengths in this study (such as channel width and Mach stem height), it seems admissible to compress that layer into alone surface of the shock $j_3$ and apply the Chapman–Jouguet theory to that shock immediately.

2.2. Flow in a Rarefaction Wave behind a Reflected Shock

Relations for flow direction angle $\theta$ of flow in the Prandtl–Meyer wave two with rectilinear acoustical characteristics of the first family

$$\theta ={\theta }_0+\nu \left(M_{20}\right)-\nu \left(M_2\right);$$

here, $\nu \left(M\right)=1/\sqrt{\epsilon }\arctan \sqrt{\epsilon \left(M^2-1\right)}-\arctan \sqrt{M^2-1}$ is the Prandtl–Meyer function, together with the conditions,

$$y/y_T=q\left(M_{30}\right)/q\left(M_3\right) \mathrm{and}$$

$q\left(M\right)=M{\left[1+{\epsilon }_3\left(M_{}^2-1\right)\right]}^{-1/2{\epsilon }_3}$, ${\epsilon }_3=\left({\gamma }_3-1\right)/\left({\gamma }_3+1\right)$ is the isentropic flow rate function, and the equality of pressures across the slipstream $\tau$ are written in the form

$$\pi \left(M_2\right)/\pi \left(M_{20}\right)={\pi }_3\left(M_3\right)/{\pi }_3\left(M_{30}\right);$$

here, $\pi \left(M\right)={\left[1+\left(\gamma -1\right)M^2/2\right]}^{-\gamma /\left(\gamma -1\right)}$ and ${\pi }_3\left(M\right)={\left[1+\left({\gamma }_3-1\right)M^2/2\right]}^{-{\gamma }_3/\left({\gamma }_3-1\right)}$ are the isentropic pressure functions written for gases above and below the slipstream, leading to the following equations for the shape $y\left(x\right)$ of the slipstream $\tau$ and the change in flow properties on its sides:

$$\begin{gathered} \frac{dy}{dx}=\tan \theta , \frac{d\theta }{dx}=-\frac{{\gamma }_3{M}_3^2\sqrt{M_2^2-1}\tan \theta }{\gamma M_2^2\left(M_3^2-1\right)y}, \frac{d{M}_2}{dx}=\frac{{\gamma }_3{M}_3^2\left[1+\epsilon \left(M_2^2-1\right)\right]\tan \theta }{\left(1+\epsilon \right)M_2\left(M_3^2-1\right)y}, \mathrm{and} \\ \frac{d{M}_3}{dx}=\frac{M_3\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]\tan \theta }{\left(1-{\epsilon }_3\right)\left(M_3^2-1\right)y}. \end{gathered}$$

(14)

Considering the Mach number $M_3$ below the slipstream as an independent variable, we use Equation (14) obtain the forms

$$\begin{gathered} \frac{dx}{d{M}_3}=\frac{\left(1-{\epsilon }_3\right)\left(M_3^2-1\right)y}{M_3\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]\tan \theta }, \frac{dy}{d{M}_3}=\frac{\left(1-{\epsilon }_3\right)\left(M_3^2-1\right)y}{M_3\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]}, \\ \frac{d\theta }{d{M}_3}=-\frac{\left(1+{\epsilon }_3\right)M_3\sqrt{M_2^2-1}}{\gamma M_2^2\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]}, \mathrm{and} \frac{d{M}_2}{d{M}_3}=\frac{\left(1+{\epsilon }_3\right)M_3\left[1+\epsilon \left(M_2^2-1\right)\right]}{\left(1+\epsilon \right)M_2\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]}. \end{gathered}$$

(15)

With the same chemical composition of flows on both sides of slipstream $\tau$, relations (14) and (15) take the simplified form known from [13,14]. Equations (14) or (15) are integrated up to the point D of intersection of the slipstream $\tau$ with the boundary characteristic BD of the second family, which falls from the exit point of the reflected shock $j_2$ (TB) on the jet boundary (Figure 1b). The flow below the slipstream must remain subsonic ($M_3<1$); the achievement of critical sound velocity indicates that the proposed triple point height is significantly underestimated.

Calculations show that the angle of turn of a slipstream in the TD section is usually small in its absolute value, but it can be 1.5–2 times greater than the initial angle ${\theta }_0$ of the slipstream’s inclination at the triple point. Therefore, ignoring the curvature of the slipstream [7] leads to a significant underestimation of the obtained height $y_T$ of the main shock.

An analysis of the differential characteristics of the flow in the Prandtl–Meyer wave, similar to that carried out in [13,14] at $\gamma ={\gamma }_3$, leads to the following equations for determining the shape $y_F\left(x_F\right)$ of the curvilinear boundary acoustic characteristic BD of the second family and flow parameters along the direction $\zeta$ of its incidence:

$$\begin{gathered} \frac{d{x}_F}{d\zeta }=\cos \left(\theta -{\mu }_2\right), \frac{d{y}_F}{d\varsigma }=\sin \left(\theta -{\mu }_2\right), \frac{d\theta }{d\zeta }=-\frac{2{\gamma }_3{M}_3^2\left(M_2^2-1\right)\sin \theta }{\gamma M_2^3\left(M_3^2-1\right)y_{F_1}}\cdot {\left[\frac{y_F-y_{F_1}}{\Delta y}+1\right]}^{-1}, \\ \frac{d{M}_2}{d\zeta }=\frac{2{\gamma }_3{M}_3^2\left[1+\epsilon \left(M_2^2-1\right)\right]\sqrt{M_2^2-1}\sin \theta }{\left(1+\epsilon \right)M_2^2\left(M_3^2-1\right)y_{F_1}}\cdot {\left[\frac{y_F-y_{F_1}}{\Delta y}+1\right]}^{-1}, \mathrm{and} \\ \frac{d{M}_3}{d\zeta }=\frac{2M_3^{}\left[1+{\epsilon }_3\left(M_3^2-1\right)\right]\sqrt{M_2^2-1}\sin \theta }{\left(1-{\epsilon }_3\right)M_2^{}\left(M_3^2-1\right)y_{F_1}}\cdot {\left[\frac{y_F-y_{F_1}}{\Delta y}+1\right]}^{-1}. \end{gathered}$$

(16)

Here, ${\mu }_2=\arcsin \left(1/M_2\right)$ is the Mach angle; $M_2$ is flow Mach number at an arbitrary point on the characteristic BD; and $M_3$ is the Mach number of the flow on the other side of the slipstream $\tau$ at the point $F_1$ of its intersection with the corresponding rectilinear characteristic $F{F}_1$. The coordinate $y_{F_1}$ obeys the formula

$$y_{F_1}=y_{T}q\left(M_{30}\right)/q\left(M_3\right),$$

and the distance $\Delta y$ from this point to the envelope of the family of rectilinear characteristics can be calculated by such a way:

$$\Delta y=\frac{\left(1+\epsilon \right)\sqrt{M_2^2-1}\left(M_3^2-1\right)y_{F_1}\sin \left(\theta +{\mu }_2\right)}{{\gamma }_3{M}_2{M}_3^2\sin \theta }.$$

The integration of the system (16) starts from point B of the intersection of the reflected shock with the jet boundary and ends at point D of the intersection of the characteristic BD with the slipstream. If we take the distance L between the curved characteristic BD and the slipstream as an independent variable, then Equation (16) obtains the forms

$$\begin{gathered} \frac{d{x}_F}{dL}=\frac{d{x}_F/d\varsigma }{dL/d\varsigma }, \frac{d{y}_F}{dL}=\frac{d{y}_F/d\zeta }{dL/d\zeta }, \frac{d\theta }{dL}=\frac{d\theta /d\zeta }{dL/d\zeta }, \frac{d{M}_2}{dL}=\frac{d{M}_2/d\zeta }{dL/d\zeta }, \mathrm{and} \\ \frac{d{M}_3}{dL}=\frac{d{M}_3/d\zeta }{dL/d\zeta }. \end{gathered}$$

(17)

The value

$$\frac{dL}{d\zeta }=\sin \left({\mu }_2-\theta \right)+\frac{2\sqrt{M_2^2-1}\sin \theta }{M_2}\cdot {\left[\frac{y_F-y_{F_1}}{\Delta y}+1\right]}^{-1}$$

characterizes the decrease in the distance L as the curvilinear characteristic BD approaches the slipstream. Equation (17) should be integrated starting from the value L = |B₁B|, which is the initial distance between the characteristic BD and the slipstream, and the flow parameters correspond to point B behind the reflected shock to the value L = 0 at point D of intersection of the slipstream and the characteristic BD.

When taking into account finite-rate kinetics (for example, in classical ZND theory), corresponding elements should be introduced into model of flow in zone III (especially just after the Mach stem). Though the zone of chemical reactions is very thin, its introduction into the proposed model can be a subject of further studies.

2.3. Shape and Parameters of the Reflected Curved Shock

It has been repeatedly shown [29,30] that the reflection coefficient of rarefaction or compression disturbances that overtake the previous oblique shock is very small (if the flow behind the shock is not transonic, which is not the case now). In addition, the intensity of the rarefaction wave acting on the reflected shock $j_2$ in the triangle TBB₁ (Figure 1b) is small in itself. Therefore, we apply the method first proposed in [31], and the shape of the reflected shock is determined based on the condition of flow direction conjugation (i.e., the flow angles behind each point of the reflected shock $j_2$ and at the corresponding points of the Prandtl–Meyer expansion wave 2 must coincide). The analysis of the differential characteristics of the flowfield in wave 2 [13,14], together with the conditions on the oblique shock, determine the shape of the reflected shock in polar coordinates ($r$,$\xi$):

$$\frac{dr}{d\xi }=r\cot \left({\omega }_2-{\theta }_1-\xi \right)$$

(18)

$$\mathrm{and} \frac{d{\theta }_2}{d\xi }=\frac{K_s{M}_{2}r\sin \left({\mu }_2+{\theta }_2-{\omega }_2\right)}{\sin \left({\omega }_2-{\theta }_1-\xi \right)}.$$

(19)

Here, $r$ is the distance from the nozzle edge to the considered point on the shock, $\xi$ is the polar angle measured from the horizontal direction (see Figure 1b), and ${\mu }_2=\arcsin \left(1/M_2\right)$ is the Mach angle. Flow deflection angle ${\theta }_2$ on the reflected shock depends on shock strength $J_2$:

$$\tan {\theta }_2=\sqrt{\frac{\left(1+\epsilon \right)M_1^2-J_2-\epsilon }{J_2+\epsilon }}\cdot \frac{\left(1-\epsilon \right)\left(J_2-1\right)}{\left(1+\epsilon \right)M_1^2-\left(1-\epsilon \right)\left(J_2-1\right)},$$

and the shock slope angle ${\omega }_2$ obeys the dependence

$$J_2=\left(1+\epsilon \right)M_1^2{\sin }^2{\omega }_2-\epsilon .$$

The Mach number $M_2$ at the corresponding point of the rarefaction wave behind the shock is determined by the relation

$$\nu \left(M_2\right)=\nu \left(M_{20}\right)+{\theta }_0-{\theta }_2.$$

Here, $\nu \left(M\right)$ is the Prandtl–Meyer function. The corresponding Mach number $M_3$ on the other side of the slipstream obeys the relation

$$\pi \left(M_3\right)=\pi \left(M_2\right)\cdot {\pi }_3\left(M_{30}\right)/\pi \left(M_{20}\right).$$

The curvature of the streamline $K_S$, which is present in (19) behind an arbitrary point on shock surface, is defined as follows:

$$K_S=\frac{{\gamma }_3{M}_3^2\sqrt{M_2^2-1}\sin \theta }{\gamma M_2^2\left(M_3^2-1\right)y_A}\cdot \frac{y_A-y_{A_1}}{y-y_{A_1}}.$$

Here, $y=H+r\sin \xi$ is the ordinate of a given point behind the shock, $\theta ={\theta }_0+\nu \left(M_{20}\right)-\nu \left(M_2\right)$ is the flow angle at this point, $y_A=y_{T}q\left(M_{30}\right)/q\left(M_3\right)$, $y_{A_1}=y_A\left[1-\left(1+\epsilon \right)\sqrt{M_2^2-1}\left(M_3^2-1\right)\sin \left({\mu }_2+\theta \right)/\left({\gamma }_3{M}_2{M}_3^2\sin \theta \right)\right]$, and $H$ is the half-width of the nozzle exit section.

Equations (18) and (19) are to be integrated from the value $\xi =-{\omega }_1$ at the triple point T to the value $\xi =-{\theta }_1$ at point B of the shock $j_2$ exit to the jet boundary with subsequent reflection of the rarefaction wave ${\psi }_4$.

Our calculations, according to (18–19), demonstrate that the geometrical curvature of the reflected shock $j_2$ is usually almost negligibly small (its slope angle diminishes on 0.1–0.3°). When applying any finite-rate kinetics model (as in classical ZND approach), a small segment of relatively large curvature of the reflected shock $j_2$ (TB) is present in the vicinity of the triple point T. At this segment, reflected shock strength quickly diminishes from “Hugoniot” parameters to “Chapman-Jouguet” ones under the influence of pressure drop in von Neumann peak after the Mach stem. In our model, the initial parameters of the reflected shock are directly predicted by Chapman–Jouguet theory, and the influence of the finite rate of the chemical reaction is subject to further studies and comparisons.

2.4. Incidence of a Rarefaction Wave on a Slipstream

Many modern approximate models of flows with Mach reflection [15,16,17] actually replace a fast analytical assessment of the interaction of a rarefaction wave ${\psi }_4$ (Figure 1a,b) with a slipstream $\tau$ with a calculation by the method of characteristics. However, on the other hand, the reduction in the DCE region of this interaction to a single point D, ignoring the finite length of the interval DC of the slipstream turn in the horizontal direction [7,10], also leads to significant errors in determining the parameters of the shock-wave structure.

By analogy with the model [11,12,13], the analysis of the interaction of a wave with a slipstream is carried out as follows. It is assumed that the Prandtl–Meyer flow is realized in the wave, and the Mach number $M_{2BD}$ of the flow in front of this wave corresponds to the average slope of the characteristic BD in region 2:

$$\nu \left(M_{2BD}\right)+\mu \left(M_{2BD}\right)=\left[\nu \left(M_{2B}\right)+\mu \left(M_{2B}\right)+\nu \left(M_{2D}\right)+\mu \left(M_{2D}\right)\right]/2$$

(20)

The initial Mach number $M_{3D}$ in flow region 3 and the width $y_D$ of this “virtual nozzle” are determined from the previously obtained calculations of the shape of the slipstream in the section TD, the boundary characteristic BD, and the variation in the flow parameters along them.

An arbitrary Mach number $M_{3H}\in \left[M_{3D}; 1\right]$ on the lower side of the slipstream corresponds to the Mach number $M_{2H}$ on its upper side according to the relation

$$\pi \left(M_{2H}\right)/{\pi }_3\left(M_{3H}\right)=\pi \left(M_{2BD}\right)/{\pi }_3\left(M_{3D}\right),$$

(21)

and the width $y_H$ of area 3 is such that

$$y_H/y_D=q\left(M_{3D}\right)/q\left(M_{3H}\right).$$

(22)

In this case, the abscissa $x_H$ of the point H is approximately determined from the condition of straightness of the voluntary incident characteristic BH₁H:

$$\left(y_H-y_B\right)/\left(x_H-x_B\right)=\tan \left[\zeta +\nu \left(M_{2BD}\right)+\mu \left(M_{2BD}\right)-\nu \left(M_{2H_1}\right)-\mu \left(M_{2H_1}\right)\right]$$

(23)

Here, $\tan \zeta =\left(y_D-y_B\right)/\left(x_D-x_B\right)$. Slope $y_H^{’}\left(x_H\right)$ of the slipstream at this point obeys the condition

$$y_H^{’}\left(x_H\right)=\tan \left[{\theta }_D+2\nu \left(M_{2H_1}\right)-\nu \left(M_{2BD}\right)-\nu \left(M_{2H}\right)\right],$$

(24)

which takes into account the flow deflection in the waves incident on the slipstream $\tau$ $({\psi }_4)$ and $({\psi }_5)$ reflected from it. Here, $M_{2H_1}$ is the Mach number on the BH₁ characteristic, and ${\theta }_D$ is the previously determined slope of the slipstream at point D.

The integration of the system (21-24) determines the shape $y_H\left(x_H\right)$ of the slipstream $\tau$. It is carried out until Mach number value below the slipstream reaches unity.

If the rarefaction wave ${\psi }_4$ experiences refraction on a reflected shock $j_2$ (see, for example, Figure 1a) before falling onto a slipstream, we apply the analytical solutions for the oblique shock–rarefaction fan interaction, obtained in [32,33].

2.5. Algorithm of Application of Our Generalized Approximate Analytical Model

To quickly estimate the parameters of the shock-wave structure of a supersonic flow with a Mach reflection, a pulsed energy supply, and a change in the chemical composition on the main shock, first of all, the flow parameters of the initial gas mixture (values $M$ and $\gamma$), as well as the adiabatic index ${\gamma }_3$ of the combustion products, the dimensionless pulsed energy supply $\phi$, and the strength $J_1$ of the incident shock, which satisfies inequality (12), should be given. Further, the flow parameters are determined as follows:

• Using the relations (1–7), we calculate the flow parameters in the vicinity of the triple point. Relations (2), (8), (9), and (13) determine the Mach numbers $M_{2T}=M_{20}$ and $M_{3T}$ on both sides of the slipstream, the initial Mach number $M_{30}$ of the flow along “virtual nozzle” 3, and the angle ${\theta }_0={\theta }_3$ of the slipstream inclination in the vicinity of the triple point.
• The value of the height $y_T$ of the triple point in the first approximation should be set.
• Relations (18–19) establish the shape of a slightly curved (upwards–convex) reflected shock wave $j_2$ (TB), as well as the coordinates of point B and the flow parameters behind the shock at this point.
• Simultaneously, the shape of the boundary characteristic BD (16–17), the slipstream $\tau$ (15), and the flow parameters along them should be determined. The integration of Equations (15)–(17) ends at point D of intersection of the boundary characteristic BD with a slipstream. The flow along “virtual nozzle” 3 can reach the critical velocity ($M_3=1$) up to this point. It indicates a significant decrease in the value $y_T$ in the accepted approximation.
• The problem (20–25) of the reversal of the slipstream $\tau$ under the influence of the rarefaction wave ${\psi }_4$ is being solved. If it is found that the critical flow velocity in area 3 is reached before slipstream $\tau$ turns to the horizontal direction ($M_{3H}=1$ at $y_H^{’}\left(x_H\right)<0$), the proposed height $y_T$ is too small; otherwise, it is excessive.
• According to the results of the iteration, we refine the value $y_T$ of Mach stem height, and the calculations should be made in the next approximation, returning to point 2.

3. Results of the Application of the Approximate Analytical Model

As an example of the application of the presented model, the outflow of a uniform (in the exit section) planar jet of a fuel–air gas mixture is calculated. The adiabatic indices of upstream flow and the gas mixture behind the main shock were taken equal, correspondingly, to $\gamma =1.396$ and ${\gamma }_3=1.290$. They approximately corresponded to the stoichiometric methane–air mixture (the mass fractions were 4.856% CH₄, 74.212% N₂, 19.980% O₂, 0.951% of other impurities) and the products of its complete combustion, including carbon dioxide gas and water vapor. The specific heat of combustion of the fuel was assumed to be 55.266 MJ/kg, which corresponded to the value of $\lambda =2.684$ MJ/kg in terms of the entire gas mixture and the dimensionless value of the pulsed energy supply (6) $\phi =30.045$ at T = 300 K. This value of the dimensionless specific energy supply was much greater than that considered in [19].

Due to the large sudden energy supply, the solution of system (1-2) with restrictions (1–7), which describes the stationary Mach reflection, existed only at $M>5.436$. According to condition (11), which limited the range of the possible strength of the main shock at $M=5.436$, the detonation polar III (Figure 3) is compressed to the point

$$J_{3\min }=J_{3\max }=\frac{\gamma +1}{{\gamma }_3+1}\cdot \frac{J_m\left(M\right)+1}{2}=\frac{\gamma +1}{{\gamma }_3+1}\cdot \frac{\left(1+\epsilon \right)M^2+\left(1-\epsilon \right)}{2}=\frac{\gamma M^2+1}{{\gamma }_3+1},$$

and it disappears at lower Mach numbers. In such a case, the solution describing the triple configuration of the Mach reflection should be sought, admitting the propagation of the main shock (detonation) wave $j_3$ upstream from the incoming flow.

Examples of calculation of the dimensionless (referred to the half-width of the exit section) height of the triple point $y_T$, depending on the incident shock strength (provided that the considered strongly overexpanded jet flows out without separation), is shown in Figure 4. Curves I–VI correspond to Mach numbers M = 6, 6.5, 7, 8, 9, and 10 in the entire range (12) of theoretically possible incident shock strengths.

For comparison, curves I’–VI’ show the values of $y_T\left(J_1\right)$ calculated for the same flow parameters but without the energy release ($\phi =0$) and changes in chemical composition of the gas mixture (${\gamma }_3=\gamma =1.396$), i.e., according to the algorithm [11,12,13], a generalization of which is presented in this paper. From a comparison of curves I–VI and I’–VI’, it is obvious that a significant pulsed energy supply significantly shifts the conditions of transition to the Mach reflection (from points 1′–6′, which correspond to “classical” von Neumann criterion, to points 1–6). In addition, the energy release leads to a significant increase in the size of the main (Mach) shock and the width of the subsonic flow region behind it.

The next example of the application of the proposed approximate analytical model relates to Figure 2. The planar flow of the stoichiometric methane–air gas mixture in the narrowing channel between two wedges was studied numerically using the ANSYS Fluent 2020 R2 CFD software package (academic version). It was admitted that the wedge angle ${\theta }_1={24}^{\circ }$, which is also flow deflection angle on the incident shock; the distance between two front edges is equal to 200 mm, and the minimal width of the channel is equal to 100 mm. The inflow had normal atmospheric parameters (static pressure and temperature), its dynamic viscosity obeyed the Sutherland law, and ordinary $\kappa -\omega$ SST model of turbulence was applied. In one approach, the finite-rate model was applied to the calculations of the kinetics of chemical reactions initiated by high static temperature downstream from the Mach stem. In another approach, chemical reactions were artificially “switched off”, and no variations in chemical composition of the mixture occurred. The unstructured numerical grid contained about 2·10⁶ cells. Calculated dimensionless (referred to the width of the entrance section) heights $y_T$ were approximately equal to 0.377 and 0.207 in flows with present and artificially absent chemical reactions, correspondingly. The first calculated example (at presence of chemical reactions) is demonstrated in Figure 2a,b. It can be also seen in Figure 2b that chemical reactions are initiated only downstream from the main shock but not after the incident shock and the reflected one (generally, products of combustion, such as carbon dioxide, are present only in flow region after the Mach stem).

Our approximate approach to this flow, which includes also the model of oblique shock–expansion fan interaction derived in [32,33], gives us $y_T=0.412$ in the first case and $y_T=0.215$ in the second one. Thus, the comparison between numerical and analytical results seems quite satisfactory.

At the next stage, we provided a comparison of our analytical model with results of numerical experiment in [34]. At first, the flow of a lean hydrogen air mixture (0.4 H₂ + 4.772 Air) was studied from the upstream Mach numbers M = 3.0 and M = 3.15. Such a problem had no steady numerical solution with the Mach reflection in [34] (transition to regular reflection occurred, or auto-oscillations of Mach stem with detonation downstream started). We presumed the full combustion of hydrogen (its specific heat of combustion could be estimated as 55.266 MJ/kg, thus $\lambda =0.821$ MJ/kg was the energy supply per unit mass of gas mixture. It caused the energy efflux to exceed the critical value, which was equal to 0.538 MJ/kg at M = 3.15 and 0.477 MJ/kg at M = 3.0. At those values, the detonation polar for the Mach stem compressed into a single point and disappeared, thus there was no solution for stationary detonation (at least at $\gamma =1.403$ and ${\gamma }_3=1.322$, which were adiabatic indices of gas mixtures upstream from the shock and downstream from it). Thus, the coincidence between analytical and numerical results was also obtained, though it was a negative one. A more positive coincidence for that gas mixture can be obtained, for example, at M = 4.0 (stationary solutions exist there, and reasonable results can be provided).

At the last stage, the flow of a stoichiometric hydrogen–oxygen mixture (2 H₂ + O₂) was considered. We analyzed a supersonic jet flow of a gas mixture with chemical parameters and a Mach number similar to the number considered in [34]. At M = 4.0, the energy supply also exceeded the critical value, thus it is not a wonder that only unsteady flow regimes were obtained in [34] at ${\theta }_1={25}^{\circ }$.

At M = 7.0, the stationary detonation was computationally obtained in [34] only within a small range of flow deflection angles at the incident shock (between ${\theta }_1={20}^{\circ }$ and ${\theta }_1={23.5}^{\circ }-{24}^{\circ }$). The upper limit completely obeys our theory (maximal value $J_1=J_{1P}$ that corresponds to shock polar IIb in Figure 3 is to be obtained at ${\theta }_1={23.824}^{\circ }$). According both to the numerical experiment in [34] and to our model, the Mach stem occupies almost all of the cross-section of the nozzle or the narrowing channel at $J_1=J_{1P}$.

The lower limit (${\theta }_1={20}^{\circ }$) in [34] quite corresponds to the “classical” von Neumann criterion of shock reflection transition (${\theta }_1={19.716}^{\circ }$). Our theory diminishes the lower limit of Mach reflection to ${\theta }_1={17.104}^{\circ }$ (one can compare points 1 and 1′ in Figure 5). It is difficult to determine which criterion is more correct in practice because that range of flow deflection angles belongs to so-called dual solution zone (up to ${\theta }_1={30.245}^{\circ }$, according to classic theory for perfect gas, [26]), and we cannot answer whether the Mach reflection or the regular one is realized. At least, our computations demonstrate the regular reflection even at ${\theta }_1={22}^{\circ }$, though the Mach stem size at the Mach reflection of such an incident shock should be sufficiently large. Moreover, chemical reactions began at M = 7.0 and M = 8.0, even downstream from the incident shock and not only after the Mach one.

Dependencies $y_T\left({\theta }_1\right)$ of the Mach stem height on flow deflection angle on the oblique shock that falls onto the symmetry plane in supersonic jet flow are shown in Figure 5 for both “detonation” and “non-detonation” cases at M = 7.0 (curves I and I’) and M = 8.0 (curves II and II’).

4. Conclusions

The temperature of the flow behind the Mach stem, which is formed in high supersonic gas flows, is much larger than behind the incident shock and the reflected one. For this reason, detonation effects (impulse release and change in the chemical composition of the reactive gas mixture) are initiated, first of all, behind the main (Mach) shock. An analysis of the resulting triple configurations using the Chapman–Jouguet stationary detonation model at the main shock shows that a significant pulsed energy release leads to a shift in the criterion of shock reflection transition. The shock waves, which are regularly reflected according to the classical theory, can be reflected with formation of a triple point if sufficient energy release at the Mach stem is allowed.

This paper presents an approximate analytical model for fast calculation of the parameters of the shock-wave structure of the flow of a reactive gas mixture with a Mach reflection. The proposed model for the first time takes into account the change in the chemical composition and the pulsed energy supply at the main shock. The primary results obtained in the calculation of a supersonic jet flow of a stoichiometric composition of a fuel–air gas mixture show not only an earlier occurrence of Mach reflection compared to a similar flow without chemical reactions but also a significant increase in the geometric dimensions of the main shock (Mach stem).