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. With numerical 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 contain 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.

its scale effects 1 . A large deal of research focuses on the scale effects of the test model and the actual aircraft 1,2 , but little research has been published on the scale effects of the high enthalpy nozzle 3 , especially the thermochemical non-equilibrium scale effects on the high enthalpy tunnel simulation ability. There are many high enthalpy tunnels in the world, such as TCM2 4 , HEG 5 ,T4 6 , T5 6 , HIEST 6 , JF-10 7 , LENS I 8 , and FD-21 9 . Due to the complexity and difficulty of the hypervelocity simulation, a great number of national and international space and hypersonic flight projects have been carried out extensively, such as integrated scramjet configuration 10 , complex hypersonic flight configuration 10 , and hypersonic boundary layer transition 11 . The experimental data obtained from the same test model are some differences under the similar stagnation gas and boundary condition. Except for the factors of the flowfield quality of the high enthalpy tunnel, the diameters of the high enthalpy nozzles range from 0.4m to 2.0m.
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 effects 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 nozzles scale when halved normal scale (0.5NS), normal scale (NS), doubled normal scale (2NS) and quintupled normal scale (5NS) are numerically calculated. This paper also calculates the 0.5NS halved normal scale and NS 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 of 5-species(O 2 , N 2 , NO, O, N), 7-species(O 2 , N 2 , NO, O, N, NO + , e -), and 11-species(O 2 , N 2 , NO, O, N, NO + , O + , N + , O 2 + , N 2 + , e -) are used 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 5-species 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 two-temperature 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 ER 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 ER 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 5species and 7-species air in the thermochemical non-equilibrium state can be expressed in the vector equation 18 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 S is species mass fraction. u and v are the components of the velocity V. r is the coordinate in the radial direction. Les is the Lewis number of the species. Pr is the Prandtl number of the mixture. k is the Coefficient of heat conduction. Cp is the specific heat. τ xx , τ rr , τ xr , and τ θθ are the components of the viscous stress tensor. w s 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: 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 17 . 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 18 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 18 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 18 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 300K and 8000K, 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 reported by Park 20 . τ i,MW is TVERT for the molecular species s from the correction by Millikan and White 20 .
For a perfect gas, the speed of sound is given by . For the non-equilibrium sound speed 1 a n , the corresponding expression can be written as where subscripts denote partial derivatives of specific enthalpy 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 =600K.

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 1a. Cases 1 and 2 in Table 1b 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 200mm 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 251250. 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 (125100; 150120 and 251250), 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 Case1 and Case2 in Table 2, respectively. The pressure, temperature and the species mass concentrations in stagnation conditions are obtained by equilibrium assumption and are reported in Table 2. The initial condition of FD-21 is also obtained from case2 in Table 2.

Verification of Calculation Procedure
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 2, and the normal scale is selected. Fig.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 2400K, and T V,O2 is the vibrational temperature of O 2 and is about 1000K. The unified vibrational temperature Tv is about 2290K 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 good agreement with the reference. Fig.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  3a 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 Case2 in Table 2. The calculation results are compared with the results for the reference 5 , as shown in Table 3b.     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 percent or less for all the cases considered. The variation can be explained. There are three main reasons for these differences. First, by the different vibrational modeling. Second, by the different 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 6000K, there are trace species such as electrons,NO+. It is more appropriate to use 7-species, so the following calculation uses 7-species.

Freezing point position
Numerical results of the different scale nozzles are now presented. The centerline evolution of the vibrational temperature T V and translational temperature T are plotted in Fig.5. Not far downstream of the nozzle throat, the larger-scale nozzle has a slightly higher translational temperature. The vibrational temperature modifies the translational temperature, as deduced from the kinetic energy in Eq. (10). Compared with the translational temperature, the vibrational temperature is more affected by the absolute scale. The larger the scale, the further away from the throat the freezing point, and the lower the vibrational temperature. The vibrational temperature of 0.5NS, NS, 2NS, 5NS is 2613, 2298, 2036, 1755, respectively. In this initial condition, the absolute scale of the nozzle is increased tenfold and the vibrational temperature is reduced by 858K, as shown in Fig.5. The freezing point moves down the nozzle throat as the nozzle scale increases. At the centerline vibrational temperature in the different nozzle scales (from 0.5NS to 5NS), the distances x f from the freezing point to the throat are 451mm, 1241mm, 3091mm, and 9712mm. The freezing point diameters d f of their corresponding nozzles are 15.7mm, 39.3mm, 92.2mm, and 276.1mm. The ratio d f /d t of the freezing point diameter to the throat diameter is 2.86, 3.57, 4.19, and 5.02, respectively. The thermochemical non-equilibrium scale effects are also suitable for the nozzle of FD-21. Fig.7 is a comparison of 0.5NS and NS flow field temperatures. At the centerline vibrational temperature in the 0.5NS and NS nozzles, the distances x f are 281mm and 875mm. The freezing point diameters d f of their corresponding nozzles are 98.1mm and 254.1mm, and the ratios d f /d t are 9.81 and 12.71, respectively. The absolute scale not only affects the translational temperature and vibrational temperature of the nozzle flow field but also the freezing point. Those phenomena can be explained in the following aspects.
(1) As the nozzle scale increases, the pressure gradient around the throat decreases (see Fig.6). 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 25 , 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. Fig.8 shows the distribution of vibrational temperature and translational temperature along the centerline of nozzle. For the throat radius d t are 20mm, 40mm, and 60mm, the freezing temperatures are 2442K, 2207K, and 2140K, respectively. The distances x f from the freezing point to the throat are 0.717m, 1.568m, and 3.030m. The freezing point diameters d f of their corresponding nozzles are 0.254m, 0.536m, and 1.018m. 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 decrease(See Fig.10) and the flow time increase. 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.  Fig.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 increase. 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 4a), resulting in the decrease in the Mach number. Fig. 12 can also prove this. The larger the size, the larger the ratio of the exit uniform area radius to the exit radius, 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 Fig.11 and Fig.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 1b 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. Figs. 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 non-equilibrium before the freezing point. Tables 4a and 4b 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, the larger the mass fraction of the species N 2 and O 2 , but the smaller the mass fraction of the species NO, O, N, NO + . The species N and NO + are trace species, which are not listed here.
As can be seen from Table 4a, 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 non-equilibrium scale effects, which has a certain influence on the flow field parameters of the nozzle exit. For a given HCA half cone angle and ER expansion ratio, the thermochemical non-equilibrium scale effects of the high enthalpy nozzle mainly have the following consequences: (1) As the nozzle scale and ER expansion ratio increases, the degree of the thermal non-equilibrium decreases, and the frozen vibrational energy decreases, and translational temperature increases. The larger nozzle scale and nozzle ER expansion ratio can effectively suppress the thermochemical non-equilibrium flow, and is more suitable for simulating the flight environment.
(2) The absolute scale affects the freezing point of the high enthalpy nozzle. The larger absolute scale and ER expansion ratio cause the freezing state of the vibrational dynamics to be delayed. The freezing point moves down the nozzle throat as the nozzle scale and ER expansion ratio increase.
(3) As the nozzle scale increases, the static pressure, and translational temperature at the nozzle exit increase slightly, but the Mach number and vibrational temperature decrease significantly. The larger the nozzle scale and ER expansion ratio, the larger the mass fraction of the species N 2 and O 2 , but the smaller the mass fraction of the species O, N, NO, etc.