The thermochemical non-equilibrium scale effects of the high enthalpy nozzle

The high enthalpy nozzle converts the high enthalpy stagnation gas into the hypervelocity free flow. The flow region of the high enthalpy nozzle consists of three parts: an equilibrium region upstream of the throat, a non-equilibrium region near the throat, and a frozen region downstream of the throat. Here we propose to consider the thermochemical non-equilibrium scale effects in the high enthalpy nozzle. By numerically solving axisymmetric compressible Navier-Stokes equations coupling with Park’s two-temperature model, the fully non-equilibrium solution is employed throughout the entire nozzle. Calculations are performed at different stagnation conditions with the different absolute scales and expansion ratio. The results of this study are twofold. Firstly, as the absolute scale and expansion ratio increase, the freezing position is delayed, and the flow approaches equilibrium. Secondly, the vibrational temperature and Mach number decrease with the increase in the nozzle scale and expansion ratio, while the speed of sound, static pressure, and translational temperature increase as the nozzle scale and expansion ratio increase.

The flow regimes in hypervelocity nozzle flow exhibit two types of scale effects [1,12]. One of these scale effects, caused by viscosity, is confined to the boundary layer near the nozzle wall and can be avoided when simulating on the ground. The other is the thermochemical non-equilibrium scale effects, which must be taken seriously as it affects the entire flow field including the core flow region; it is this scale effect that is the focus of this article. The high enthalpy convergent-divergent nozzle flow is a complex process, especially nozzle divergence, because of the large range of Reynolds number scale and the high-temperature effects [13]. As evidenced in previous studies, there are three distinct flow regions: an upstream near the equilibrium region, a non-equilibrium region near the throat, and frozen flow not far downstream of the throat [12]. When the flow is nearly frozen, there are no longer the non-equilibrium scale effects in the expansion zone. The flow parameters and chemical species depend on the dissociation of non-equilibrium before the freezing point.
This paper investigates the thermochemical non-equilibrium scale effects of the high enthalpy nozzle by the numerical method. To verify the reliability of the numerical method, the HEG and TCM2 nozzles are used to calibrate. Because they have undergone many tests and numerical validation. Based on the HEG nozzle, the flow field of the halved normal scale, normal scale, doubled normal scale and quintupled normal scale is numerically calculated. This paper also calculates the halved normal scale and normal scale of the FD-21 nozzle and the different expansion ratio based on the FD-21 nozzle. FD-21 was built in China Academy of Aerospace Aerodynamics in 2016 and is the largest-scale free-piston shock tunnel in the world [14,15].
At present, the three models to calculate the high enthalpy flow field [1,6,8,11]. The calculation results show that the flow field pressure, temperature, velocity and distribution of major chemical species are broadly coincident [4,5,11]. The ions of 7-species and 11-species, more than 5species belong to the microspecies, and the chemical energy contained in them is negligible relative to the total gas flow, so they do not affect the main flow and the thermochemical parameters [12]. When it comes to questions such as flow field ionization and radiation characteristics, it is necessary to use the 7-species or 11-species models. The 5-species and 7-species are used to calculate the high enthalpy flowfield in this paper. The axisymmetric compressible Navier-Stokes equations with Park's twotemperature model [16] are solved with a multi-block finite volume method. The results show that the freezing point moves downstream as the nozzle scale and expansion ratio increase. The flowfield characteristics, such as static temperature, Mach number, and species mass concentration distribution, are related to its absolute position, namely the nozzle scale. The larger nozzle scale and expansion ration can effectively suppress thermochemical non-equilibrium flow, which are more suitable for simulating the flight environment.

Governing equations and numerical method
The unsteady dimensionless axisymmetric form of the Navier-Stokes equation for the 5-species and 7-species air in the thermochemical non-equilibrium state can be expressed in the vector equation [17] as: where the vector U is the vector of conservative quantities. The vectors F and G represent the convective fluxes. The viscous terms F v and G v contain the transport term.
The vectors M and M v are the inviscid and viscous source terms. The vector W represents chemical and vibrational terms. The term Re is the Reynolds number. These quantities are given by the following expression: where ns is the number of species, which is equal to 5 or 7. ρ is the density of the mixture. C i is species mass fraction. u and v are the components of the velocity V. r is the coordinate in the radial direction. Le i is the Lewis number of the species. Pr is the Prandtl number of the mixture. k is the Coefficient of heat conduction. C p is the specific heat. τ xx , τ rr , τ xr , and τ θθ are the components of the viscous stress tensor. w i is the chemical source terms. H is the total enthalpy. ρE is the total energy per unit volume, containing translational energy ρe t , rotational energy ρe r , vibrational energy ρe v , electronic energy ρe e , free energy ρe 0 , and the kinetic energy of the gas ρV 2 /2. ρE is expressed as: ρC i e ti þ e ri þ e ei þ e 0i ð Þ þ where nt is the number of species that remove electrons from ns. q x , q r and q vx , q vr are the components of the total heat flow and the vibrational heat flow in the x and r coordinates, respectively. They are expressed as: S v is the vibrational source terms of the species. It is expressed as:

Non-equilibrium model
For mixed gas, each species is considered as a perfect gas. The chemical source terms are derived from several reactions. The following equations are the 5-species air with 17-reaction scheme and 7-species air with 21-reaction scheme [18]. The coupling between chemical and vibrational processes is taken into account. The forward reaction rate k f,r , and backward reaction rate k b,r of the dissociation equations are expressed in the modified Arrhenius form [17] as where A f,r and A b,f are the forward and backward reaction rate constant. T Df,r and T Db,r are the characteristic reaction temperatures for the forward and backward reactions, respectively. Values for the reaction rates k f,r and k b,r are taken from reference [17] for the 5-species and 7-species air model. T av is the geometric average temperature [16]: where T is the translational temperature, and T V is vibrational temperature. The translational temperature T is equal to the static temperature T.
The mixture viscosity and thermal conductivity of the mixture are computed by applying Gupta-Yos viscosity fits [17] for each species.
where A μ,s , B μ,s , C μ,s , A k,s , B k,s , C k,s , D k,s and E k,s are constants. Wilke's Law [19] is then used to compute the bulk quantities. For temperatures between 300 K and 8000 K, the vibrational relaxation time τ i is given by the following equation [20].
where τ i,P is translational-vibrational energy relaxation time (TVERT) for the species s from the report by Park [20]. τ i,MW is TVERT for the molecular species s from the correction by Millikan and White [20].
Where P is pressure. N i is number density. M si is reduced molecular weight of species. A s is constant.
For a perfect gas, the speed of sound is given by a 2 = (∂p/∂ρ) s . For the nonequilibrium sound speed [1] a n , the corresponding expression can be written as where subscripts denote partial derivatives of specific enthalpy h = h(p, ρ, C i ) for the variable: the pressure, density and specific mass fraction.
For the frozen sound speed a f , the gas assumes no reactions. a f is the frozen sound speed defined by The unsteady flow equations are integrated in time, which are solved for the convergent-divergent axisymmetric high enthalpy nozzle with given stagnation conditions and wall conditions. Jameson [21] fourth-order Runge-kutta method [21] is employed for time integration. The inviscid numerical fluxes are computed using a low-dissipation AUSM scheme [22], with a MUSCL scheme for a third-order spatial accuracy. For the physics of the strongly favorable pressure gradient around the nozzle throat, the Spalart-Allmaras one-equation turbulence model [23] with the Catris-Aupiox compressibility correction is evaluated for the task. The chain rule is used to calculate the viscosity flux [24]. The exit boundary is dominated by hypersonic flow by diverging free stream conditions generated by the nozzle. All conservative variables are extrapolated. No-slip conditions for the solid walls. The wall conditions are given for a fully catalytic wall and an isothermal wall where T w = 600 K.

Physical model and grid generation
As mentioned above, the TCM2 conical nozzle [4], the HEG conical nozzle [5], and the FD-21 conical nozzle [9] are studied. The TCM2 nozzle and HEG nozzle are calculated to verify the accuracy of the program. The HEG nozzle is halved normal scale (0.5NS), normal scale (NS), doubled normal scale (2NS) and quintupled normal scale (5NS). For the conical nozzle, the main dimensions contain the half cone angle (HCA), throat diameter (d t ), exit diameter (d e ) and the expansion ratio (ER, defined as the ratio of exit area to throat area). The normal scale TCM2 and HEG nozzles are recorded in Table 1. Cases 1 and 2 in Table 2 is the NS and 0.5NS of FD-21 nozzle, and Cases 3 and 4 change ER by changing the nozzle throat diameter.
A structured grid is used for all the numerical simulations. The grid consists of two blocks: one of the blocks is the shock tube part, which is 200 mm from the end of the shock tube in the axial direction; the other is the nozzle part. The structured grid is shown in Fig. 1. The total grid number is 251,250. The grid system of the shock tube is 51 × 201 and the grid system of the nozzle part contains 1200 × 201. The system is used with Δx min = 1 × 10 − 4 m located at the throat section, and Δy min = 2 × 10 − 8 m at the wall. Exponential compression is used in the radial direction. A grid convergence study is conducted with three different grid resolutions (125,100; 150,120 and 251,250), and a small change in pressure, species, velocity, and Mach number distribution is registered. For each grid system, the y + value is less than one.

Conditions of numerical simulation
The initial conditions are chosen from the data obtained by the TCM2 nozzle [4] and the HEG nozzle [5], which correspond to Case 1 and Case 2 in Table 3, respectively. The pressure, temperature and the species mass concentrations in stagnation conditions are obtained by equilibrium assumption and are reported in Table 3. The initial condition of FD-21 is also obtained from Case 2 in Table 3.

Precision estimates
Smirnov et al. [25,26] point out that estimating precision and errors accumulation is necessary for large scale simulations of complex chemically reacting gas dynamics. A definite error occurs in integration at each step, which is dependent on spatial resolution and numerical solver. This error can be expressed in absolute or relative values. Accumulation of error takes place for successive time steps. The relative error of integration in one-dimensional case S 1 is proportional to the mean ratio of cell size ΔL to the domain size L 1 in the direction of integration in power depending on scheme accuracy [25]: Where k is the order of accuracy of numerical scheme. In the presence of several directions in integration, the errors [25] are written as The allowable value of total error S max [25] is defined by the user of the simulator. For initial and boundary conditions are usually not known with higher accuracy, S max should be between 1% and 5%. Then the following inequality should be satisfied.
Where n is the number of time steps in Navier-Stokes equations integration. Another important characteristic is introduced for the results of supercomputer simulations. The ratio of the maximal allowable number of time steps to the actual number of time steps is defined as [25]: Table 4 shows the accumulation of error for the grid resolution and physical time simulations. The allowable error was assumed at 5%. For present numerical simulations, all results demonstrate high reliability, but it is not always appropriate to simulate longer steps.

Verification of the calculation method
To verify the calculation method, the comparisons with the TCM2 [4] and HEG [5] results are made. For the TCM2 case, the initial condition is based on Case 1 in Table 3, and the normal scale is selected. Figure 2 shows the centerline vibrational temperature T V and translational temperature T distribution. The T V,N2 is the vibrational temperature of N 2 and is about 2400 K, and T V,O2 is the vibrational temperature of O 2 and is about 1000 K. The unified vibrational temperature Tv is about 2290 K in this paper. Considering the mass ratio of N 2 to O 2 in the flow, this value is reasonable. The centerline evolution of vibrational temperature and translational temperature are plotted, as shown in Fig. 3. The results of the freezing position and temperature are in good agreement with the reference. Figure 4 shows the centerline Mach number distribution, and the variation of about 5% is predicted. The discrepancy can be explained by the different vibrational modeling, resulting in different translational temperature T (static temperature). Table 5 shows the comparison of the Mach number and species mass    fraction in the nozzle exit. The accuracy of the static pressure P and translational temperature T is verified by the HEG nozzle. For the HEG case, the initial condition is based on Case 2 in Table 3. The calculation results are compared with the results for the reference [5], as shown in Table 6. The calculation results of the nozzles of TCM2 and HEG in this paper agree well with the reference values [4,5]. The discrepancy calculated by comparing the main reference values is nominally 5% or less for all the cases considered. The variation can be explained. There are three main reasons for these differences. First, the vibrational modeling modifies the translational temperature, as deduced from Eq. (10). Second, there are the different expressions of the transport coefficients. The species NO + and e − belong to the micro-species. The results verify the positive determination of the relevant parameter settings and the accuracy of the procedure. Third, the turbulence energy could decrease some other energy in an isolated system. When the total temperature exceeds 6000 K, there are trace species such as electrons, NO+. It is more appropriate to use 7-species, so the following calculation uses 7-species.      (1) As the nozzle scale increases, the pressure gradient around the throat decreases (see Fig. 7). The pressure downstream of the throat increases, which shortens TVERT, as deduced from Eq. (28). Additionally, the decrease of pressure gradient makes the flow velocity decrease and the flow time increase. Both lengthen the transition from equilibrium to freezing, and the flow is close to the equilibrium flow. (2) Since the nozzles use the same inlet conditions, the velocity u at the nozzle exit is almost identical [27], which is consistent with h ≈ u 2 /2. Because of the freezing of vibrational energy, the increase of the kinetic energy of the airflow in nozzle mainly comes from the translational energy, which leads to the decrease of translational temperature T.   As can be seen from the above conclusion, the flow is frozen not far downstream of the high enthalpy flow throat. Cases 1,3 and 4 in Table 1 investigate the thermochemical nonequilibrium scale effects by changing the diameter of the throat alone. Figure 8 shows the distribution of vibrational temperature and translational temperature along the centerline of nozzle. For the throat radius d t is 20 mm, 40 mm, and 60 mm, the freezing temperatures are 2442 K, 2207 K, and 2140 K, respectively. The distances x f from the freezing point to the throat are 0.717 m, 1.568 m, and 3.030 m. The freezing point diameters d f of their corresponding nozzles are 0.254 m, 0.536 m, and 1.018 m.
The ratio d f /d t of the freezing point diameter to the throat diameter is 12.7, 13.4, and 16.9, respectively. As the throat radius increases, the freezing point moves away from the throat. They can be explained in the following aspects.
(1) The effects of the throat diameter are the same as the scale effects of the nozzle.
As the nozzle diameter increases, the pressure gradient around the throat decreases (see Fig. 9), and the flow velocity around the throat decreases (See Fig. 10) and the flow time increases. Both lengthen the transition from equilibrium to freezing, and the flow is close to the equilibrium flow. (2) In thermal equilibrium, the increase of the kinetic energy of the gas comes from both molecular translational energy and molecular vibrational energy. Therefore, thermal equilibrium causes the translational temperature and vibrational temperature to decrease. When the vibrational energy freezes, the increase in the kinetic energy of the gas comes only from the translational kinetic energy. The nozzle outlet speed is almost unchanged, which makes the reduction of translational energy deeper and the translational temperature lower.   Figure 11 shows the Mach number distribution along the centerline. Before the freezing point, the value and evolution of the Mach number are the same. After the freezing point, this phenomenon has changed. It can be explained by the freezing of the vibrational levels. As the nozzle scale increases, the vibrational temperature decreases and the translational temperature increases. With the increase of the nozzle scale, the flow gradually approaches the thermodynamic equilibrium state. As the translational temperature increases, the local sound speed increases (see Table 5), resulting in the decrease in the Mach number. Figure 12 can also prove this. The larger the size is, the larger the ratio of the exit uniform area radius to the exit radius is, which will make the exit Mach number smaller. This is in contradiction with the assumption that the thicker the boundary wall of the nozzle wall and the smaller the Mach number under the same area ratio. It can be known by combining Figs. 11 and 12. The influence of thermochemical non-equilibrium scale effects on Mach number is more important than the viscous scale effects. Cases 1, 3 and 4 in Table 2 are employed to study the effect of the variation throat diameter on the Mach number, as shown in Fig. 13. When the throat diameter increases, the flow is close to the equilibrium state. More vibrational energy is transmitted to the translational energy, which makes the translational temperature rise, as deduced from the kinetic energy in Eq. (10). The change in Mach number due to ER and non-equilibrium effects is much more pronounced than the change in Mach number caused by ER alone. Figures 14 and 15 show the changes of O 2 , O and NO in the axial direction based on the HEG nozzle, respectively. The positions where these components change are the   upstream of the freezing point, which is consistent with the changing trend of the vibrational temperature. The flow parameters depend on the dissociation of nonequilibrium before the freezing point. Tables 7 and 8 show the main flowfield parameters in the centerline of the nozzles exit, and the variation trend of the flowfield parameters can be seen. The nozzle flowfield parameters are monotonous to the absolute scale and throat diameter, such as pressure, temperature, sound velocity, Mach number, and species mass fraction. The larger the nozzle scale and throat diameter are, the larger the mass fraction of the species N 2 and O 2 is, but the smaller the mass fraction of the species NO, O, N, NO + is. The species N and NO + are trace species, which are not listed here. As can be seen from Tables 7 and 8, when the nozzle scale increases, the degree of the thermal non-equilibrium decreases and the frozen vibrational energy decreases, so that the kinetic energy of the airflow, namely the speed u ∞ , increases slightly.

Conclusions
As analyzed in the paper, the high enthalpy nozzle has the thermochemical nonequilibrium scale effects, which have a certain influence on the flow field parameters of  the nozzle exit. For a given half cone angle and expansion ratio, the thermochemical non-equilibrium scale effects of the high enthalpy nozzle mainly have the following consequences: (1) As the nozzle scale and expansion ratio increase, the degree of the thermal nonequilibrium decreases, and the frozen vibrational energy decreases, and translational temperature increases. The larger nozzle scale and nozzle expansion ratio can effectively suppress the thermochemical non-equilibrium flow, and is more suitable for simulating the flight environment. 6 Nomenclature 0.5NS halved normal scale NS Normal scale 2NS doubled normal scale 5NS quintupled normal scale 5S2T 5-species two-temperature model 7S2T 7-species two-temperature model ns the number of species, which is equal to 5 or 7 nt the number of species that remove electrons from ns JS conceived of the study, designed the study and collected the data. All authors analysed the data and were involved in writing the manuscript. All authors read and approved the final manuscript.