Numerical simulation of the thermal non-equilibrium flow-field characteristics of a hypersonic Apollo-like vehicle

Minghao Yu(喻明浩), Zeyang Qiu(邱澤洋), Bo Lv(呂博), and Zhe Wang(王哲)

Faculty of Mechanical and Precision Instrument Engineering,Xi’an University of Technology,Xi’an 710048,China

Keywords: atmospheric reentry demonstrator, reentry altitude, flow-field characteristics, two-temperature model

1. Introduction

During the process of reentry of human spacecraft, intercontinental missiles and returning satellites, strong shock waves will be formed at the front of the vehicle due to the very fast reentry speed. The gas in the shock wave layer is rapidly compressed and heated,resulting in a temperature increase in the flow-field around the vehicle, which causes chemical reactions such as dissociation, compounding and ionization of the gas components in the flow-field,forming a plasma sheath layer with a temperature of several thousand or even tens of thousands of degrees Celsius. This will lead to rapid loss of speech signals,telemetry signals and GPS positioning signals,and the vehicle will lose communication with the ground control center, causing a serious threat to the safety of the astronauts and the accurate landing of the vehicle. Thus, it is important to understand the variation of flow-field parameters around the reentry vehicle, and the aerodynamic shape of the vehicle has a significant impact on the flow-field around it.[1-3]

Current reentry vehicles are mainly divided into two types: sharp and blunt. The sharp type has a certain mitigation effect on the blackout phenomenon, reduces the thickness of the plasma sheath layer and increases the possibility of ground-to-space communication. However, in actual applications it is found that this structure not only increases the heat flow on the vehicle surface and causes severe ablation but also reduces the payload of the research target, which is not suitable for widespread use in future space transportation systems. The blunt-headed type makes up for the weaknesses of the sharp-headed type to a certain extent and has become the mainstream vehicle type in the aerospace field and deep space exploration missions today, such as China’s Shenzhou series of spacecraft, the US Apollo spacecraft and the European atmospheric reentry vehicle.[4]

Thanks to the structural advantages of the blunt-headed body, this type of vehicle is widely adopted and many researchers wordwide have conducted a lot of research on bluntheaded vehicles. Takahashiet al.[5]conducted numerical simulations of heat fluxes on the surface of a membrane reentry vehicle and the surrounding flow-field using a coupled multiphysics field solution and compared the results with the experimental data of Yamadaet al.[6]The results show that the numerical predictions are in agreement with the experimental data and that a membrane reentry vehicle can effectively mitigate aerodynamic heating during reentry. Jung[7]numerically investigated an axisymmetric and three-dimensional model of the OREX blunt-headed return vehicle by using a fourtemperature model considering the internal energy exchange and complex chemical reactions.He constructed a mathematical model that could accurately predict the complex flow-field in the shock layer based on the inlet conditions and reproduced the plasma flow around the OREX blunt-headed return vehicle. Tchuenet al.[8]used direct Monte Carlo simulation compared with solving the Navier-Stokes equation to study the chemical and thermochemical non-equilibrium of a twodimensional model of the SARA vehicle and gave the effect of thermochemical non-equilibrium on flow-field parameter variation. Papadopouloset al.[9]pointed out through simulations of NASA’s X-38 reentry vehicle that the magnitude of the vehicle’s surface temperature gradient is a key factor for predicting heat flow through the surface, and the grid Reynolds number can be used to calculate the normal height of the first grid layer of the surface. Men’shov and Nakamura[10]further found that the accuracy of heat flow simulation results can be ensured for grid Reynolds numbers less than 3. By studying the relationship between aerodynamic heat and mesh for obtuse cone and double obtuse cone models at high supersonic velocity conditions, Yanet al.[11]found that the magnitude of the heat flow is influenced by the height of the first grid layer of the surface,and the heat flow error can be effectively reduced when the number of normal grid points on the surface is more than 30. It was also pointed out that both the advection upstream splitting method (AUSM) family hybrid format and the Roe format can obtain more accurate surface heat flow distribution trends than other spatial format effects for both the axisymmetric and three-dimensional models, but the AUSM format is more accurate. Liu[12]studied the difference between the complete gas model,the single-temperature model and the two-temperature model by using an Apollo-like return vehicle. It was found that the surface pressure and friction results of the three models do not differ much, while the complete gas model has the highest calculated temperature.The chemical non-equilibrium effect of the two-temperature model is more obvious in terms of chemical reactions. Jia[13]investigated the relationship between the standing point region of a blunt-headed body and the aerodynamic thermal format by using experiments and simulations for comparison. She found that both shear-stress transport (SST)k-ωand realizablek-εsimulated standing point heat flow errors were less than 5%, which can better predict the aerodynamic thermal phenomenon of the flow.Zhou and Yan[14]investigated the effect of six improved spatial formats on thermochemically nonequilibrium flows by using a three-dimensional sphere model.The results showed that the improved spatial formats significantly improve the grid’s ability to capture shock waves, but there is little difference in the flow-field parameters for various cases. Lvet al.[15]investigated the effect of structural and nonstructural grids on the flow-field under specific operating conditions by using the double ellipsoidal model. It was found that both structures can respond to the flow-field parameter distribution but the structural grid can capture the shock boundary more clearly.

Most of these references use different types of threedimensional blunt-headed vehicle models for calculations,focusing mostly on popular research areas such as model accuracy, electromagnetic communication and aerodynamic heating, but they ignore the systematic analysis of basic properties such as the relationship between flow-field parameters and reentry altitude. The references that analyzed the flow-field parameters around the reentry vehicle are too old to match with the existing computational models.

Based on previous studies, this paper mainly takes a three-dimensional model of the European Space Agency’s Atmospheric Reentry Demonstrator(ARD)launched in 1998 as the research object by coupling the two-temperature model with the SSTk-ωturbulence model, and analyzes in detail the relationship between the surrounding flow-field including flow pressure, Mach number, temperature and other parameters and altitude when this type of vehicle is in the blackout region, so it can build the foundation for the subsequent systematic study of the blunt-headed vehicle by adding the chemical non-equilibrium model.

2. Physical models

The ARD is the European Space Agency’s first selfdeveloped Apollo-like test blunt-headed return vehicle designed to evaluate thermodynamic prediction models, collect data from the reentry flow-field,and help develop future space transport systems. The ARD is an important step towards Europe’s strategic goal of autonomous reentry, making Europe the third group,after the USA and Russia,to successfully master spacecraft launch and reentry technology.[16]The model of the ARD and its geometry are shown in Fig.1.

Fig.1. ARD model and geometry: (a)configuration of the ARD,(b)threedimensional model of the ARD.

The shape of the ARD is similar to a spherical cone and mainly consists of three parts: the front shield, the rear cone and the back cover. The front shield is a double-sided spherical structure with internal heat shields that can be subdivided into a head with a spherical radius of 3360 mm and a shoulder with an annular section radius of 140 mm. The rear cone has a half-angle of 33°and is used to mount the navigation system, attitude control system and other auxiliary electrical equipment. The back cover of the ARD is equipped with a signal receiver and a parachute system, which is responsible for both communicating with the ground base station during flight and ensuring the safety of landing of the ARD through the parachute system.

The liftoff and reentry sequence of the ARD is shown in Fig.2. In 1998,the Ariane 5 V503 rocket with the ARD was launched from the Guiana Space Center in Kourou, French Guiana. The rocket was launched on schedule,sent the target return vehicle to a predetermined orbital trajectory of 218 km above Earth and separated the ARD from the rocket by using the ejection mechanism to achieve free ballistic flight of the ARD in space. Then the ARD passed the highest point of its orbit and reentered the Earth’s atmosphere under the effect of microgravity. Eventually,the ARD achieved a soft landing in the Pacific Ocean through atmospheric friction and onboard parachute action.

Fig.2. Diagram of the ARD flight sequence.

During the ARD’s reentry process the blackout phenomenon appears at an altitude of 90 km-40 km above the ground, as shown by the black line in Fig.2. Since the flowfield parameters in the blackout area change drastically, the aerodynamic heating phenomenon is significant and the ARD may show obvious overload and vibration, so study of the flow-field in this area has great theoretical and practical significance. In this paper, we study the variation of the flowfield parameters of the ARD with an angle of attack of-20°at different flight altitudes in the blackout region. The known flow-field parameters and calculated conditions of the ARD are shown in Table 1.[7]

Table 1. Known flight parameters and calculated working conditions for the ARD.

3. Governing equations

3.1. Flow field

Since this study concerns the high-altitude flow-field characteristics of the ARD,we must pay attention to the continuity assumption of the Navier-Stokes equation and the second viscosity coefficient problem.

The Stokes assumption is still mainly used to simplify the computational process for studying shock waves.[17,18]Although the consideration of the second viscosity coefficient can describe the shock wave structure more accurately, at present we do not calculate the second viscosity coefficient because of its complexity and high computational costs.We may consider and study the second viscosity coefficient further in the future.

Whether the free flow is continuous or not in the highaltitude environment can be judged by Knudsen number(Kn).If theKnis less than 0.001, the flow-field can be assumed to be continuous.

The formula forKnis as follows:

wheredandnrepresent the diameter of the particle and number density,respectively.

To approximately evaluate the number density, the following equation for the gas state is used:

wherekis the Boltzmann constant(k ≈1.38×10-23J·K-1).

Since the free flow is in air,and N2is the main component of air,we use the particle diameter of N2to approximateKnof the high-altitude flow-field;the results are shown in Table 2.

Table 2. Knudsen number.

As can be seen from the table,the estimated values ofKnat 40 km-70 km are all less than 1.0×10-3, so the flow can be approximately deemed continuous.

3.2. SST k–ω model

The SSTk-ωturbulence model is a two-equation eddy viscosity model that adds an eddy viscosity limiting formulation to the baselinek-ωturbulence model.[19]The model combines thek-εmodel with the standardk-ωmodel by using blending functions: thek-εmodel is used to deal with the free-flow problem and the standardk-ωmodel is used to capture the viscous bottom flow. In particular,the specific expressions for the turbulence kinetic energykand the specific dissipation rateωare as follows:

whereΓkandΓωare the effective diffusivity ofkandω, respectively,Gkis the production ofkandGωis the generation ofω,YkandYωare the dissipation ofkandωcaused by turbulence,respectively,SkandSωare user-defined source terms,GbandGωbare thekandωcaused by buoyancy,respectively,andDωis the cross-diffusion term.

The specific expression of the eddy viscosity limiting formulation is as follows:

3.3. Two-temperature model

The return vehicle reenters the Earth’s atmosphere, and the high-temperature flow-field excites the internal gas molecular vibrations and electron energy modes, leading to the inability to explain the hypersonic flow phenomenon by the conventional Navier-Stokes equation. Instead, the twotemperature model.[20-22]is used to describe the thermal nonequilibrium effects of the gas during reentry.

The two-temperature model assumes that the energy is concentrated in two modes. One part is the translationalrotational mode,in which the rotational energy of the molecule is easily balanced with the translational energy of the particle in a very short time and can be described by a temperatureT.The other part is the vibrational-electron mode,where the vibrational energy of the molecule easily equilibrates with the electronic energy in a very short time,and is described by another temperatureTV. Since the relaxation times of the vibrational and electronic energies of molecules are usually very short, it is assumed that the temperatureTv, which describes the vibrational energy, is equal to the temperatureTe, which describes the electronic energy,i.e.,Te=Tv=TV. Therefore,for the two-temperature model, it is necessary to solve separately the equations for the conservation of translational energy (which also represents the equation for the conservation of total energy)and the conservation of vibrational energy(or the equation for the conservation of electronic energy)for the microscopic particles to characterize the thermodynamic state of the system.

The specific expressions of the two-temperature model are as follows:

where the specific expressions forQ,A,B,C,Av,Bv,Cv,andS,respectively,are

whereQandSare the conservation variables and source terms,respectively,A,B,C,andAv,Bv,Cvrepresent the inviscid and viscous fluxes in thex,y,zdirections, respectively,u,v,ware the velocities in thex,y,zdirections, respectively,His the total enthalpy,SVis the source term of the vibrational-electron mode energy equation, with specific expressions as follows:

wheree?v,sis the translational energy evaluated with the translational temperatureT,ev,sis the vibrational energy calculated with the vibrational temperatureTVandτsis the relaxation time between the translation and the vibration modes.

4. Grid setup

4.1. Height of the first layer grid

The height of the first layer grid of the surface has a large impact on accurate prediction of the heat flux across the ARD surface, so its height needs to be calculated. In turbulence simulation, the dimensionless surface distanceY+is usually used as the grid judging criterion for calculating the first layer grid height.[23]During the calculation, it is necessary to first make a prediction ofY+and use the predicted value to calculate the first layer grid heighty. After the simulation, we check whetherY+meets the requirements and, if not, adjust theyvalue and repeat the calculation. The expression for the calculation is as follows:

wherey,ρ,andμare the first layer grid height,free flow density and dynamic viscosity, respectively. In consideration of the SSTk-ωmodel used in this paper, it is taken asY+≈1.The skin friction velocityuτis calculated as follows:

研究中水齡最大的為25 d，這與其他水齡為10 d時(shí)的樣本中檢測(cè)出的各致病菌的數(shù)量基本在一個(gè)量級(jí)上，這可能是因?yàn)樗g為10 d時(shí)壓載艙內(nèi)的環(huán)境變化已趨于穩(wěn)定，致病菌的各生物活動(dòng)也趨于穩(wěn)定。此外，該研究對(duì)壓載艙新加裝的壓載水進(jìn)行15 d的跟蹤檢測(cè)發(fā)現(xiàn)，隨著水齡的增加，3種致病菌數(shù)量?jī)H在垂直分布上發(fā)生變化，且隨著水中其他動(dòng)植物的死亡，逐漸集中于壓載艙底部，最終呈現(xiàn)出隨水深的增加而逐漸增加，但在總數(shù)上并沒有像水生動(dòng)物和藻類一樣呈現(xiàn)出減少的趨勢(shì)。

whereτw,Cf,U∞,Re, andLare, respectively, the skin shear stress,skin friction coefficient,inflow velocity,Reynolds number, and feature length. Through repeated calculations of the above equations,the number of normal grid points on the surface and the height of the first grid layer of the ARD at four altitudes are obtained as shown in Table 3.

Table 3. The boundary layer parameters of the ARD grid at four altitudes.

4.2. Computational grid and convergence settings

The surface calculation grid of the three-dimensional model of the ARD at 50 km is given in Fig. 3. We take this height calculation grid as an example because the shoulder is the transition part connecting the head and the rear cone with a special geometric structure. At the same time,the shock wave mainly acts on the head of the ARD and there are intense variations in the flow-field around the two parts so the grid should be encrypted for these two parts. In the actual calculation,since the ARD is a left-right symmetric structure,in order to reduce the calculation and improve the calculation efficiency,the right half of the structure is taken for calculation under the condition of ensuring the height of the first layer grid.

The ARD reentry problem is part of a hypersonic calculation, and the density-based implicit solver is selected for compressible flow. At the same time, the energy equation is activated and the two-temperature model is checked to simulate the thermal non-equilibrium phenomenon at hypersonic speeds. The SSTk-ωturbulence model is selected with the option ‘compression effects’. Since the flow boundary is far away from the ARD,the inlet boundary is set as pressure-farfield,the outlet boundary is set as pressure-outlet and the angle of attack is-20°,while the surface of the ARD is set as a stationary wall with no slip at a constant temperature. The‘highspeed-numerics’ command is used to activate the high-speed numerical format, the AUSM numerical dissociation format is used and the item ‘convergence acceleration for stretched meshes’is checked to accelerate the convergence of the highly stretched and anisotropic mesh solutions. Finally, full multigrid initialization (FMG initialization) is performed to accelerate the convergence of the numerical solution by adjusting the Courant number.

Fig.3. Surface calculation grid for the ARD at 50 km.

5. Results and discussion

Figure 4 shows the distribution of the pressure around the ARD with an angle of attack of-20°at different flight altitudes. It consists of pressure clouds of the flow-field in the shock wave at altitudes of 40 km, 50 km, 60 km, and 70 km and pressure clouds of the surface of the ARD;the unit is Pa.Figure 5 shows surface pressure curves of the ARD at different altitudes. The horizontal axis represents the horizontal coordinate of ARD in the flow-field and the vertical axis represents the surface pressure. The surface pressure curve of the back cover of the ARD is shown in the partially enlarged view.

Fig.4. Flow-field pressure and surface pressure distribution around the ARD at different altitudes: (a)40 km,(b)50 km,(c)60 km,(d)70 km.

Fig.5. Surface pressure curves of the ARD at different altitudes.

During the reentry of the ARD with an angle of attack of-20°, according to the flow-field pressure distribution the flow-field high-pressure area is concentrated in the shock wave layer around the windward side of the head of the ARD,and the flow-field pressure around the rear cone and the back cover is much smaller than the flow-field pressure in the highpressure area. The peak pressure of the flow-field tends to increase and then decrease with decreasing altitude,and the peak pressure of the flow-field at 50 km is the maximum compared with the other three altitudes,which is about 21 kPa.

According to the surface pressure distribution, it can be obtained that the high-pressure area of the surface is concentrated on the windward side of the head. It is distributed in a toroidal shape, and the pressure value is inversely proportional to the radius of the torus. After passing the shoulder of the ARD,the surface pressure of both windward and leeward sides of the rear cone area decreases rapidly. Because the connection between the rear cone area and the back cover is designed as a step structure, the surface pressure here suddenly decreases and then rises slightly,and its peak is inversely proportional to the reentry altitude,which means that the surface pressure is maximum at 40 km,at about 1000 Pa.

The stagnation point pressure and the distance of the stagnation point from the center of the ARD sphere at four altitudes are given in Table 4,whereRrepresents the distance between the stagnation point and the center of the sphere. It can be found that with a constant angle of attack, the stagnation point pressure of the ARD first increases and then decreases with decreasing reentry altitude, but the reentry altitude has less influence on the position of the stagnation point on the spherical surface.

Table 5 shows the stagnation point pressure and the distance of the stagnation point from the center of the ARD sphere at 50 km. For reentry capsules, the angle of attack is generally negative. For the convenience of illustrating the regularity of variation,Table 5 shows the absolute value of the angle of attack, expressed as the magnitude of the angle between the capsule centerline and the flight trajectory. It can be seen that the stagnation point pressure increases when the angle of attack continues to increase at a certain altitude, but the increment is not very large,which indicates that the effect of changing the angle of attack only is not very obvious for peak pressure of the ARD.However,the distance between the stagnation point and the center of the ARD sphere becomes more and more distant with increase in the angle of attack,and constantly deflects to the shoulder of the windward side.

Table 4. Stagnation point position of the ARD at four altitudes.

Table 5. Stagnation point position of the ARD at 50 km.

Figure 6 shows distributions of the Mach number around the ARD at different altitudes. Comparing the value of the Mach number for the four altitudes,the maximum Mach number gradually decreases with decreasing altitude. The low-Mach region of the flow-field is mainly concentrated in the shock layer around the head of the ARD and the recirculation zone at the back cover. After the shoulder flow expansion plays a major role and the gas flow velocity gradually increases and the Mach number of the flow-field at the windward side of the rear cone is relatively high.

Figure 7 shows the distribution of the surface heat flux of the ARD with an angle of attack of-20°at different flight altitudes(e.g.,4 km,5 km,6 km,and 7 km)in units of W·m-2.With decreasing altitude the surface heat flux first increases and then decreases. Note that a negative is added to the heat flux scale to indicate direction, and the absolute value of the number is usually used to compare magnitudes. The region of maximum heat flux is mainly distributed on the windward side of the shoulder. This is a result of airflow heat transfer when the accelerated airflow due to the expansion wave passes across the shoulder,predicting that the region will suffer a higher heat load. After the shoulder,the surface heat flux decreases rapidly; however, there is a small increase in the windward side of the back cover so thermal protection measures are needed for both parts during vehicle design.

Fig.6. Mach number(Ma)distribution around the ARD flow-field at different altitudes: (a)40 km,(b)50 km,(c)60 km,(d)70 km.

Fig.7. Distribution of surface heat fluxes of the ARD at different altitudes: (a)40 km,(b)50 km,(c)60 km,(d)70 km.

Figure 8 shows the curves for variation of the translational-rotation temperature and vibration-electron temperature around the ARD standing point line at different altitudes. The horizontal axis refers to the vertical distance from the shock layer to the head of the ARD, which is parallel to the flight path of the ARD, and the vertical axis refers to the temperature in units of K.

Fig. 8. Temperature variation curves around the ARD standing point line at different altitudes: (a)translational-rotation temperature,(b)vibrationalelectronic temperature.

Comparing the temperature variation curves for the four altitudes, it can be seen that the peak values of both temperatures between the ARD and the shock wave gradually decrease with decrease in the ARD reentry altitude, and the variation trends of the translational-rotation temperature and vibration-electron temperature corresponding to the four altitudes are consistent. The two energy modes of the gas are not fully excited because there is not enough energy outside the shock wave layer, so both temperatures are low to start with. After a certain point,when the fluid passes through the shock wave layer,both temperatures increase rapidly,and the translational-rotation temperature always increases before the vibration-electron temperature. This indicates that the fluid is in thermal non-equilibrium here, because the vibrationelectron mode requires more energy to be fully activated and the energy transfer rate between the two modes of the particles is limited. Since the surface of the ARD is an isothermal surface and its temperature is much lower than the shock layer temperature, a larger temperature gradient exists around the surface of the ARD,causing a rapid decrease in both temperatures near the surface.

6. Conclusion

In this paper, the ARD with an angle of attack of-20°is used as the research target, and a combination of the twotemperature model and turbulence model is used to simulate the flow-field of the ARD in the blackout region at altitudes of 40 km,50 km,60 km,and 70 km. The changes of parameters including pressure, Mach number and temperature are discussed in detail. In summary,the following conclusions are drawn:

(i)During the reentry of the ARD with an angle of attack of-20°the high-pressure region of the flow-field is concentrated in the shock wave layer around the windward side of the head, while the flow-field in the low-Mach region is mainly concentrated in the shock wave layer around the head and the recirculation zone at the back cover. After the shoulder of the ARD, flow expansion plays a dominant role, the Mach number of the windward side of the rear cone increases and the flow-field pressure decreases rapidly.

(ii) During reentry of the ARD with an angle of attack of-20°, the surface high-pressure region is concentrated in the windward side of the head and is distributed in a toroidal shape;the value of the pressure is inversely proportional to the radius of the torus. Meanwhile,the stagnation point is less affected by the reentry altitude. After passing the shoulder, the surface pressure decreases rapidly,but there is a small increase on the windward side of the back cover,and its value is much smaller than the surface pressure on the head.

(iii) The translational-rotational temperature increases before the vibrational-electron temperature in the shock wave layer around the head, and a thermal non-equilibrium effect occurs. Due to the existence of this temperature gradient,both temperatures will again decrease rapidly near the surface of the ARD.At the same time,the shoulder,as well as the windward side of the back cover, will suffer a large thermal load and require appropriate thermal protection measures.

Acknowledgments

Project supported by the National Natural Science Foundation of China (Grant No. 12175177) and the China Postdoctoral Science Foundation (Grant No. 2021M693889). All the calculations involved in this study were carried out on the Tianhe-2 Supercomputer at the National Supercomputer Center in Guangzhou,China.