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

The high enthalpy shock tunnel can simulate the free-flow speed above 3km/s. The characteristic of the flow is that the kinetic energy of the high enthalpy stagnation gas is high enough to effectuate high-temperature effects such as dissociation even ionization of fluid molecules. The high enthalpy nozzle converts the high enthalpy stagnation gas into hypervelocity free flow. The flow of the high enthalpy nozzle consists of three distinct flow regions: an equilibrium region upstream of the throat, a non-equilibrium region near the throat, and a frozen region downstream of the throat. The study focuses on the conical nozzle, testing thermochemical non-equilibrium expansion effects under the different expansion angle of the expansion section, the curvature radius of the throat, the throat radius, and the convergence angle of the convergent section. A multi-block solver for axisymmetric compressible Navier-Stokes equations is applied to simulate the thermochemical non-equilibrium flow in several high enthalpy conical nozzles. The significant conclusions of this study contain tripartite. Firstly , the thermochemical non-equilibrium effects are sensitive to the maximum expansion angle and throat radius, but not to the radius of throat curvature and the contraction angle. Secondly, as the maximum expansion angle decreases and the throat radius increases, the flow approaches equilibrium. Finally, the maximum expansion angle and the throat radius not only affect the position of the freezing point but also impacts the flow field parameters, such as temperature, Mach number, and species mass concentration.


Introduction 1
During the re-entry flight of the space vehicle in the earth's atmosphere or the interplanetary atmospheric entry of space vehicles, the flow surrounding the aircraft is characterized by the high-temperature effects such as the dissociation even ionization of the fluid molecules in the flow past hypersonic vehicles 1,2,3 . The high enthalpy tunnel can duplicate the high flow velocities and subsequently the high-temperature effects 4,5 . For the high enthalpy shock tunnel, the stagnation gas passes through a convergent-divergent nozzle and converts into the hypervelocity free flow, where the stagnation gas is highly dissociated 5,6 . During the expansion of the nozzle, the velocity of the airflow increases rapidly, but the temperature and density of the airflow decreases rapidly. The flow can't reach the equilibrium under the local temperature and pressure, and it tends to be thermochemically frozen not far from the throat 6,7 . The vibrational relaxation and recombination reaction can't return it to the state of the natural atmosphere. At the exit of the nozzle, the vibrational temperature of the airflow is very different from the translational temperature, and the chemical species has not reached the level of the natural atmosphere. In other words, the free flow of the high-enthalpy wind tunnel has a certain degree of thermochemical non-equilibrium 8,9 . This situation is undesirable, as it makes the test data of the nozzle exit airflow inconsistent with the flight environment, and can restrict the simulation capabilities of the high enthalpy nozzle. It is necessary to reduce the thermochemical non-equilibrium flow of the high enthalpy nozzle.
There are many high enthalpy shock tunnels in the world, such as T4 5 , T5 5 , HIEST 5 , TCM2 10 , HEG 11 , LENS I 12 , T6 13 , JF-10 14 , and FD-21 15 . The high enthalpy nozzles of these tunnels can be divided into the contour nozzles and the conical nozzles. The high enthalpy nozzles have different maximum expansion angle , ranging from 4 to 12 degrees. If the maximum expansion angle is too small, the nozzle will be too long, resulting in high cost. When the maximum expansion angle is too large, it is easy to cause the flow separation downstream of the throat, which affects the quality of the flow field. There are many studies on the high enthalpy nozzles themselves 6,7,8 , but few focuses on the thermochemical non-equilibrium expansion effects based on the different maximum expansion angle. Therefore, it is necessary to analyze the influence of the maximum expansion angle and another design parameters on the nozzle, such as the throat curvature radius, the throat radius, and the convergence angle. The high enthalpy contour nozzle has many factors to consider. If the design parameters are unreasonable, it will produce a large error 9,10,11 . Conical nozzles have fewer design parameters and are easier to design. Therefore, conical nozzles are suitable for studying the thermochemical non-equilibrium expansion effects.
At present, the high-temperature gas models are usually used to calculate the high enthalpy flow field. They mainly includes 5-species model(O 2 , N 2 , NO, O, N), 7-species model (O 2 , N 2 , NO, O, N, NO + , e -) and 11-species model (O 2 , N 2 , NO, O, N, NO + , O + , N + , O 2 + , N 2 + , e -). The calculation results 16,17,18 show that the pressure, temperature, density, velocity and distribution of the major chemical species are broadly coincident. The ions of the 7-species and 11-species, more than the 5-species, belong to the micro-species: the chemical energy contained in them is negligible relative to the total gas flow, so they do not affect the main flow and thermochemical parameters 16,17,18 . It is only necessary to use the 7-species or 11-species model when it comes to problems such as flow field ionization characteristics and radiation characteristics.
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 high enthalpy nozzles of the different maximum expansion angle of the expansion part, the throat curvature radius, the throat radius, and the convergence angle of the convergent part are studied. By comparing the nozzles with the different conditions, it can be concluded that the small maximum expansion angle and the large throat radius can effectively suppress the thermochemical non-equilibrium effects, but the throat curvature radius and the convergence angle have little influence on non-equilibrium flow. These conclusions can provide some references for the design of high enthalpy conical nozzles, especially high enthalpy contour nozzles.

Governing equations
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 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. 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:

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 .
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 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 . 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 . For the non-equilibrium sound speed a n and the frozen sound speed a f , 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. 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, Candler 23 and Wang 24 calculated the typical hypersonic nozzle flows and showed good agreement with measurable nozzle flow field using the of Spalart-Allmara one-equation turbulence model 25 with the Catris-Aupiox compressibility correction. The outflow boundaries are dominated by hypersonic flow by diverging free stream condition. 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 nozzle 10 , the HEG nozzle 11 , and the FD-21 nozzle 15 are studied. The TCM2 nozzle and HEG nozzle are calculated to verify the accuracy of the calculation method. Based on the FD-21 nozzle, the thermochemical non-equilibrium expansion effects are studied by changing the maximum expansion angle(MEA) of the expansion part, the throat curvature radius(TCR), the throat radius(TR), and the convergence angle(CA) of the convergent part.
For the conical nozzles, the main dimensions, such as CA θ c , TCR r c , TR r t , exit radius r e , MEA θ e , ER(expansion ratio) are listed in Tables 1 and 2. The parameters of the TCM2 nozzle and the HEG nozzle are shown in Table 1. Table 2 shows a total of 10 cases. Cases 1 to 4 are to verify the effect of the MEA θ e . Cases 3, 5 and 6 are used to verify the influence of changes in TR r t on thermochemical non-equilibrium effects. Cases 3, 7, and 8 are used to verify the influence of changes in CA θ c on thermochemical non-equilibrium effects. Cases 3, 9, and 10 is used to verify the influence of the change of TCR r c on the thermochemical non-equilibrium effects. For all cases based on the FD-21 nozzle, the inlet radius is 145mm and the exit radius is 1000mm.   A structured grid is used for all the numerical simulations, and the structured grid is shown in Fig.1. The grid consists of two blocks(see Fig.2): one of the blocks is the shock tube part, which is 300mm from the end of the shock tube. The other is the nozzle part.. 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 10 and the HEG nozzle 11 , which correspond to Case1 and Case2 in the Table 3, respectively. The pressure, temperature, and species fraction in reservoir conditions are obtained by equilibrium assumption and are reported in Table 3. The wall conditions are given for the fully catalytic wall and the isothermal wall where T w =600K. Table 3 Initial conditions

Precision estimates
For large scale simulations of complex chemically reacting gas dynamics, it's necessary to estimate precision and errors accumulation 26 . 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 6 .
The relative error of integration in one-dimensi onal case S 1 is proportional to the mean ratio of c ell size ΔL to the domain size L 1 in the direction of integration in power depending on scheme accur acy 27 : The allowable value of total error S max 27 is def ined by the user of the simulator. For initial and b oundary conditions are usually not known with hig her accuracy, S max should be between 1% and 5%. Then the following inequality should be satisfied.  Table 4 shows the accumulation of error for t he grid resolution and physical time simulations 6 . T he allowable error was assumed at 5%. For present numerical simulations, all results demonstrate high reliability, but it is not always appropriate to simu late longer steps.

Verification of the calculation method
The authors 6 use the 5-species 2-temperature(5S2T) model and the 7-species 2-temperature(7S2T) model to calculate the flow field of the TCM2 and HEG nozzle. For the TCM2 case and HEG case, the initial condition is based on Case1 and Case2 in Table 3 Tables 5 and  6 show the comparison of the Mach number and species mass fraction in the TCM2 and HEG nozzle exit.   Table 6 Comparison of flow parameters in the HEG nozzle 6 Results P (Pa) The calculation results 6 of the TCM2 nozzle and HEG nozzle in this paper agree well with the reference values 10,11 . The discrepancy calculated by comparing the reference values is nominally 5 percent or less for all the cases considered. The variation can be explained. There are three main reasons 6 for these differences. First, the vibrational modeling modifies the translational temperature, as deduced from Eq. (10). Second, there are the different expression of the transport coefficients.. The species eand NO + 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 e -, NO+. It is more appropriate to use 7-species, so the following calculation uses 7-species.

Discussion of Results
The different geometries in Table 2 are calculated, which are based on the FD-21 nozzle. The reservoir conditions for Case 2 are listed in Table 3.

Freezing point position
Comparisons of Cases 1, 2 3, and 4 reveal the influence of the MEA θ e on thermochemical non-equilibrium effects. The translational temperature T and vibrational temperature T V of the nozzle centerline are exhibited in Fig.6. It can be seen that the vibrational temperature T V and the translational temperature T are the same before the throat, and the gas is in the thermodynamic equilibrium state. After passing through the throat x=0m, the vibrational temperature T V and the translational temperature T begin to separated, and the vibrational non-equilibrium phenomenon occurs. Then, as the gas velocity gradually increases, the vibrational temperature freezes. Not far downstream of the nozzle throat, the larger MEA of the nozzle has higher vibrational temperature but slightly lower translational temperature. For the MEA are 4°, 6.5°, 9°, and 12°, freezing temperatures are 2018K, 2113K, 2442K, and 2623K. The distances x f from the freezing point to the throat are 3.071m, 1.576m, 0.717m, and 0.473m. The freezing point diameters d f of their corresponding nozzles are 0.438m, 0.376m, 0.254m, and 0.220m. The ratio d f /d t of the freezing point diameter to the throat diameter is 21.9, 18.8, 12.7, and 11, respectively. The freezing point approaches the nozzle throat as the MEA increases. Those phenomena can be explained in the following aspects.
(1) As the MEA increases, the pressure gradient near the throat increases (see Fig.7), resulting in the flow speed increases, as illustrated in Fig.8. As the pressure gradient increases, the pressure downstream of the throat decreases, which growths TVERT, as deduced from Equation (25). Additionally, the increase in pressure gradient makes the flow velocity increase but the flow time decrease. Both shorten the transition from equilibrium to freezing, and the flow is close to the non-equilibrium flow.
(2) The pressure p decreases with the increase of the MEA. According to the formula p=nkT, the translational temperature is directly related to the pressure p, and the vibrational freedom has no direct influence on the pressure. So the influence of the non-equilibrium on the pressure distribution is consistent with the trend of the influence on the translational temperature T.
(3) Not far downstream of the nozzle throat, the larger MEA θ e has a slightly larger vibrational temperature, resulting in a higher vibrational energy level. The vibrational model modifies the translational temperature, as deduced from the kinetic energy in Equation (10). Due to 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. Since the vibration energy is frozen, it can also be obtained that the contribution of vibrational energy to the speed is extremely small, which means that the non-equilibrium effect has almost no effect on the velocity u of the nozzle outlet. The velocity u is almost the same, which is consistent with h≈u 2 /2, and the velocity u is only related to the total enthalpy h of the nozzle stagnation gas. The thermochemical non-equilibrium scale effects are studied by changing the TR r t . Cases 3, 5, and 6 in Table 2 are selected, and the translational temperature T and vibrational temperature Tv are shown in Fig.9. In the downstream of the throat, the vibrational temperature and translational temperature are separated. The separation point of the nozzle is different. With the increase of throat radius, the separation point and freezing point gradually move to the downstream of the nozzle. For the throat radius r t are 10mm, 20mm, and 30mm, freezing temperatures are 2442K, 2207K, and 2140K.
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) As the TR increases, the pressure gradient near the throat decreases (see Fig.9), which raises the pressure downstream of the throat, resulting in decreased TVERT, as deduced from Equation (28). Additionally, as the TR increases, the velocity near the throat decreases and the flow time increases. The overall effect is that the freezing point is far from the throat and the flow is close to the equilibrium flow.
(2) For the equilibrium flow, the flow parameters on any cross-section are determined only by the nozzle stagnation conditions and the local area ratio. The throat x=0mm and the upstream of the throat is the equilibrium flow, and the translational temperature and the vibrational temperature are not separated, as shown in Fig.10.
(3) The increase in the kinetic energy of the equilibrium flow comes from the translational kinetic energy and vibrational energy. Since the vibrational energy of the non-equilibrium flow freezes, the increase of kinetic energy can only come from the translational energy. As the TR increases, the velocity u at the throat decreases. The velocity u of the nozzle outlet is only related to the total enthalpy h of the nozzle stagnation gas, and the nozzle outlet velocity is almost the same (see Fig.11). The contraction section and the throat are upstream of the freezing point of the high enthalpy nozzle. We change the CA θ c and TCR r c to study the thermochemical non-equilibrium effect. In the subsonic region upstream of the throat, the flow is the equilibrium flow. Changing the nozzle geometry in the equilibrium flow area has almost no effect on the non-equilibrium effect of the flow field. Cases 3, 7, and 8 in Table 2 are selected to study the influence of the CA on the nonequilibrium effects. The freezing position and vibrational temperature of the flow field are almost unchanged (see Fig.12). The change of the TCR has little effect on the non-equilibrium effect, and the freezing position and vibrational temperature of the flow field are almost unchanged. The calculation results of cases 3, 9, and 10 are shown in Fig.13. As the airflow passes through the throat, the non-equilibrium effects begin to increase, and the flow freezes quickly downstream of the throat.

Mach number
Although the high enthalpy nozzle is mainly used to simulate the speed and pressure of the flying environment, the Mach number is also a key parameter. It has been mentioned in the previous research that changes in the MEA and TR have an effect on the non-equilibrium effect. Now we study the effects of changes in the MEA and TR on the Mach number.
To study the influence of the MEA on the Mach number, cases 1, 2, 3, and 4 in Table 2 are selected. The centerline Mach number distribution is shown in Fig.14.
With the increase of the MEA, the Mach number gradually decreases. This is because the translational temperature increases, the local sound speed gradually increases (see Fig.15). Equation (29) shows that the sound speed has some dependence on the pressure gradient in the stream direction.
Cases 3, 5 and 6 in Table 2 are used to study the effect of the variation of the TR on the Mach number (see Fig.16). The variation of the TR changes the expansion ratio of the nozzle. The change of the expansion ratio changes the outlet Mach number, which is similar to the calorically perfect gas. When the throat radius becomes larger, the freezing point is farther from the throat, and the flow is closer 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 Equation (10). The change in Mach number due to changes in expansion ratio and non-equilibrium effects is much more pronounced than the change in Mach number caused by changes in expansion ratio alone.

Species mass concentration
Cases 1, 2, 3, and 4 in Table 2 are selected to study the influence of the MEA on the species mass fraction. Fig.17 shows the centerline evolution of NO mass fraction. It can be seen that the flow freezes not far downstream of the throat. The MEA variations also bring changes in the dissociation of NO. Considering the influence of the MEA on vibrational temperature, the variation of NO mass fraction can reflect the effect of vibrational temperature on dissociation. As the MEA increases, the freezing point approaches the throat, and the NO mass fraction increases. The non-monotonic change of NO along the centerline upstream of the throat is caused by the non-monotonic change of NO with temperature. For example, in the equilibrium air composition at one atmospheric pressure, the temperature range of NO exists is 2000-6000K, and the corresponding peak temperature is about 3500K Error! Reference source not found. . Under the conditions of this paper, the pressure at the upstream of the throat is high, and the temperature range and peak value of NO are corre-sponding to the temperature rise. The peak value of NO is about 5000K. Fig.18 shows that the mass fraction distribution of O 2 and O is also significantly affected by the different MEA. The remaining flowfield parameters are reported in Table 7. As the MEA increases, the flow field parameters tend to be non-equilibrium. The mass fractions of dissociating species such as NO, N, and O increase, while the mass fractions of N 2 and O 2 decrease. The influence of the different TR on the flow field parameters is also shown in Table 8. As the TR increases, the nozzle flow field parameters are biased to the equilibrium state.   Table 8 The flowfield parameters in the centerline at nozzle outlet with the different throat radii rt Case r t (mm)

Conclusions
As this paper has demonstrated, the high enthalpy nozzle has the thermochemical non-equilibrium expansion effects, which has appreciable influence on the flowfield parameters of the nozzle outlet. The following conclusions can be drawn: (1) Under the same expansion ratio, as the MAE increases, the thermochemical non-equilibrium effects of the flow field become more obvious. As the TR of the nozzle increase, the flow approaches the equilibrium state. The TCR and CA have little effect on the non-equilibrium effect. The smaller MEA and the larger TR can effectively suppress thermochemical non-equilibrium effects.
(2) With the increase of the MEA, the vibrational temperature increases, but the translational temperature and the speed of sound decrease, which makes the Mach number increase. As the TR increases, the vibrational temperature decreases, but the translational temperature and the speed of sound increase, making the Mach number decrease.
(3) The thermochemical non-equilibrium effects change the species mass fraction. The larger the MEA or the smaller throat radius is, the larger the dissociation species O, N, NO and NO + at the nozzle outlet is, but the smaller the mass fraction of the species O 2 and · 11 · N 2 is.