High-Temperature Effects on Supersonic Flow around a Wedge

Abstract

The development of the flow pattern that arises during the interaction of a shock wave with a wedge is discussed. Mathematical modeling of the flow around the wedge is carried out with Euler equations. These equations describe the unsteady flow of an inviscid compressible fluid around a wedge in a two-dimensional domain. To take into account high-temperature effects on super- and hypersonic flows, the model developed takes into account equilibrium chemical reactions in the air, ionization, and dissociation processes. The initial parameters of the flow are set equal to the parameters of the flow behind the shock wave in accordance with the Rankine–Hugoniot relations. The solutions to the problem obtained with the model of ideal perfect gas and the model taking into account high-temperature effects in the air are compared. The influence of high-temperature effects on the distribution of flow quantities is discussed.

1. Introduction

The complexity of various technical problems associated with the design and development of hypersonic aircraft leads to the need for research in the field of aerodynamics and heat transfer using mathematical modeling. Hypersonic aircraft is characterized by flat aerodynamic shapes of streamlined surfaces [1]. As a propulsion system, it is supposed to use a ramjet engine with supersonic combustion integrated with the body [2]. A characteristic feature of aircraft with an air intake formed by the planes of the aircraft and elements of the power plant is the presence of extended structural elements in the form of a wedge with a small opening angle and a blunting of a small radius, subject to the most intense aerodynamic heating [3, 4].

The supersonic flow around a wedge is one of the well-studied problems in computational fluid dynamics (CFD) [5]. This interaction results in the formation of an impact configuration with two triple points [6]. The flow disturbance starts from the wedge tip when the shock wave touches it. Depending on the Mach number and the angle of the wedge, regular reflection (shock is attached to wedge) or non-regular reflection (a reflected shock wave has a curvilinear front) is generated. Depending on inlet flow conditions (angle of wedge and inlet Mach number), the main and reflected waves meet at a triple point. Between this point and the wedge, two waves merge into one wave, forming a Mach configuration [7, 8]. Downstream of the triple point, a tangential discontinuity is formed on which pressure and normal velocity remain continuous, while density and tangential velocity suffer a discontinuity.

CFD tools are applied to simulate regular and non-regular shock reflection in [9, 10]. In other studies, shock reflection from a line of symmetry is discussed when two wedges are placed in a supersonic flow [11]. This formulation of the problem allows to exclude boundary layer effects [12].

The inviscid compressible flow around a wedge is carried out in [13]. When the angle of flow turning, equal to the angle of inclination of the wedge, is less than the maximum, the problem has two solutions. In the solution with an oblique shock of lesser intensity (a “weak” shock), the uniform flow between the shock and the wedge is almost always supersonic [14]. The exception is a small neighborhood of the maximum angle of rotation. For a perfect gas, this neighborhood does not exceed fractions of a degree for all Mach numbers of the oncoming flow. After a shock of greater intensity (a “strong” shock), the flow of a perfect gas is always subsonic.

For supersonic and hypersonic flows of an inviscid gas, there is the possibility of the existence of two solutions (strong and weak). Depending on the angle of incidence and the input value of the Mach number, two typical configurations are formed: two-hop (regular) and three-hop (Mach), and in a certain range of parameters, both options. The possibility of hysteresis with a change in the angle of incidence of the shock has been shown in both physical and numerical experiments [12].

The experimental results of the flow around aircraft model in a hypersonic shock tube are presented in [15]. The results of the numerical simulation of the thermal state of a sharp wedge with a blunt edge in a high-speed airflow are presented in [16]. Accounting for viscosity and turbulence complicates the situation due to such effects as separation and reattachment of the boundary layer [17].

When a strong shock wave moves and interacts with a body, the temperature and pressure of the gas behind its front increase. The perfect gas model does not provide the required accuracy of the numerical solution, since the molecular weight and heat capacities are not constant. In high-temperature air, these quantities are functions of pressure and temperature. In this case, it is necessary to take into account the processes of dissociation and ionization taking place in a gas. In practice, various models are used that take into account high-temperature processes in gases, as well as analytical dependencies and interpolation of tabular values. From a computational point of view, the model proposed in [18] (Kraiko model) for air and taking into account the reactions between 13 components are interesting and successful. The main advantage of this model is that it takes into account the dissociation and ionization of air at high temperatures. In the temperature range up to 20,000 K and pressures from 0.001 to 1000 atm, the error of the model does not exceed ±2%, usually falling within the band of ±1%. Accounting for non-equilibrium chemical reactions is discussed in [19] based on a one-temperature model.

In this study, calculations are performed with the model of perfect gas and the model taking into account high-temperature effects. The results of numerical calculations obtained with different models are compared with each other and with theoretical values from the theory of oblique shock waves.

2. Formulation of problem

The governing equations are solved in the domain shown in Figure 1. The wedge angle is β. The length of the domain is 3.2, and its height is 2.2. The wedge is located at a distance of 0.2 from the origin. The shock is at the point x = 0 at time t = 0. There is no gas motion ahead of the shock wave (u = v = 0), and the pressure and temperature are p = 101325 Pa and T = 288.2 K. The flow velocity is fixed at 1000 and 3000 m/s. The flow quantities behind the shock wave are found from the Rankine-Hugoniot conditions. The gas does cross the centerline. Free outflow conditions are applied to the upper and right boundaries.

Figure 1.

Computational domain.

For this problem, there is an exact solution in the framework of the theory of oblique shock waves [1]. With the help of this problem, the capabilities of the scheme for reproducing shock waves are checked, which make it possible to assess the discrepancy between CFD calculation and the exact solution. As a parameter characterizing the proximity of the numerical solution to the analytical one, the Mach number is considered, which changes abruptly when passing through the shock wave.

3. Oblique shock wave

When a supersonic flow flows around a wedge, under certain restrictions on the half-angle of the wedge and the Mach number, an oblique shock occurs. When flowing around a cone, the shock front has a conical surface.

The angle of the wedge, β, is equal to the angle of flow turning at the shock. The angle between the shock front and the direction of the undisturbed flow is the shock slope angle, σ. The velocity of the undisturbed flow, v₁, is decomposed into the normal and tangential components to the shock surface, v_n1 and v_τ1. Therefore

$$v_{n1}=v_{1}sin\sigma ,v_{\tau 1}=v_{1}cos\sigma .$$

Flow quantities before and behind shock wave are interconnected by relationships following from the laws of conservation of mass, momentum, and energy:

conservation of mass

$$\begin{array}{l}{\rho }_1{v}_{n1}={\rho }_2{v}_{n2};\end{array}$$E1

conservation of momentum in the direction normal to the shock

$$\begin{array}{l}{\rho }_1{v}_{n1}^2+p_1={\rho }_2{v}_{n2}^2+p_2;\end{array}$$E2

conservation of momentum in the direction tangential to the shock (the pressure gradient in the direction tangential to the surface of the shock is equal to zero)

$$\begin{array}{l}{\rho }_1{v}_{n1}{v}_{\tau 1}={\rho }_2{v}_{n2}{v}_{\tau 2};\end{array}$$E3

conservation of energy

$$\begin{array}{l}{\rho }_1{v}_{n1}\left({\epsilon }_1+\frac{v_1^2}{2}+\frac{p_1}{{\rho }_1}\right)={\rho }_2{v}_{n2}\left({\epsilon }_2+\frac{v_2^2}{2}+\frac{p_2}{{\rho }_2}\right).\end{array}$$E4

Here, ε is specific internal energy. In relation (4), the terms in parentheses represent the sum of the specific internal and specific kinetic energy before and behind the shock. The change in this quantity is associated with the work performed on a given mass of gas by external forces, of which only surface pressure forces are taken into account. Taking into account the mass conservation condition (1), the requirement of equality of the tangential velocity components on the shock is v_τ1=v_τ2. Dividing both parts of Eq. (4) on the shock by the mass flux through the shock, ρ₁v_n1=ρ₂v_n2, the relation has the form

$$\begin{array}{l}{h}_1+\frac{v_1^2}{2}=h_2+\frac{v_2^2}{2},\end{array}$$

where $h=\epsilon +p/\rho$ is specific enthalpy.

The given conservation laws are valid for any gas model (perfect, real, dissociating, or ionized) when passing through an oblique shock wave since they express the general relations of the conservation laws without reference to any relations connecting thermodynamic variables to each other, and relations determining the form thermodynamic functions.

To close the conditions of dynamic compatibility at the shock, it is necessary to give specific dependencies that determine the specifics of the thermodynamic state of the gas. The enthalpy and molar mass are functions of pressure and temperature, h=h(p, T) and μ=μ(p, T). The equation of state of an ideal gas is taken, which has a composition corresponding to this state

$$\begin{array}{l}p=\frac{\rho R_{0}T}{\mu \left(p,T\right)},\end{array}$$

where R₀ is the universal gas constant. The thermodynamic quantities before the shock are known (h₁=c_p1T₁, μ₁=0.029 kg/mol for air), and the relations relating to the thermodynamic parameters behind the shock are introduced into the model from additional conditions describing the thermodynamic model of high-temperature air [3, 4].

To determine thermodynamic quantities of high-temperature flow behind a shock wave (v_τ1=v_τ2), the following equations are applied

$$h_2=h_2\left(p_2,T_2\right),{\mu }_2={\mu }_2\left(p_2,T_2\right),p_2=\frac{{\rho }_2{R}_0{T}_2}{{\mu }_2}.$$

For a perfect gas, simple transformations allow one to find the thermodynamic quantities behind the shock [5]

$$\begin{array}{l}\frac{p_2}{p_1}=\left(\frac{2\gamma }{\gamma +1}{\mathrm{M}}_1^2{sin}^2\sigma -\frac{\gamma -1}{\gamma +1}\right);\\ \frac{{\rho }_2}{{\rho }_1}={\left[\frac{2}{\left(\gamma +1\right){\mathrm{M}}_1^2{sin}^2\sigma }+\frac{\gamma -1}{\gamma +1}\right]}^{-1};\\ \frac{T_2}{T_1}=\frac{{\rho }_1}{{\rho }_2}\frac{p_2}{p_1};\\ \frac{{\mathrm{M}}_2}{{\mathrm{M}}_1}={\left\{\frac{T_1}{T_2}\left[{cos}^2\sigma +{sin}^2\sigma {\left(\frac{\gamma -1}{\gamma +1}+\frac{2}{\left(\gamma +1\right){\mathrm{M}}_1^2{sin}^2\sigma }\right)}^2\right]\right\}}^{1/2}.\end{array}$$

Here, σ is the shock angle.

Using the replacement ${\rho }_2=p_2{\mu }_2/\left(R_0{T}_2\right)$, the conditions of dynamic compatibility on the shock take the form

$$\begin{array}{l}{p}_2{\mu }_2\left(p_2,T_2\right)v_{n2}-C_1{R}_0{T}_2=0;\\ C_1{v}_{n2}+p_2-C_2=0;\\ h_2\left(p_2,T_2\right)+\left(\frac{v_{n2}^2+v_{\tau 2}^2}{2}\right)-C_3=0.\end{array}$$

Here, $C_1={\rho }_1{v}_{n1}$, $C_2={\rho }_1{v}_{n1}^2+p_1$, $C_3=h_1+v_1^2/2$, $v_{\tau 2}=v_{\tau 1}$.

To solve a nonlinear system of equations, Newton’s method is used, which consists of solving a sequence of linear problems

$$\begin{array}{l}{U}^{n+1}=U^n-J^{-1}F,\end{array}$$

where

$$\begin{array}{l}U=\left(\begin{array}{l}{p}_2\\ T_2\\ v_{n2}\end{array}\right),F=\left(\begin{array}{c}{p}_2{\mu }_2\left(p_2,T_2\right)v_{n2}-C_1{R}_0{T}_2\\ C_1{v}_{n2}+p_2-C_2\\ h_2\left(p_2,T_2\right)+0.5\left(v_{n2}^2+v_{\tau 2}^2\right)-C_3\end{array}\right).\end{array}$$

Jacobian has the form

$$\begin{array}{c}J=\frac{\partial F}{\partial U}=\left(\begin{array}{ccc}{v}_{n2}\left[{\mu }_2\left(p_2,T_2\right)+p_2\frac{\partial {\mu }_2}{\partial p_2}\right]& p_2{v}_{n2}\frac{\partial {\mu }_2}{\partial T_2}-C_1{R}_0& p_2{\mu }_2\\ 1& 0& C_1\\ \frac{\partial h_2}{\partial p_2}& \frac{\partial h_2}{\partial T_2}& v_{n2}\end{array}\right).\end{array}$$

The inverse matrix is written as

$$\begin{array}{l}{J}^{-1}=\frac{1}{\mid J\mid }\left(\begin{array}{lll}{J}_{11}& J_{12}& J_{13}\\ J_{21}& J_{22}& J_{23}\\ J_{31}& J_{32}& J_{33}\end{array}\right)\end{array}.$$

Here

$$\begin{array}{l}{J}_{11}=-C_1\frac{\partial h_2}{\partial T_2},\\ J_{12}=C_1{R}_0{v}_{n2}+p_2{\mu }_2\frac{\partial h_2}{\partial T_2}-p_2{v}_{n2}^2\frac{\partial {\mu }_2}{\partial T_2},\\ J_{13}=p_2{v}_{n2}{C}_1\frac{\partial {\mu }_2}{\partial T_2}-C_1^2{R}_0,\\ J_{21}=C_1\frac{\partial h_2}{\partial p_2}-v_{n2},\\ J_{22}=v_{n2}^2\left({\mu }_2+p_2\frac{\partial {\mu }_2}{\partial p_2}\right)-p_2{\mu }_2\frac{\partial h_2}{\partial p_2},\\ J_{23}=p_2{\mu }_2-C_1{v}_{n2}\left({\mu }_2+p_2\frac{\partial {\mu }_2}{\partial p_2}\right);\\ J_{31}=\frac{\partial h_2}{\partial T_2},\\ J_{32}=\frac{\partial h_2}{\partial p_2}\left(p_2{v}_{n2}\frac{\partial {\mu }_2}{\partial T_2}-C_1{R}_0\right)-v_{n2}\left({\mu }_2+p_2\frac{\partial {\mu }_2}{\partial p_2}\right)\frac{\partial h_2}{\partial T_2},\\ J_{33}=C_1{R}_0-p_2{v}_{n2}\frac{\partial {\mu }_2}{\partial T_2}.\end{array}$$

The matrix determinant is found from the relation

$$\begin{array}{l}\mid J\mid =p_2{v}_{n2}{C}_1\left(\frac{\partial h_2}{\partial p_2}\frac{\partial {\mu }_2}{\partial T_2}-\frac{\partial {\mu }_2}{\partial p_2}\frac{\partial h_2}{\partial T_2}\right)+{\mu }_2\frac{\partial h_2}{\partial T_2}\left(p_2-C_1{v}_{n2}\right)+\\ +C_1{R}_0\left(v_{n2}-C_1\frac{\partial h_2}{\partial p_2}\right)-p_2{v}_{n2}^2\frac{\partial {\mu }_2}{\partial T_2}.\end{array}$$

The initial approximation is taken from the solution for a perfect gas.

To find partial derivatives of enthalpy and molecular weight with respect to pressure and temperature, finite difference formulas are used

$$\begin{array}{l}\frac{\partial h}{\partial p}=\frac{h\left(p+\mathrm{\Delta }p,T\right)-h\left(p-\mathrm{\Delta }p,T\right)}{2\mathrm{\Delta }p};\\ \frac{\partial h}{\partial T}=\frac{h\left(p,T+\mathrm{\Delta }T\right)-h\left(p,T-\mathrm{\Delta }T\right)}{2\mathrm{\Delta }T};\\ \frac{\partial \mu }{\partial p}=\frac{\mu \left(p+\mathrm{\Delta }p,T\right)-\mu \left(p-\mathrm{\Delta }p,T\right)}{2\mathrm{\Delta }p};\\ \frac{\partial \mu }{\partial T}=\frac{\mu \left(p,T+\mathrm{\Delta }T\right)-\mu \left(p,T-\mathrm{\Delta }T\right)}{2\mathrm{\Delta }T}.\end{array}$$

Here,$\mathrm{\Delta }p=k_p{p}^n$ и$\mathrm{\Delta }T=k_T{T}^n$ . The coefficients k_p and k_T are in the range of 0.005–0.01.

4. Governing equations

The Euler equations are used to describe the unsteady flow of inviscid compressible gas. The Euler equations are written as

$$\begin{array}{l}\frac{\partial }{\partial t}\underset{V}{\oint }‍\mathbf{U}\mathit{dV}+\underset{\partial V}{\oint }‍\mathbf{F}·d\mathbf{S}=0,\end{array}$$E5

where t is time, U is vector of conservative flow quantities at point x at time t, F is vector of fluxes, V is control volume with boundary ∂V, dS=ndS is vector of area element dS to the boundary ∂V with external normal n. The vector of conservative flow quantities and vector of fluxes in equation (5) have the form

$$\begin{array}{l}\mathbf{U}=\left(\begin{array}{c}\rho \\ \rho \mathbf{v}\\ \mathit{\rho e}\end{array}\right),\mathbf{F}=\left(\begin{array}{c}\rho \mathbf{v}\\ \rho \mathbf{vv}+p\mathbf{I}\\ \left(\mathit{\rho e}+p\right)\mathbf{v}\end{array}\right).\end{array}$$

The specific total energy is equal to the sum of the internal energy ε (it includes the energies of translational motion, rotational, vibrational, and electronic excitation of atomic and molecular components of the gas mixture), and kinetic energy

$$\begin{array}{l}e=\epsilon +\frac{1}{2}{\left|\mathbf{v}\right|}^2.\end{array}$$

Here, ρ is density, p is pressure, v is velocity vector, I is unit tensor. For a thermally perfect gas, specific heat capacity and molar mass remain constant. In high-temperature air, when the processes of dissociation and ionization become noticeable, the dependences of these characteristics on the parameters of the medium become more complicated; in particular, the dependence on pressure and temperature begins to play a role (as a result of a change in the molecular composition).

Eq. (6) is supplemented by the equation of state

$$\begin{array}{l}p=\rho \frac{R_0}{M_{\mathrm{\Sigma }}\left(p,T\right)}T,\end{array}$$

where R₀ is universal gas constant, $M_{\mathrm{\Sigma }}$ is molecular weight. The specific enthalpy is $h=\epsilon +p/\rho$ . The specific total energy is defined as

$$\begin{array}{l}e=h+\frac{1}{2}{\left|\mathbf{v}\right|}^2-\frac{p}{\rho }.\end{array}$$

The specific enthalpy is

$$\begin{array}{l}h=h_0+\underset{T_0}{\overset{T}{\int }}‍c_p\mathit{dT},\end{array}$$

where h₀ is enthalpy of formation of a substance at temperature T₀.

The molecular weight, $M_{\mathrm{\Sigma }}$, the gas constant $R=R_0/M_{\mathrm{\Sigma }}$, the heat capacity at constant pressure, c_p, and the heat capacity at constant volume, c_v, of perfect gas is constant. Therefore, Mayer’s relationship ($c_p-c_v=R$) is applicable. The ratio of specific heat capacities is $\gamma =c_p/c_v$, and the relations $h=c_{p}T$ and $\epsilon =c_{v}T$ are applied to find the enthalpy and internal energy.

5. Numerical method

The finite volume method is used to find numerical solutions to governing equations. Eq. (6) is solved numerically with the finite volume method. The flow domain is divided into closed control volumes. The mesh value found at the center of the control volume, V_i, is the average integral value

$$\begin{array}{l}{\mathbf{U}}_i=\frac{1}{V_i}\underset{V_i}{\int }‍\mathbf{U}\mathit{dV}.\end{array}$$

The integral over the boundary of the control volume i is a sum of the products of flux vector at the centers of face j and face area, S_ij. Then, Eq. (6) is written in semi-discrete form

$$\begin{array}{l}\frac{d{\mathbf{U}}_i}{\mathit{dt}}+\frac{1}{V_i}\sum _j^{N_i}‍{\mathbf{F}}_{\mathit{ij}}{S}_{\mathit{ij}}=0,\end{array}$$E6

where V_i is the volume of the control volume i, F_ij is the flow vector from cell i to cell j at the center of the face of the control volume, S_ij is the face area j of the control volume i.

To discretize time derivative in Eq. (7), an explicit third-order Runge–Kutta scheme is used. There are different approaches to calculate convective flows on the edge of the control volume. In this case, standard schemes for calculating flows, for example, the widely used Roe scheme, lead to a loss of accuracy and divergence of the computational procedure. The Godunov scheme and the Rusanov scheme are used to discretize convective flows. The second order of approximation in space is applied. Flow quantities are interpolated from the center of the cells to the edge of control volume, and limiters are applied to restrict the gradient of the solution and to ensure the monotonicity of the scheme.

The fluxes in finite volume scheme (7) are calculated in the direction of the normal to the boundary. The flux through the edge of control volume is found as

$$\begin{array}{l}{\mathbf{F}}_{j+1/2}=\frac{1}{2}\left[\mathbf{F}\left({\mathbf{U}}_L\right)+\mathbf{F}\left({\mathbf{U}}_R\right)\right]-\frac{1}{2}\left|A\right|\left({\mathbf{U}}_R-{\mathbf{U}}_L\right),\end{array}$$

where $\left|A\right|=R\left|\mathrm{\Lambda }\right|L$ and $\mathrm{\Lambda }=diag\left\{v_n-a,v_n,v_n+a\right\}$ is diagonal matrix with Jacobian eigenvalues on the main diagonal.

For an approximate accounting of complex physical and chemical processes in real gases, a methodology has been developed for the effective adiabatic exponent, which makes it possible to decompose the complete problem of modeling high-speed flows into stages. This ensures the creation of a universal computing complex, structured into a number of autonomous segments, with the possibility of independent modification of their functional content, improvement of algorithms, and computer implementation.

The heat capacity at constant pressure is calculated using numerical differentiation (second-order central difference discretization)

$$\begin{array}{l}{c}_p\left(p,T\right)={\left(\frac{\partial h}{\partial T}\right)}_p=\frac{h\left(p,T+\mathrm{\Delta }T\right)-h\left(p,T-\mathrm{\Delta }T\right)}{2\mathrm{\Delta }T},\end{array}$$

where $\mathrm{\Delta }T=0.01T$ .

6. Transformation of variables

The fluid flow described by the Euler equations is determined by the vector of conservative variables

$$\begin{array}{l}\mathbf{U}=\left\{\rho ,\rho u_x,\rho u_y,\mathit{\rho e}\right\},\end{array}$$

where ρ is density, $\rho u_x$, $\rho u_y$ are momentums in x and y directions, $\mathit{\rho e}=\rho \left[\epsilon +\left(u_x^2+u_y^2\right)/2\right]$ is total energy.

In addition to the vector of conservative variables, the vector of physical (primitive) variables is used

$$\begin{array}{l}\mathbf{Q}=\left\{\rho ,u_x,u_y,p\right\},\end{array}$$

where ρ is density, $u_x$, $u_y$ are velocities in x and y directions, p is pressure.

Physical variables are expressed in terms of conservative variables. The direct transformation operation is not straightforward. However, the reverse transformation requires the solution of equations

$$\begin{array}{l}\left\{\begin{array}{l}\mathit{\rho e}=\mathit{\rho \epsilon }\left(p,T\right)+\rho \frac{v^2}{2}\\ p=\rho \frac{R_0}{M_{\sum ‍}\left(p,T\right)}T\end{array}\right.\Rightarrow \left\{\begin{array}{ll}& \epsilon \left(p,T\right)=\frac{\left(\mathit{\rho e}\right)}{\rho }-\frac{v^2}{2}\\ & p{M}_{\sum ‍}\left(p,T\right)-\rho R_{0}T=0\end{array}\right.\end{array}$$

Unknown quantities are p, T, ε=h–p/ρ. These equations are non-linear. To solve non-linear equations, Newton’s method is applied. The variables at the new iteration are found from values at the previous iteration

$$\begin{array}{l}{p}^{n+1}=p^n+\mathrm{\Delta }p,T^{n+1}=T^n+\mathrm{\Delta }T,\end{array}$$

where $\mathrm{\Delta }p$, $\mathrm{\Delta }T$ are increments of pressure and temperature. The superscript n denotes the iteration number.

The quantities ε and $M_{\sum ‍}$ are expanded in a Taylor series keeping terms of first order

$$\begin{array}{l}{\epsilon }^{n+1}\left(p,T\right)={\epsilon }^n+\frac{\partial \epsilon }{\partial p}\mathrm{\Delta }p+\frac{\partial {\epsilon }^n}{\partial T}\mathrm{\Delta }T;\\ M_{\sum ‍}^{n+1}\left(p,T\right)=M_{\sum ‍}^n+\frac{\partial M_{\sum ‍}^n}{\partial p}\mathrm{\Delta }p+\frac{\partial M_{\sum ‍}^n}{\partial T}\mathrm{\Delta }T.\end{array}$$

Then, equations take the form

$$\begin{array}{l}\frac{\partial {\mathrm{\epsilon }}^n}{\partial p}\mathrm{\Delta }p+\frac{\partial {\mathrm{\epsilon }}^n}{\partial T}\mathrm{\Delta }T=\frac{\left(\mathit{\rho e}\right)}{\rho }-\frac{v^2}{2}-{\mathrm{\epsilon }}^n;\\ \left(M_{\sum ‍}^n+p^n\frac{\partial M_{\sum ‍}^n}{\partial p}\right)\mathrm{\Delta }p+\left(p^n\frac{\partial M_{\sum ‍}^n}{\partial T}-\rho R_0\right)\mathrm{\Delta }T=\rho R_0{T}^n-p^n{M}_{\sum ‍}^n.\end{array}$$

The increments of the pressure and temperature, Δp and ΔT, are obtained. The iterations finish when Δp and ΔT satisfy the accuracy conditions. The derivatives of ε and $M_{\sum ‍}$ are found with the central difference formulas of the second order

$$\begin{array}{l}\frac{\partial {\mathrm{\epsilon }}^n}{\partial p}=\frac{\mathrm{\epsilon }\left(p^n+{\delta }_p,T^n\right)-\mathrm{\epsilon }\left(p^n-{\delta }_p,T^n\right)}{2{\delta }_p};\\ \frac{\partial {\mathrm{\epsilon }}^n}{\partial T}=\frac{\epsilon \left(p^n,T^n+{\delta }_T\right)-\epsilon \left(p^n,T^n-{\delta }_T\right)}{2{\delta }_T};\\ \frac{\partial M_{\sum ‍}^n}{\partial p}=\frac{M_{\sum ‍}\left(p^n+{\delta }_p,T^n\right)-M_{\sum ‍}\left(p^n-{\delta }_p,T^n\right)}{2{\delta }_p};\\ \frac{\partial M_{\sum ‍}^n}{\partial T}=\frac{M_{\sum ‍}\left(p^n,T^n+{\delta }_T\right)-M_{\sum ‍}\left(p^n,T^n-{\delta }_T\right)}{2{\delta }_T}.\end{array}$$

Here, ${\delta }_p=k_p{p}^n$, ${\delta }_T=k_T{T}^n$. The coefficients k_p and k_T are chosen from the interval 0.005–0.01.

7. Results and discussion

A supersonic flow around a wedge with a half-angle β=30° by a perfect and high-temperature gas is simulated. The inlet pressure and inlet temperature (flow quantities before the shock) are 10⁵ Pa and 290 K. The working substance is air (γ=1.4 and μ=0.029 kg/mol). The inlet Mach number varies from 2 to 16. For the ideal perfect gas, the solution is available in a tabular form. The results presented in the study correspond to two cases. The difference between them is the velocity behind the shock wave front. It equals 10³ m/s (Case 1) and 3∙10³ m/s (Case 2). In Case 1, density is ρ=5.5 kg/m³, pressure is p=16.78 bar, temperature is T=1068 K for a perfect gas, and density is ρ=5.7 kg/m³, pressure is p=16.59 bar, temperature is T=1012 K for high-temperature air. In case 2, density is ρ=7.0 kg/m³, pressure is p=133.5 bar, temperature is T=6649 K for perfect gas and density is ρ=10.0 kg/m³, pressure is p=127.1 bar, temperature is T=4251 K for high-temperature air. In Case 3, the flow velocity is 3∙10³ m/s, and the inlet conditions for real and perfect gas are identical (density is ρ=7.0 kg/m³, pressure is p=134.5 bar, temperature is T=6649 K).

The mesh contains 110×160 nodes. Mesh nodes are clustered near the solid boundaries and shock wave front to take into account gradient regions of flow (Figure 2). The minimum residual level is used as a criterion for the convergence of the difference solution to the stationary solution of the problem. Approximately 2200 time steps are taken to achieve the specified residual level (in the calculations $R=10^{-10}$).

Figure 2.

Mesh.

The pressure distributions found from the perfect and real gas models are shown in Figure 3 at different times. In this case, the shock-wave structure for both models is similar. However, the compressed region for a real gas is slightly smaller than in the case of the perfect gas model. The temperature in case 1 does not exceed 1900 K (Figure 4). Temperatures are low, and there are no chemical reactions. Therefore, the molar mass of air remains constant. It should be noted that narrow regions with high temperatures exist where the temperature reaches 2480–3100 K, but this has little effect on the flow pattern.

Figure 3.

Case 1. Pressure contours at time 0.697 (a, b), 1.384 (c, d), 2.054 ms (e, f) for perfect (a, c, e) and real (b, d, f) gas.

Figure 4.

Case 1. Temperature contours at time 2.054 ms for a perfect gas.

Large flow velocity leads to significant differences in flow quantity distributions computed with the perfect gas model and real gas model (Figures 5 and 6). The pressure computed with a real gas model exceeds the pressure computed with a perfect gas model. The shock-wave structure in a real gas has a flattened shape in comparison with a perfect gas model. Dissociation and ionization processes in high-speed flow lead to different distributions of temperature. The maximum temperature in a real gas (it is about 11,000 K) is two times lower than the temperature computed with a perfect gas model (it is about 23,000 K). Temperature distributions computed with two models are compared in Figure 7. Density distribution is similar to pressure distribution, however, density computed with a perfect gas model is two times smaller than those computed with a real gas model.

Figure 5.

Case 2. Pressure contours at time (a, b) 0.242), (c, d) 0.476, and (e, f) 0.708 ms for perfect (a, c, e) and real (b, d, f) gas.

Figure 6.

Case 2. Temperature contours at time 0.708 ms for perfect (a) and real (b) gas.

Figure 7.

Case 2. Temperature contours at time 0.749 ms for a perfect gas (a) and at time 0.723 ms (b) for a real gas.

The distributions of the flow characteristics along the lower wall of the computational domain are shown in Figure 8 shows flow quantities distributions along x=0 line. When high-temperature effects in the air are taken into account, pressure undergoes relatively small changes. Dashed lines correspond to a perfect gas model, and solid lines correspond to a real gas model.

Figure 8.

Case 2. Pressure (a), temperature (b) and density (c) distributions along x = 0 line.

Distributions of flow quantities computed at the same inlet conditions show that the shock wave structures computed with different models of air are similar to each other. However, the distribution flow quantities are different (Figure 9). In a region of shock, the difference in the values of the parameters is small and similar to that observed in a perfect gas. At the same time, the temperature distributions computed with perfect and real gas models are different (Figure 10). Temperature distributions computed with perfect and real gas models are compared in Figure 11.

Figure 9.

Case 3. Pressure contours at time 0.245 (a), 0.478 (b), 0.711 ms (c) for a real gas.

Figure 10.

Case 3. Temperature contours at time 0.711 ms for a real gas.

Figure 11.

Case 2. Temperature contours at time 0.749 ms for a perfect gas (a) and at time 0.723 ms for a real gas (b).

Flow quantities distributions along line x = 0 are presented in Figure 12. High-temperature effects in the air have a significant impact on density and temperature distributions. At the same time, pressure distributions are relatively weakly affected by physical and chemical processes. Dashed lines correspond to a perfect gas model, and solid lines correspond to a real gas model.

Figure 12.

Case 3. Pressure (a), temperature (b), and density (c) distributions along x = 0 line for perfect and real gas models.

The influence of wedge angle and inlet Mach number is shown in Figure 13 (β=30°). Pressure distribution is not affected by high-temperature effects in air. At the same time, temperature distributions computed with different gas models are different. Dashed lines correspond to a perfect gas model, and solid lines correspond to a real gas model.

Figure 13.

Pressure (a), temperature (b), density (c), and shock angle (d) distributions computed with the perfect gas model and high-temperature model based on the Kraiko model [18].

For comparison, Figure 14 shows the distributions of flow characteristics behind a normal shock as a function of the inlet Mach number. The flow velocity behind a normal shock is subsonic. Therefore, a difference between the flow quantities computed with perfect and real gas models exceeds the mismatch of the flow quantities observed behind an oblique shock wave. Dashed lines correspond to a perfect gas model, and solid lines correspond to a real gas model.

Figure 14.

Pressure (a), temperature (b), and density (c) distributions behind a normal shock computed with the perfect gas model and Kraiko model [18].

The use of the perfect gas model at high Mach numbers of the oncoming flow leads to the inaccurate solutions. The relative error of computations of flow quantities obtained with the perfect gas model and real gas model is shown in Figure 15. An increase in the angle of the wedge leads to an increase in the error between solutions computed with various models of air.

Figure 15.

Relative error in determining pressure (a), temperature (b) and density (c) at β=10° (line 1), 20° (line 2), 30° (line 3). Line 4 corresponds to a normal shock wave.

8. Conclusion

At high intensity of the shock wave, which occurs at high supersonic and hypersonic flow velocities, the properties of the gas differ from the properties of a perfect gas. This leads to significant differences in the distributions of flow characteristics behind the shock wave front, corresponding to the models of a perfect and real gas. An approach and a calculation module have been developed that allow taking into account equilibrium chemical reactions in the air at high temperatures. To demonstrate the capabilities of the developed model, the problem of supersonic flow around a wedge with an attached shock wave is used. A comparison is made of the space-time distributions of flow characteristics calculated using the perfect and real gas models.

The developed computational module allows the inclusion in the design systems of advanced aircraft shapes, as well as integration with both commercial and open-source CFD packages.