Numerical simulation of hypersonic thermochemical nonequilibrium flows using nonlinear coupled constitutive relations

Shuhu ZENG, Zhenyu YUAN, Wenwen ZHAO,*, Weifng CHEN

a School of Aeronautics and Astronautics, Zhejiang University, Hangzhou 310027, China

b School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, China

KEYWORDS Hypersonic flow;Nonlinear coupled constitutive relations;Rarefied gas;Thermochemical nonequilibrium effect;Vibrational excitation

Abstract To predict aeroheating performance of hypersonic vehicles accurately in thermochemical nonequilibrium flows accompanied by rarefaction effect, a Nonlinear Coupled Constitutive Relations (NCCR) model coupled with Gupta’s chemical models and Park’s two-temperature model is firstly proposed in this paper.Three typical cases are intensively investigated for further validation, including hypersonic flows over a two-dimensional cylinder, a RAM-C II flight vehicle and a type HTV-2 flight vehicle.The results predicted by NCCR solution,such as heat flux coefficient and electron number densities, are in better agreement with those of direct simulation Monte Carlo or flight data than Navier-Stokes equations,especially in the extremely nonequilibrium regions,which indicates the potential of the newly-developed solution to capture both thermochemical and rarefied nonequilibrium effects.The comparisons between the present solver and NCCR model without a two-temperature model are also conducted to demonstrate the significance of vibrational energy source term in the accurate simulation of high-Mach flows.

1.Introduction

During the reentry process of hypersonic aircraft, the gas molecules in the upstream regions will greatly transfer their kinetic energy into internal energy1due to the strong compression and viscous effects inside shock waves.2Thus, a sharp increase in temperature occurs behind the shock wave and around the Thermal Protection System (TPS) of hypersonic vehicles.3This will result in a series of typical hightemperature gas effects, such as vibrational energy excitation,molecules dissociation and ionization,which means the calorimetric perfect gas model will lose its accuracy.For example,thermochemical nonequilibrium flows with a twotemperature model past a reentry Apollo Command Module were simulated by Hassan et al.,4and it showed that the trim angle and lift-to-drag ratio including real gas effects were less than those from wind tunnel experimental results with nonreacting gas, but they were in good agreement with the flight experimental data.The high-temperature gas effect is still one of the vital problems for the accurate prediction of the aerothermodynamics of hypersonic vehicles.However, it is quite difficult to reveal the complex fluid mechanisms in ground test,5and flight data are too expensive to obtain plentifully.Hence,numerical simulation has become an effective approach to simulate these high-temperature flows.

Till now, a variety of CFD solvers,6such as Langley Aerothermodynamic Upwind Relaxation Algorithm7(LAURA), viscous upwind algorithm for complex flow analysis8(VULCAN) and NNW-HyFLOW,9have been developed to solve the well-known Navier-Stokes(NS)equations coupled with chemical kinetic models or thermochemical nonequilibrium models.Meanwhile, much effort has been brought into the development and validation of different thermal and chemical models proposed by Blottner,10Park,11,12Gupta,13Dunn,14etc.All these reaction models are deeply applied for the description of hypersonic flows under the atmospheric conditions, but it is noteworthy that different models may yield significant differences in the thermochemical nonequilibrium processes and distributions of the gas species,as demonstrated in many existing literature.15–18In addition to the thermochemical nonequilibrium effect, another problem faced by the hypersonic aircraft at high altitudes is the rarefied nonequilibrium effect.The flows of perfect gas, multi-species chemically reacting gas and rarefied gas through an Apollo-style capsule were calculated by Pezzella and Votta19using the NS solvers named H3NS20and DSMC solver named DS2V.21The results underlined that coupling with the rarefaction and high-temperature effects had a great impact on both aerodynamic and aerothermodynamic characteristics of capsule.Furthermore, it revealed that NS equations based on small deviation from local thermodynamic equilibrium become questionable and inadequate in high-Knudsen-number flows where rarefied nonequilibrium effects become dominant.

The Direct Simulation Monte Carlo(DSMC)22–25is unified for the simulation of reaction flows from the continuum region to the rarefied regime.This is because DSMC method has a stronger physical background and can simulate the chemical reaction process at the molecular level.There are many chemical reaction models developed based on DSMC,among which the Total Collision Energy(TCE)model proposed by Bird21,26is the most popular one.The core success of the model is to calculate the chemical reaction probability between molecules through the total collision energy.However, the TCE model excessively relies on the chemical reaction rate coefficient obtained mainly from experimental data.For the hightemperature flows, the completeness and accuracy of data still need to be further corrected.In recent years, the Quantum Kinetic (QK) model26,27applied to DSMC method has attracted extensive attention of many researchers around the world.By directly establishing the chemical reaction model from the micro-scale level, the previous dependence on the chemical reaction rate is skillfully avoided.However, the efficiency of the existing DSMC method for the simulation of continuum and transition flows is relatively low28compared with the solution of partial differential equations.

To make up for the deficiency of NS equations in the transition regime and overcome the low computational efficiency of DSMC method, Eu proposed Generalized Hydrodynamic Equations29,30(GHE),which is strictly consistent with the second law of thermodynamics.In Eu’s theory, the connection between the dissipative evolution of macroscopic nonconserved variables and the entropy production is tactfully built up with the help of the exponential nonequilibrium distribution function,which is proposed by himself.However,GHE has been merely utilized to calculate the shock structure,31and it is very difficult to be extended to multi-dimensional problems because of the highly nonlinear coupled terms involved in the equations.Recently,it was noted that a shock structure singularity problem32,33would emerge similarly in GHE without introducing Eu’s adiabatic assumption.31Following Eu’s work, Myong implemented the adiabatic assumption and balance closure32on GHE.Then, he developed an efficient nonlinear algebraic system in the physical sense, namely,Nonlinear Coupled Constitute Relations34(NCCR).Subsequently, a decomposed algorithm35was established to solve NCCR model numerically, which has been employed successfully in one-dimensional shock wave structure and twodimensional blunt body36flow problems to validate the potential and capability of NCCR model in high-speed and rarefied nonequilibrium flows.More researches37–40have been carried out extensively for NCCR model in single-component gas,including a Discontinuous Galerkin(DG)method on unstructured grid,41,42an undecomposed algorithm43and a simplified analytical method44in the framework of Finite Volume Method (FVM) and a new wall boundary condition for micro-Couette flow,45and so on.Jiang et al.46proceeded to a linear stability analysis on multi-dimensional GHE and NCCR model in the equilibrium rest state as well as uniform-moving state and it was demonstrated that they were unconditionally stable for all wavenumbers and frequencies.The connection and discrepancy between Eu’s and Grad’s moment equations were also deeply clarified with the truncated distribution function in this research.To reveal the rotational nonequilibrium effect in the transition regime, a multitemperature NCCR model with a modified bulk-to-shear viscosity was derived by Yuan et al.47for the accurate simulation of diatomic rarefied gas flows.Moreover,it was found that the bulk viscosity ratio is significant to determine the topological aspects of diatomic and polyatomic gas flows based on NCCR model.48Mankod and Myong49firstly proposed a new NCCR model including the vibrational mode for single-component diatomic and polyatomic gases, and the classical shock structure problem was accordingly solved with special emphasis on the interaction between the vibrational-translational relaxation and the second-order effects in NCCR model.On the basis of GHE, Ahn and Kim48built up an axisymmetric Generalized Hydrodynamic (GH) model to investigate hypersonic chemical nonequilibrium problems, in which the nonlinear mass diffusion flux of mixed gas was considered.At the same time, five-species GH model including O2, N2, NO, O and N was applied to one-dimensional shock wave structure as well as hypersonic rarefied flows past a sphere and a reentry body.The results showed that multi-species GH model yielded closer results to the DSMC data than NS solvers.However, only two-dimensional axisymmetric problems were studied in this paper due to the computational cost of multi-species GH equations.Recently, the hypersonic reaction flows in the transition regime were investigated by Yuan et al.50using NCCR equations coupled with chemical reaction source terms, and the comparisons between NCCR solver with perfect gas and multi-species reacting gas were also conducted,which indicates the importance of chemical reaction source terms.

The goal we pursue in this paper is focused on the disclosure and analysis of flow mechanisms for high-speed flows coupled with the thermochemical nonequilibrium and rarefied gas effects based on the contribution of NCCR model.Thus, the hypersonic reacting flows over two typical configurations,RAM-C II and type HTV-2,are numerically investigated using Gupta’s finite rate chemical reaction models13with twotemperature model.51The rest of this paper is structured as follows.In the next section, a comprehensive mathematical description of the governing equations and nonlinear constitutive models for multi-species gas is provided.In Section 3, we mainly depict the details of the chemical kinetic model and transport properties, as well as the vibrational energy source term.Subsequently, the numerical results of the two typical test cases are compared and analyzed in detail in Section 4,whereas a hypersonic reacting flow around a cylinder with the existing DSMC results is employed as a validation case to assess the performance of the present numerical models before those.Conclusions are then presented in Section 5.

2.Governing equations and nonlinear coupled constitutive relations for multi-species gas

The high speed and high temperature of hypersonic reentry flow will cause extremely thermal nonequilibrium, vibrational and electronic excitation, and chemical nonequilibrium phenomena.The detailed depiction of these phenomena is not sufficiently modeled by perfect gas and requires the usage of a real gas model.In addition,the classical macroscopic NS equations based on the local thermodynamic equilibrium approximations may not be an ideal approach for the accurate description of hypersonic rarefied gas flows.To accurately calculate the farfrom-equilibrium flow problems coupled with the thermochemical nonequilibrium effects and obtain accurately engineering data to support the designs of the flight trajectory and the TPS for hypersonic aircraft, this paper is committed to extending the Nonlinear Coupled Constitutive Relations(NCCR) for a single gas to multi-species gases with Park’s two-temperature model.Thus, the effects of hypersonic and high-temperature chemically reacting processes can be taken into account felicitously.

2.1.Governing equations

For the multi-species formulation of the conservation laws of species mass, mixture momentum, total energy and total vibrational-electronic energy are given as follows:

The mass conservation for each species i is governed by.

in which I means the unit second-rank tensor;Π represents the shear stress;Δ is associated with the bulk viscosity and denotes the excess normal stress.The work conducted by Kustova et al.52shows considering the bulk viscosity effects in the description of hypersonic flows can result in more accurate shock layer around the aircraft and the heat transfer to the surface of the vehicles.The mixture pressure p is given as the sum of partial pressure of each species and evaluated by Dolton’s law and the perfect gas law for each species,53

where R is the universal gas constant; Mimeans the molar mass of species i; T denotes the translational-rotational temperature and is used to describe the distribution of translational and rotational energies under the two-temperature model assumption of Park.Furthermore, it is assumed that the vibrational,electronic excitation and electron translational energies are described by the other temperature named vibrational-electronic temperature Tvib.This is a reliable assumption which is based on the fact that the rotational energy equilibrates with the translational energy in just several collisions and the energy transfer between the translational mode of electrons and the vibrational mode of molecules is very fast in air.54

The total energy conversion equation is.

where evibis the total vibrational energy per unit mass of mixture, eivibrepresents the species vibrational energy, and ˙ωvibdenotes the vibrational energy source term.

Together with the conservation Eqs.(1), (2), (4) and (5),and as mentioned above, the governing equations of the multi-species can be summarized as.

It is observably noted that Eq.(6)is still a multidimensional open system of partial differential equations without the help of these constitutive relations, which may directly determine or control the accuracy of hydrodynamic equations.Hence,we introduce the nonlinear constitutive models to close the mathematical system in the following for the high-Knudsennumber and hypersonic flows where nonequilibrium effects become significant and dominant.

2.2.Nonlinear coupled constitutive relations

The NCCR model has been intensively studied in plentiful work,55–58including the computational method and the boundary conditions.Meanwhile, it also has recently gained great success in engineering applications57and stable numerical predictions of hypersonic rarefied flow and even chemical reaction flow over three-dimensional complicated aircraft under various flight conditions.50However, there is an urgent demand for the NCCR theory to consider vibrational excitations in chemically reacting gases.For completeness and computational efficiency, we still implement multi-species NCCR model in conjunction with the linear vibrational heat flux Qvib= λvib?Tvibsimilar to the form of the NS equations,which is not only theoretical but also logical.Here, λvibis vibrational-electron-electronic thermal conductivity.The NCCR model for multi-species gas is derived and simplified from the Generalized Hydrodynamic Equations (GHE) proposed by Eu30in 1992.Different from the Chapman-Enskog(CE) method59expanding the velocity distribution function with small Knudsen number, the development of GHE theory was based on the Boltzmann-Curtiss kinetic equation60and started from a nonequilibrium distribution function in exponential format compared with maximum entropy distribution.61Another point is noteworthy that the collision term was established by the employment of the cumulantexpansion method.62,63The exponential distribution function fican be expressed as.

where nistands for the number density of i species, ^hiis the enthalpy density per unit mass, and the bracket symbol [A](2)means the traceless symmetric part of the second-rank tensor A.

If one wants to get a more detailed derivation and description of the nonlinear dissipation term with great emphasis on Eu’s primary theory, uncovering the previous investigations conducted by Eu,29,63Myong,32Xiao,64etc.should be a pleasant way to do.Using the closing-last balanced closure32or Eu’s closure,31adiabatic approximation31and Myong’s simplification,34,65a set of nonlinear algebraic equations for multi-species flows,namely the NCCR model, is developed as follows:

In Eq.(12), the shear stress, excess normal stress and heat flux of each species are modeled individually.If the above nonlinear algebraic system is applied to describe the high-speed flows without any simplification, the expensive computational consumption and great storage are so difficult to support in terms of the current limited resource.For the sake of calculational efficiency and accuracy, the concept of the mixed gas is adopted to consider each species uniformly.Thus, the NCCR model can be reduced to.

2.3.Numerical methods to solve the governing equations and constitutive relations

Based on the framework of the finite volume method,the governing equations in partial differential form and the NCCR model in an algebraic form are attempted to be solved numerically,which has been validated in many published literatures.50,65,57,58The Advection Upstream Splitting Method by pressure-based Weight function66(AUSMPW+)is extensively applied in the simulation of hypersonic flows due to its robustness and strong capability in capturing the shock and is also used in computing the convective flux of governing equations in this work.Meanwhile, a highresolution Van Albada limiter with less dissipation is adopted together with Monotonic Upstream-centred Scheme for Conservation Laws(MUSCL)to reconstruct the left and right values of the interface.The spatial derivatives of viscous flux are discretized by the second-order central difference scheme.To obtain an efficient steady convergence, the implicit Lower-Upper Symmetric Gauss-Seidel relaxation (LU-SGS) is employed for the temporal discretization.

Note that the only difference between the solution of traditional NS equations and the NCCR model lies in the treatment of viscous stresses and heat flux appearing in governing equations.These non-conserved variables can be computed explicitly in NS Eq.(15),but implicitly in the NCCR model(13)due to the existence of the nonlinear factor q(κ ).Therefore, an additional iterative process is necessary to numerically solve the NCCR model.This paper focuses on the solution of the NCCR model using an undecomposed solver developed by Jiang et al.,43where the fixed-point iterative and Newton iterative formulations are coupled together to obtain a steady convergence.More details about the undecomposed algorithm and its convergent condition can be found in Ref.43.

3.Thermochemical nonequilibrium models

3.1.Chemical process and vibrational energy source term

A chemical kinetic model is a crucial factor of numerical simulation for hypersonic flows coupled with thermal nonequilibrium phenomena and chemical reacting processes.In this paper, Gupta’s chemical reaction model is adopted because it has been widely utilized in the literature of hypersonic flows simulations.15–17There are 20 reaction equations for eleven species in the Gupta model,13where the first six reaction formulas and the first seven reaction formulas are used for five species and seven species, respectively.A general finite-rate reaction model can be expressed as.

where Af,rand Ab,rrepresent the frequency factor for forward and backward reaction,respectively;Bf,rand Bb,rare named as the temperature exponent;Ef,rand Eb,rdenote the ratio of activation energy to the universal gas constant.More details about these constants appearing in the above Arrhenius expressions are given in Zhao’s dissertation.68Tcdenotes the controlling temperature of chemical reaction and can be obtained under Park’s two-temperature model by

where a and b=1-a are the exponential weight factors of translational-rotational and vibrational-electronic temperatures, respectively.The weight factors of the forward and backward rates in detail are arranged in Zhao’s dissertation.68Due to the endeavor of these controlling temperatures, the influences of thermal nonequilibrium on chemical reactions can be considered and described at large.

According to the law of mass reaction,the mass rate of production of species i per unit volume ˙ωiis computed by.53

where τimeans the vibrational relaxation time,and is modeled by Millikan-White’s semi-empirical expression69and Park’s correction.11

To simplify the calculation, the simplified form of translational-vibrational energy relaxation term is employed in this work according to Ref.70.

3.2.Transport properties

The transport properties of gases include viscosity coefficient,thermal conductivity and diffusion coefficient, and are generally calculated using the fitting formula.For singlecomponent gas, the viscosity coefficient can usually be obtained by the inverse power law or Sutherland viscosity model, while the coefficient of heat conduction is given by the Prandtl number relationship.Furthermore, Wilke’s semiempirical mixing rule71needs to be used to approximate the mixture transport coefficient of mixed gas.The specific expressions are written as.

where T represents the local temperature.Ai,Biand Ciare constants determined for viscosity curve fit of each species,and their values for an eleven-species air model are listed in Table 1.

Different from the perfect gas, there are both monatomic gases and diatomic gases in the mixed gas, and excess normal stress does not exist in the monatomic gases.Meanwhile, the viscosity coefficient ratio fb=ηb/η of different diatomic gases is different.Therefore, the calculation of the bulk viscosity ηbfor multi-species gases needs to be corrected on the basis of single-component gas.Referring to the modified method in Ref.73, the viscosity coefficient ratio is determined by the weighted average according to the mole fraction of diatomic gas in the mixed gas as follows:

Table 1 Constants for viscosity curve fit in Gupta’s model.

4.Results and discussion

In this section, three numerical cases are selected to test the capability of NCCR model in conjunction with thermochemical nonequilibrium models:(A)a half-cylinder without a wake at Mach number 27.8;(B)the RAM-C II flight vehicle at Mach number 23.9 and 28.3; (C) the type HTV-2 flight vehicle at Mach number 20.Among these representative cases, case (A)is employed to validate the accuracy and physical consistency of the numerical models conducted in this paper.The others are applied to assess the performance of NCCR model with the Gupta chemical process and two-temperature model, and then analyze the effect of thermal nonequilibrium in hypersonic reaction flows.For the two-temperature model, it is assumed that a single temperature controls the translationalrotational modes,and the vibrational-electronic mode of molecules is described by a single temperature.Both flow field and surface properties obtained from the above numerical simulations are analyzed and compared with DSMC or NS results to reveal more details of the thermochemical nature of high-speed flows.

4.1.Code validation

The hypersonic reacting flows around a half-cylinder without a wake at Mach number 27.8 are simulated by NCCR model with a coupled algorithm43in the transition regime.The radius of the cylinder is 25.4 mm.To validate the present numerical solver against the existing DSMC results reported by Fan and Shen,75Gupta’s five-species chemical kinetic model with thermal-vibrational nonequilibrium is used in this case.The computational conditions are listed in Table 2.The wall boundary condition is assumed at a constant temperature of 1000 K, and the interaction between gas flow and solid wall is the first-order Maxwell-Smoluchowski (M/S) slip and jump boundary condition, which has been widely used and is also applied in the next two cases.The inflow gas contains O2and N2, with a mass fraction of 76.3 % and 23.7 %,respectively.

The hexagonal structured mesh in the symmetrical plane is illustrated in Fig.1 with 81 and 59 points set in circumferential and radial directions of the cylinder, respectively.The grid space in the near-wall region is very important to achieve heating convergence and reliable heating prediction.Therefore, a grid independence study for this numerical case is carried out as an example with three sets of grids, in which the first grid spacing in the wall normal direction is 1 × 10-4m (Mesh 1),1 × 10-5m (Mesh 2) and 5 × 10-6m (Mesh 3), respectively.The refined mesh 2 and mesh 3 yield almost the same solution,whereas the coarse mesh 1 shows obvious deviation in Uvelocity distributions along the stagnation line, as illustrated in Fig.2.By taking efficiency into consideration, we selected refined mesh 2 for further investigation.Furthermore, the mesh 2 refined with a normal spacing of 1 × 10-5m in wallnormal directions can rightly ensure that the cell Reynolds number meets the accurate numerical requirements.It is noted here that the demand for aero-thermodynamic convergence and reliable aero-thermodynamic prediction can be achieved with Recell? (0,1 ),which is verified by many numerical experiences50,76and obeyed by the two other cases.The cell Reynolds number is defined as.

where ρ∞,u∞,η∞are density, velocity, and shear viscosity of free stream, respectively, and Δx denotes the first grid spacing in the wall-normal directions.

Fig.3 shows the convergent history for the simulation of flows around a half-cylinder.The total iteration includes 150,000 time steps and the final residuals indicate that the calculated results have converged finely and the accuracy of results can be ensured.Likewise,the latter two numerical cases defer to the same convergence criterion.From Table 2, the Knudsen number of free streams equals 0.16, and the flow is identified as the transition flow, where the conventional NS equations may not provide an accurate prediction of flow details.As shown in Fig.4,the distribution curves of heat flux coefficient Chalong the solid surface calculated by the in-house NCCR solver and NS solver are compared, and at the same time,the DSMC results from Fan and Shen75are also included as a comparison.It is seen that the heat flux coefficient obtained by NCCR model is closer to DSMC results than NS equations.This is the strong support for the fact that NCCR model can obtain more accurate results than NS equations in transition flows with thermochemical nonequilibrium.On the other hand,we can draw a preliminary conclusion that the NCCR solver in conjunction with chemical kinetic models and Park’s two-temperature model developed in this paper is feasible to capture properly the thermochemical nonequilibrium in hypersonic flows.The contours of Mach number,pressure, translational temperature (T), and vibrational temperature (Tvib) are given in Fig.5.One can observe that there is a distinct discrepancy between the flow field properties predicted by NCCR model and NS equations.NCCR model not only obtains larger shock wave thickness before the stagnation point, which is a representative phenomenon induced by the rarefaction gas effect,but also shows significantly lower temperature peak after shock wave than NS results, for both translational temperature and vibrational temperature.Fig.6 compares the temperature profiles along the stagnation line.It is found that there are strong nonequilibrium phenomena in the shock as well as in the boundary layer.Meanwhile,the value of translational temperature and vibrational temperature exhibit a large difference as depicted in Fig.6, so that taking the thermal-vibrational nonequilibrium into consideration in high-temperature gas flows is significant to evaluate accurately the physical properties of nonequilibrium hypersonic flow.

Fig.7 presents the comparison of Molar fraction contours of all five species evaluated by the two solutions.Clearly, the numerical results illustrate that chemical reactions take place intensively behind the bow shock accounting for the booming temperature.It is noteworthy that the dissociations of N2and O2are extensively weakened in NCCR model compared with NS solvers due to the difference of temperature illustrated in Fig.5 and Fig.6.More details about the Molar fraction of different species along the stagnation line are plotted in Fig.8.Because the controlling temperatures of chemical reactions will directly affect the chemical reaction process, it can be clearly observed from Fig.8 that the mole fraction distributions of NO,N and O simulated by NCCR model are much lower than those of NS solvers.In other words, the chemical nonequilibrium effect captured by NCCR model is significantly weakerthan that of NS equations under the flow conditions depicted in Table 2 where extremely rarefied nonequilibrium effect exists.The rarefaction effect causes less collision between the molecules and a longer time is taken to reach the equilibrium state.All the above comparisons suggest that NS equations may overrate the aerothermal characteristics in the hypersonic flows coupled with the thermochemical nonequilibrium and rarefied nonequilibrium effect.However, the present NCCR model possesses the potential and capacity to describe properly the nonequilibrium nature and flow chemistry to some extent.

Table 2 Free stream conditions for cylinder.

Fig.1 Computational mesh for cylinder.

Fig.3 Convergence curve of residual for simulation of a halfcylinder.

Fig.4 Comparison of heat flux coefficients on cylinder surface.

4.2.RAM-C II vehicle at Mach number 23.9 and 28.3

Fig.5 Comparison of contour between NCCR model and NS equations.

Fig.6 Comparison of temperature along stagnation line.

For further assessing the performance of the present numerical model before employing it in the engineering problems of actual vehicles, a hypersonic thermochemical nonequilibrium flow past a reentry blunt cone at two flight altitudes of 61 km and 81 km is introduced.The blunt cone dates from the end of the 1960 s when a series of the Radio-wave Attenuation Measurement(RAM)experiments concerned with the reentry process of hypersonic vehicles were conducted to study the communications blackout problems.Among them, there are detailed experimental results obtained in the second flight test,77including the electron number density distributions near or on the surface of the blunt cone named as a RAM-C II vehicle.This case has been comprehensively investigated in many pieces of literature.16,68,78As depicted in Fig.9,the head radius(R) of RAM-C II vehicle is equal to 0.1524 m, the half cone angle (α/2) and overall length (L) are 9° and 1.295 m, respectively.The calculation conditions for the two altitudes are listed in Table 3.The completely non-catalytic and isothermal wall boundary conditions are extensively used in the simulation of high-temperature and hypersonic gas flows, as well as in this paper.The mass of O2accounts for 79 percent in the inflow gas, and the rest is N2.Park’s two-temperature model and Gupta’s eleven-species reaction model are employed in this subsection.

Fig.8 Molar fraction of different species along stagnation line.

Fig.10 plots translational temperature and vibrational temperature at an altitude of 61 km simulated by NS solver and NCCR model with thermochemical nonequilibrium models.The more detailed comparisons about the translational temperature, vibrational temperature, Mach number and density along the stagnation line are illustrated in Fig.11.These figures suggest that the prediction results of the two solutions considering the vibration nonequilibrium basically agree well with each other at the height of 61 km where the flow meets qualitatively the continuum assumption, and only some slight discrepancies appear in the shock wave layer while the peak value of translational temperature obtained by NS equations has a small lead over NCCR solver.We need to emphasize that the simulation of the NCCR model without the vibrational excitation is also displayed with the red dashed line in Fig.11 for the analysis of the thermal nonequilibrium effect.It is well known that the vibrational energy excitation process extensively exists as far as the atmospheric air when the temperature is over 800 K.This particular comparison,we expect,would indirectly help us realize the fundamentality of the vibrational nonequilibrium model in the accurate description of hypersonic flow and design of TPS for a high-velocity vehicle.

Fig.9 Geometry of RAM-C Ⅱvehicle.

As shown in Fig.12, the contours of translational temperature and vibrational temperature at 81 km have a significant difference.Especially, it is obvious that a larger shock wave stand-off distance and stronger rarefied nonequilibrium effect are captured by NCCR solution than NS solver.From Table 3,the Knudsen number of the free stream is 0.0328, where the rarefied nonequilibrium effect is extremely pronounced and the ability of NS equations to predict nonequilibrium phenomena is limited.Nevertheless, the NCCR model can give full play to its advantages in this case, which has been demonstrated in many studies.55–58Fig.13 gives the distributions of the temperature, Mach number and density along the stagnation line at 81 km.From the temperature distribution curves along the stagnation line presented in Fig.13(a), it is found that the temperature peak predicted by NCCR model is lower than that of NS equations when the vibrational energy relaxation source term is taken into account, but the vibrational temperature near the wall region is slightly higher than the results of NS solution.These results indicate that the vibrational nonequilibrium effect predicted by NCCR solver is a little weaker than NS equations in a sense.At the same time,we have to admit that the linear vibrational heat flux still conducted in this work is not appreciably matched with the nonlinear high-order system in NCCR model.For completeness,the nonlinear vibrational heat flux will be developed and employed in the simulation of hypersonic thermochemical nonequilibrium flows in the next step.Additionally, note that the temperature computed by NCCR model just with chemical reactions is also presented with a red dashed line in Fig.13(a).It should be emphasized the translational temperature simulated by NCCR solver with vibrational excitation is significantly higher than the temperature where a one-temperature model is adopted.This is because a great deal of vibrational excitation energy of molecular gas is released and converted into translational kinetic energy to make up for the expense resulting from the chemical reaction process.Interestingly,comparing the flow properties at the altitudes of 61 km and 81 km as plotted in Fig.11 and Fig.13, we can see that both the rarefaction gas effect and thermal-vibrational equilibrium effect gradually enhance as the flight height and velocity of vehicles increase.

To further estimate the performance of NCCR solver with a two-temperature model in hypersonic flows in conjunction with the rarefied nonequilibrium and thermochemical nonequilibrium effect,Fig.14 shows the electron number densities simulated by NS, NCCR, Simplified Conventional Burnett(SCB),68DSMC78and flight data77measured by a series of reflectometers at R1,R2,R3,R4and L1at 81 km.The five locations (R1-R4and L1) are 4.47 cm, 23.16 cm, 70.05 cm,106.04 cm and 123.42 cm from the top point of RAM-C II vehicle.The simulations of NCCR algorithm considering the thermal-vibrational nonequilibrium show a better quantitative and qualitative agreement with the experimental data or DSMC results than the rest solutions, while the slight difference appearing near the wall region suggests eagerly us to develop proper wall boundary conditions that fulfill the physical consistency of real flows in due course.Most importantly,if the vibrational excitation is not taken into account in NCCR solution, the results depicted in Fig.15 are much higher than those of other numerical models as well as the flight data,which indicates the necessity of the vibrational energy source term for expositing the extremely nonequilibrium nature of high-temperature and high-speed flows using NCCR model.The surface coefficients of pressure, friction and heat flux calculated by three solutions at 81 km are compared with each other,as seen in Fig.15.It is obviously found that the pressure coefficients obtained by NCCR solver,whether the vibrational excitation is considered or not, are relatively lower than NS equations.However, the other surface properties given by NCCR vibrational model are larger than NS solution near the head of the blunt cone.When not considering thethermal-vibrational nonequilibrium effect, the surface friction and heat flux coefficients are slightly lower than those of NCCR vibrational model.

Table 3 Computational conditions for hypersonic flow over RAM-C Ⅱvehicle.

Fig.10 Contours of translational temperature and vibrational temperature at 61 km.

Fig.11 Comparison of temperature, Mach number and density profiles along stagnation line at 61 km.

4.3.Type HTV-2 flight vehicle at Mach number 20

The last case investigated is the type HTV-2 flight vehicle at an altitude of 80 km based on NCCR model and NS equations considering the vibrational nonequilibrium effect.The performance of NCCR solution to simulate three-dimensional complex engineering configuration in a slip-transition region with extremely thermochemical nonequilibrium properties is preliminarily assessed in this part.The basic geometry of the type HTV-2 with an aerodynamic layout of wing-body-fusion is presented in Fig.16,which is firstly proposed in Xia’s dissertation.79Note that the aerodynamic characteristics of the waverider shape depicted in Fig.16 have been optimized on the foundation of the disclosed HTV-2.The Mach number of incoming flows is 20, as seen in Table 4.The wall boundary is assumed to be a non-catalytic and isothermal wall with a constant temperature of 1000 K.The M/S/slip boundary conditions with full tangential momentum and thermal accommodation coefficients are also applied at the solid surface.The working medium is air containing 79 % O2and 21 %N2.The two-temperature model and Gupta’s seven-species reaction model are used in the simulation.

Fig.12 Contours of translational temperature and vibrational temperature at 81 km.

Fig.13 Comparison of temperature, Mach number and density profiles along stagnation line at 81 km.

Fig.14 Comparison of electron number density at 81 km.

Fig.15 Comparison of surface properties at 81 km.

Fig.16 Schematic of type HTV-2 vehicle.

Fig.17 Temperature along stagnation line of type HTV-2.

The translational temperature and vibrational temperature along the stagnation line of the type HTV-2 flight vehicle obtained by two solutions show a remarkable difference,which suggests there are intense thermal-vibrational nonequilibrium effects behind the strong shock wave, as plotted in Fig.17.Meanwhile, a considerable local rarefaction phenomenon is found especially in the shock wave region of hypersonic vehicles, while the flaw of NS solution is commendably remedied by NCCR model in a manner.Fig.18 presents the surface properties simulated by NS equations and NCCR solver on the symmetrical plane (Z = 0) at 80 km.It is markedly seen that the results of NCCR algorithm are higher than those of NS solver in the sharp leading edge of type HTV-2.Subsequently, to give a global impression of the flow field of type HTV-2 predicted by NS and NCCR model more clearly,Fig.19 shows the Mach number, translational temperature and vibrational temperature contours of three sections(including X =0.5 m,1.2 m,1.98 m)at a height of 80 km.It is worth mentioning that an obvious bow shock wave is formed around the aircraft and the shock layer captured by NCCR solution is thicker than that by NS at 80 km,suggesting a strong rarefaction phenomenon at 80 km.To provide another perspective on Fig.17 and Fig.19.It can be discovered that the hightemperature area is dominantly concentrated on the head of high-velocity aircraft while the distributions of translational temperature and vibrational temperature exhibit a significant difference.The peaks of translational temperature are extensively bunched in the shock layer while the peaks of vibrational temperatures are mainly concentrated near the surface, so the vibrational excitation of molecules gives full play to its influence on the prediction of aerodynamic heating of hypersonic vehicles.

For further quantitatively comparing the results predicted by NS equations and NCCR model at 80 km,the distributions of physical variables of interest at the intersection lines of three different cross-sections(x=0.5 m,1.2 m,1.98 m)and symmetrical plane (Z = 0)along the Y-direction are given in Fig.20-Fig.22.As depicted in Fig.20, the Mach number curves near the wall surface computed by NS solver are basically consistentwith those of NCCR solution,while the difference between the two solutions is prominent in the area far away from the surface.By comparing the shock wave stand-off distance at different cross-sections, it is found that the shock wave structure simulated by NCCR algorithm is weaker than that of NS equations, which indicates the existence of rarefaction nonequilibrium effect and the potential of NCCR solver in nonequilibrium regimes.Fig.21 compares the translational temperature curves at different cross-sectional positions.It is found that the peaks of translational temperature predicted by NCCR solutions are lower than NS results.Unlike the translational temperature, the peak value of the vibrational temperature only appears near the surface, as depicted in Fig.22.From a comprehensive analysis in Figs.20–22,the following point should be emphasized here:the results of NCCR model and NS solutions are relatively similar in the zone close to the head of the type HTV-2 flight vehicle, which is given with red dashed line and blue solid line, respectively; in contrast, the difference between these two methods is larger in the tail regions.These results show that the extremely nonequilibrium effect extensively exists in the sharp leading edge and rear edge of type HTV-2 flight vehicle.To keep up with the higher requirements of powered high-speed vehicles, the present NCCR solution will be employed for the simulations of a three-dimensional hypersonic air-breathing cruise vehicle with wave-ride airframe for a challenge in the future work80,81.

Table 4 Computational conditions for hypersonic flow past type HTV-2.

Fig.18 Comparison of surface properties between NS equations and NCCR model on symmetrical plane (Z = 0) at 80 km.

Fig.19 Comparison of contours between NS equations and NCCR model at 80 km.

Fig.20 Mach number curves at different cross-sections.

Fig.21 Translational temperature curves at different crosssections.

Fig.22 Vibrational temperature curves at different crosssections.

5.Conclusions

This paper mainly studies the coupling effect of thermochemical nonequilibrium and rarefaction nonequilibrium in hypersonic flows by nonlinear coupled constitutive relations combined with Gupta’s chemical kinetic model and Park’s two-temperature model.To further validate the efficiency and accuracy of the proposed solutions about capturing the nonequilibrium properties, three typical numerical experiments including hypersonic flows past a cylinder, a RAM-C II flight vehicle and a type HTV-2 flight vehicle are carried out.From the numerical results obtained, the following conclusions can be drawn:

(1) Through comparing and analyzing the simulations at high altitude obtained by NCCR model and NS equations with thermochemical nonequilibrium models, the results predicted by NCCR solver, containing the heat flux coefficient and electron number densities,are in better agreement with DSMC results or flight data than NS solution in extremely nonequilibrium regions.This indicates the potential and physical consistency of NCCR model to capture nonequilibrium nature of hightemperature and hypersonic flows in a sense.

(2) When considering the vibrational excitation of molecules in NCCR solver, the translational temperature,surface friction and heat flux coefficients are higher than those of NCCR solution merely with chemical reaction processes.Most importantly, electron number densities near the wall surface of vehicles predicted by NCCR algorithm with two-temperature model yield a good agreement with DSMC as well as the flight data, while NCCR solution with the one-temperature model does not.Therefore, the thermal-vibrational nonequilibrium effect makes a big difference in the accurate description of high Mach number flow problems.

(3) To overcome the low computational efficiency of the DSMC method in hypersonic transition regime, an improvement over the conventional solution of NS equations is presented, with a nonlinear closure applied to replace the Stokes and Fourier linear relations for the stress tensor and heat flux,respectively.Due to the additional iteration process required for the solution of nonconserved variables, it can be concluded from the study that the computational time of the NCCR model is 1.5 to 2 times as traditional NS solver does.In other words,the NCCR solver is almost as efficient as the NS equations compared with the particle method, and more accurate with the results of DSMC and experimental data than the NS solver in nonequilibrium flows.

In conclusion, NCCR solver coupled with thermochemical nonequilibrium models is a relatively complete and effective approach to predict the coupling effect between rarefied gas and real gas in hypersonic flows.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

This work was financially co-supported by the National Natural Science Foundation of China (Nos.12002306, U20B2007,11572284 and 6162790014) and National Numerical Wind Tunnel Project, China (No.NNW2019ZT3-A08).Furthermore, the authors are grateful to Dr.Jiang Zhongzheng and Dr.He Zhiqiang for their discussions.