DSMC modeling of rarefied ionization reactions and applications to hypervelocity spacecraft reentry flows

The DSMC modeling is developed to simulate three-dimensional (3D) rarefied ionization flows and numerically forecast the communication blackout around spacecraft during hypervelocity reentry. A new weighting factor scheme for rare species is introduced, whose key point is to modify the corresponding chemical reaction coefficients involving electrons, meanwhile reproduce the rare species in resultants and preserve/delete common species in reactants according to the weighting factors. The resulting DSMC method is highly efficient in simulating weakly inhomogeneous flows including the Couette shear flow and controlling statistical fluctuation with high resolution. The accurate reliability of the present DSMC modeling is also validated by the comparison with a series of experimental measurements of the Shenzhou reentry capsule tested in a low-density wind tunnel from the HAI of CARDC. The obtained electron number density distribution for the RAM-C II vehicle agrees well with the flight experiment data, while the electron density contours for the Stardust hypervelocity reentry match the reference data completely. In addition, the present 3D DSMC algorithm can capture distribution of the electron, N+ and O+ number densities better than the axis-symmetric DSMC model. The introduction of rare species weighting factor scheme can significantly improve the smoothness of the number density contours of rare species, especially for that of electron in weak ionization case, while it has negligible effect on the macroscopic flow parameters. The ionization characteristics of the Chinese lunar capsule reentry process are numerically analyzed and forecasted in the rarefied transitional flow regime at the flying altitudes between 80 and 97 km, and the simulations predict communication blackout altitudes which are in good agreement with the actual reentry flight data. For the spacecraft reentry with hypervelocity larger than the second cosmic speed, it is forecasted and verified by the present DSMC modeling that ionization reactions will cover the windward capsule surface, leading to reentry communication blackout, and the communication interruption must be considered in the communication design during reentry in rarefied flow regimes.


Introduction
With rapid development of manned spaceflight, lunar and deep space exploration, recyclable spacecraft reentry speed has been increasing from the first to the second cosmic velocity, even higher. Therefore, typical hypervelocity flow phenomena such as thermo-chemical nonequilibrium, ionization and communication blackout occur also in rarefied flows [1,2]. For instance, the Chinese lunar exploration capsule launched in 2014 returned with a speed of 11.2 km/s, while the Stardust recorded the fastest man-made spacecraft reentry speed of 12.8 km/s [3][4][5]. Under these typical hypersonic conditions, the extremely strong bow shock wave around the vehicle leads to intensive physico-chemical nonequilibrium process, thus the atmosphere needs to be considered as an ionized gas mixture with vibrational excitation, dissociation and electronic excitation. To improve communication design of spacecraft during hypervelocity reentry, we need to accurately predict the electron density distribution at the complex three-dimensional (3D) surface of spacecraft, which has been a long-standing research challenge [6,7] of the extreme flows in the field of fluid mechanics. At low altitude, the ionization effect is usually treated as the conventional chemical reaction using a continuum flow model such as the Navier-Stokes equation [8,9]. However, the continuum flow assumption for the Navier-Stokes equation becomes invalid for rarefied gas at high altitude. Ground test facilities can hardly reproduce the crucial non-equilibrium flows of extremely highspeed reentry in non-continuum atmospheric environment [10]. Although kinetic methods including the direct solver of kinetic models [11][12][13][14][15], the discrete velocity models [16][17][18][19][20][21][22], the gas-kinetic unified algorithm [23][24][25][26][27], the unified flow solver [28,29], and the BGK-type discrete unified gas-kinetic schemes [30,31] can describe non-equilibrium gas flows, they are still under development and computationally expensive for hypersonic flows. Currently, the most widely used method for rarefied hypersonic flows is the Direct Simulation Monte Carlo (DSMC) [32][33][34] method. In 1987, Bird [35] utilized the DSMC method to simulate the rarefied air ionization for the first time, where the ionization process was treated using 11 molecular species and 41 chemical reactions, and the electrons in close proximity to the ions were kept in order to enforce charge neutrality. This pioneering work established a framework for simulation of ionization effect in rarefied flows. Carlson & Hassan [36] proposed a scheme to compute the electric field intensity, however, this scheme would suffer significant statistical fluctuation due to the scattering of charged particles. Both approaches have only been applied for one-dimensional stagnation problems because of the computational costs, thus it was difficult and unsuccessful for the DSMC results to compare with flight experiment data [35,37]. The surface geometry of aircraft will strongly influence flow structures including shock standoff distance etc. The flow along the stagnation line can be significantly different for different spacecraft configurations even if all other parameters were the same. Therefore, the real spacecraft geometries should be taken into account for accurate calculations, and it is essential to be able to simulate rarefied air ionization around 3D spacecraft with complex geometries.
Recently, multidimensional DSMC simulations of rarefied ionization have been under rapid development. Ozawa [38] used DSMC and conventional computational fluid dynamics (CFD) methods to simulate rarefied ionization process on Stardust using a simplified axisymmetric configuration. Also using axisymmetric configuration, Boyd [39] applied DSMC to study the rarefied ionization of RAMC-II and introduced a special numerical scheme to trace species and electron-ion enforced association. Fan [40] proposed a trace species separation (TSS) scheme to compute the electron density distribution around the RAMC-II which was also simplified to an axisymmetric configuration. Unlike the probabilistic models of traditional DSMC, TTS calculates the particle number densities of trace species from those of common species based on reaction rates. Shevyrin [41] compared the differences among the TCE [35], KSS [42] and QK [43] models for rarefied ionization, where the electron is not modeled as particle and its density is computed as the sum of those ions.
It is still a tough challenge to effectively capture the rarefied ionization effect for 3D hypervelocity reentry flow, as it needs to overcome two difficulties in order to enable realistic complex 3D simulations. First, restricted by computational resources, the number of simulated molecules in each cell is very limited, often less than the order of ten, thus each cell may not even have one single electron. The presence of an extremely small number of rare species can cause significant numerical fluctuations, which will affect both the determination of the number of molecular collision pair and the electron density sampling. Second, for the simulation of the 3D plasma effect, the electron-ion enforced association would lead to additional computing load and restrict the electron movement, while Carlson and Hassan's method [36] suffers great statistical fluctuation. Meanwhile, if electron is simulated as another species with realistic electron mass, restriction on the time step would be unacceptable. However, Boyd [44] proposed to amplify the mass of electron by three orders of magnitude, i.e. the electrons have almost the same order of speeds with other species after collisions, thus the restriction on the time step can be removed.
In this work, we extend our DSMC method to simulate rarefied atmosphere ionization around 3D complex spacecraft [5,33,45], in which the above artificial amplification in electron mass is adopted, and the mass of ions is adjusted to guarantee mass conservation. A weighting factor scheme for rare species is proposed for the rarefied ionization process, which will mitigate fluctuations on both the determination of the number of molecular collision pair and the electron density sampling. The chemical reactions involving ionization for 11 species air are treated to ensure reproduction of the rare species in resultants and preservation/deletion of common species in reactants according to the weighting factors. The second part of this paper is about the application of the newly developed DSMC modeling method in the hypervelocity reentry flows. The resulting DSMC code is applied to study reentry ionization of spacecraft considering their 3D complex shapes including the RAM-C II and the Stardust, and the effectiveness of the weighting factor scheme is investigated. The rarefied atmosphere ionization during the Chinese lunar exploration capsule reentry is simulated and analyzed with the actual flight test. The remaining parts of this paper are organized as follows: The DSMC algorithm for ionization reaction is introduced in Section 2. The weighting factor scheme for trace species in modelling ionization reaction is proposed in Section 3. The algorithm validation and analysis are presented in Section 4. The numerical simulation and analysis for the rarefied ionization around the Chinese lunar capsule is in Section 5, followed by conclusions and the expectation of the future work in Section 6.

DSMC modeling for rarefied ionization reaction
The DSMC method simulates the physics of the Boltzmann equation by following the motions and collisions of a large number of simulated molecules, and their movement and collisions are decoupled over a time step that is smaller than the local mean free time [5,33,45,46]. The spatial position coordinates, velocity components and internal energies of the simulated molecules are stored in the computer during computation; and each simulated molecule is assigned a sequence number of local cell index. Their values alter with time due to the movement of simulated molecules and the interactions with boundaries as well as other molecules. When a collision occurs, postcollision velocities are calculated using conservation of momentum and energy. The chemical reaction is implemented as a process of binary molecular collision, usually written as where A and B are denoted to pre-collision molecules, while C and D are the postcollision ones. The chemical reactions with ionization involve dissociation, exchange, associative ionization and its inverse, direct ionization and charge exchange. E represents atom for dissociation or electron for direct ionization. In all other cases, E is null. The reaction rate coefficients are usually expressed in the Arrhenius form where Λ and η are constants, E a is the reaction activation energy, κ is the Boltzmann constant and T is temperature. For E c > E a , the above Λ and η relate to reaction probability P by Here, σ T is the total collision cross section, ε is the symmetric factor whose value is 2 when A equals to B, or 1 otherwise. T ref is the reference temperature and Γ is the Gama function. ξ is the average internal degrees of freedom, ω AB is the average viscosity index, and m r is the reduced mass.
Due to the significant mass difference between electron and usual species, the magnitude of electron speed would be 4 to 5 orders higher than that of other species. If a unified time step is chosen for the motion of electron and molecules, the time step would be too small, resulting in unaffordable computational cost. And if different time steps are used for electron and other species, respectively, although the computational costs could be reduced, complex additional modifications are required for DSMC simulations. An efficient empirical solution is the artificial mass amplification of electron as proposed by Boyd [44]. The principle is to reduce electron speed to nearly the same as other species while the speed of other species will not be affected significantly. If the amplification is too small, the electron speed would still be too much higher than other species. On the contrary, large changes to ion speeds will occur if the amplification is too large as the summation of electron and ion mass should be equal to the original atom or molecule. To striking an appropriate balance, an amplification of electron mass by 3 orders of magnitude is chosen here. Therefore, the electron virtual mass is about 3-4% of atoms, and reduction of ion mass is also up to 3-4%. Moreover, to keep the reaction probability in Eq. (3) and to match the changes of electron mass, the chemical reaction rates involving electrons should be reduced by ffiffiffiffiffiffiffiffiffiffi 1000 p , as shown in Table 1. For other chemical reactions, the Gupta model [47,48] is adopted.
3 Weighting factor scheme for trace species in modeling ionization reaction Since each cell may not have one single electron, significant fluctuations on both the number of molecular collision pair and the electron density sampling need to be properly dealt with. Instead of using various complex tracking approaches for trace species [34,49,50], a new weighting factor scheme is developed here to model ionization reaction. Different weighting factors are used for different species: the weighting factor is 1 for N 2 , O 2 , NO, N and O, while the W s , which is no more than 1, for N þ 2 , O þ 2 , NO + , N + , O + and electron e, respectively. For convenience, a unified W s is adopted for all trace species in this paper. The weighting factor mainly affects the computation of collisions in two ways, namely, the computation of collision numbers as well as mechanism, and the realization of chemical reactions. Therefore, the following modifications will be introduced in the DSMC code for the simulation of ionization reaction.

Weighting factor model for particle collision
Based on the non-time-counter (NTC) technique of the DSMC method [51], the ionized air species are divided into two groups: group P for electrons only, and group Q for all other species [52]. Neglecting the inter-particle collisions in group P, we only need to deal with the inter-particle collisions in group Q, and the collisions between group P and Q. Without weighting factor, the number of inter-particle collisions in group Q is where N Q is particle number in group Q. The number of collisions between the particles in groups P and Q is Where N P is particle number in group P, F N denotes the real particle number that a simulating particle represents, V c is the cell volume, and (σ T c r ) QQ, max and (σ T c r ) PQ, max are the maximum value of the product of cross section and relative speed updated at each time step.
The introduction of weighting factors will increase N P and N Q as well as N Q in Eq. (4) and Eq. (5). Traditionally, an array CS(2, i, j) is used to record the particle numbers of species j in cell i. Due to the weighting factor, the N Q in Eq. (4) should now be modified as where W j is the weighting factor for species j, and N P in Eq. (5) becomes where W e is the weighting factor for electron. As binary collisions are considered in DSMC simulations, the post collision velocities of two particles with identical weight are where U ! m is the velocity of center of mass of the collision pair, and g ! is their relative velocity. If the two collisional particles have different weighting factors, the velocity of the particle with smaller weight is determined by Eq. (8a) or Eq. (8b). However, the velocity of the particle with the larger weight factor calculated by Eq. (8a) or Eq. (8b) needs to be modified with probability which means the velocity of the particle still keeps unchanged by the probability of 1 − P.
The above approach cannot exactly conserve momentum and energy in every collision step. However, violence of conservation laws is negligible even after a sufficiently large number of collisions.

Weighting factor scheme for chemical reactions
For the 11-species air model, the chemical reactions process involving rare species can be classified into four types: the reaction between the common species such as associative ionization N + O → NO + + e which produces the rare species; the reaction between the common and rare species such as direct ionization N + e → N + + e + e which produces the rare species; the reaction between the rare species such as association ionization inverse NO + + e → N + O which produces the common species; and the reaction between the rare and the common species such as N 2 þ N þ →N þ N þ 2 which exchanges charge and produces new species. These four types of reactions are treated separately.
For N + O → NO + + e, it reproduces the products of NO + and e by 1/W s copies. And for the inverse association reaction NO + + e → N + O, the products are reserved by a probability of W s . For N + e → N + + e + e, whose reactant is N and the product e is also copied by a factor of 1/W s times. However, for N 2 þ N þ →N þ N þ 2 , in which usual and trace species are involved in reactants and products, N 2 in the reactants is eliminated and N in the products is reserved by a probability of W s .
The purpose of trace species weighting factor is to increase the number of simulated particles and reduce the statistical fluctuation for the species participated in ionization reaction. Obviously, according to the above weighting factor scheme, the numbers of simulated particles of those species are amplified by a factor of 1/W s .

Computational test of DSMC modeling on controlling statistical fluctuation
In order to validate the present DSMC modeling rules for simulating weakly inhomogeneous flow in simple gases, the planar Couette shear flow from two parallel infinite plates with relative motion placed at a distance of H is simulated and compared with the results of the gas-kinetic unified algorithm (GKUA) [23][24][25][26][27] where we directly solve the Boltzmann-type velocity distribution function equation. The Couette flow is a simple and benchmark boundary-value flow problem. They have the same surface temperature T W and move at a speed U W in opposite directions, as shown in Fig. 1. In this simulation, we choose Ar with temperature T ref = 273K, T W = T ref and pressure P = 0.01atm. The variable soft sphere (VSS) model [32] with ω = 0.81 and α = 1.4 is applied to describe the interaction between molecules. For this Couette flow, we will investigate the change of the statistical noise of the flow velocity with the plate velocity U W , the total number of simulated particles N t and the number of sampling N. In order to verify the accuracy of the DSMC results, Fig. 2 gives a comparison of the flow velocity distributions between the present DSMC simulations and the Unified Algorithm (GKUA) [23,53] for direct solution of Boltzmann model equation, where the dimensionless flow velocity distributions U/U W along the distance y/H in the upper half of this Couette flow are presented. The DSMC results in Fig. 2 for the Knudsen number of Kn = 0.1128 are respectively corresponding to (a) U W = 0.2Ma (with N = 2 × 10 4 and N t = 1000), (b) N t = 64000 (with U W = 0.025Ma and N = 2 × 10 4 ), and (c) N = 128 × 10 4 (with U W = 0.025Ma and N t = 1000). It can be shown that since the inflow speed is generally known, in order to reduce the influence of the statistical noise on flow velocity, the  The planar Couette shear flow in rarefied transition regime with different Knudsen numbers has been simulated, and Fig. 3 gives the flow velocity distribution of the present DSMC and the GKUA as well as numerical analysis of linearized Boltzmann equation [54] related to Kn = 0.1128, 1.128 and 11.28 respectively. And the hard sphere (HS) molecular model is used in the DSMC simulation. The plate speed is 0.2Ma, the total number of simulated particles is N t = 1000 and the number of sampling is N = 12 × 10 4 . It shows that these results from three approaches agree perfectly. This figure also implies that as the Knudsen number increases, the effect of rarefied gas gradually gets stronger, and the slip velocity between the gas and the wall increases. The comparisons from Figs. 2 and 3 confirm the accurate reliability and practicality of the present DSMC modeling in solving the rarefied flows of simple gases which can control

Shenzhou reentry capsule aerodynamic characteristics of hypersonic reacting gases
To verify the validity of the present DSMC algorithm with reaction model and rare species weighting scheme for hypersonic non-equilibrium reentry flow, a parallel code has been developed using mesh self-adaption. A series of experiments for the Shenzhou reentry capsule are performed in the low-density wind tunnel as shown in Table 2, to evaluate the present algorithm. Here, M ∞ , P 0 , T 0 and Kn ∞ represent the inflow Mach number, the total pressure, the total temperature and the inflow Knudsen number, respectively. The experiment model has a length of 3.75 cm and a bottom radius of 1.888 cm. The flow media is N 2 , and the angles of attack (AOA, α(°)) are set to be 10°, 15°, 18°, 20°, 22°, 30°and 40°for each case. During the simulation, the fixed wall temperature 300 K and diffusive reflection boundary condition are adopted. We focus on the aerodynamic characteristics such as the axial force coefficient C A , the lift coefficient C L and the lift-to-drag coefficient K, as shown in Fig. 4. The axial force coefficient C A reduces as the AOA increases, with errors of 2.2% for an AOA of 10°and 7.9% for an AOA of 40°. Meanwhile, due to the increase of inflow Knudsen number, the rarefaction effect becomes more significant and leads to the increase of C A . The lift coefficient C L has negative values, which does not vary monotonously with AOA but reaches its minimum around an AOA of 30°. The behavior of the lift-to-drag coefficient K is similar as that of C L . The excellent agreement between the DSMC calculations and the experimental data demonstrates the reliability of the present algorithm.

Ionization reactions for the RAM-C II at 81 km
The RAM-C (Radio-Wave Attenuation Measurement for Study of Communication Blackout) II flight experiment [44,55] performed in the late 1960s provides unique experimental data to study ionization in rarefied atmosphere, which is usually used as a classical benchmark case to verify computational models [40][41][42]. The entry speed of the RAM-C II vehicle is 7.8 km/s, and its configuration and the probe locations are illustrated in Fig. 5. Here n ∞ , T ∞ and V ∞ represent the inflow molecule number density, the inflow temperature and speed, respectively. The maximum electron density was measured at the four probe locations. The altitude of 81 km is chosen as a typical test case. The computational setup is summarized in Table 3. The surface temperature is set to be 1000 K with diffuse reflection of wall boundary condition. Cell self-adaption and sub-cell are used during the simulation, thus the spatial discretization requirement can be satisfied naturally. The flow field around the RAM-C II vehicle is simulated and tested by the present DSMC algorithm with the weighting factor approach of trace species, in which the weighting factor W s is set to be 0.01. All simulations are performed on a workstation server, with 24 Intel cores and a memory of 192Gb on each node. Figure 6 shows the contour distribution of the flowfield macroscopic properties including number density, temperature, pressure and Mach number, around the RAM-C II vehicle for the flying state at the height of 81 km, where the "Without weighting" and "Weighting scheme" in all of the corresponding captions including the text below denote the modeling of ionization reaction "without" and "with" considering the weighting factor scheme from Section 3 for trace species, respectively. Figure 6 also shows that the weighting factor scheme has negligible effect on the distributions of macroscopic flow parameters. The air is seriously compressed around this vehicle flying at  the hypersonic reentry speed. The detached shock is formed near the surface, and the temperature behind the shock reaches as high as 15,000 K. Therefore, chemical reactions, especially ionization, occur in the surface area after the shock wave. Figure 7 presents the contour of electron number density distribution. It is observed that the electron number density concentrates on the surface area behind the shock wave, and the maximum electron number density is about 10 18 /m 3 , only about 0.01% of the molecular number density (comparing with Fig. 6a), which means the ionization is extremely weak. Moreover, Fig. 7 demonstrates the effect of the usage of rare species weighting factor scheme on the electron number density distribution, where great fluctuation exists when weighting factor scheme is not used, while the electron number density distribution is smooth enough when the weighting factor scheme is adopted. Interestingly, the electron number density distribution is almost the same as that of NO + , see Fig. 8. This is the reason why the 7-species model [56] is appropriate for moderate hypervelocity reentry, where the main ionization source is associative ionization [44]. The electron number density distributions at the four probe locations as shown in Fig. 9 are obtained by the present DSMC simulations (Symbol "◇"), which are also compared with the flight data (Symbol "○") [55] measured by onboard microwave reflectometers. The computed and measured results at all the measurement locations are in excellent agreement at the flying altitude of 81 km, which indicates that the present DSMC algorithm can capture the extremely weak electron number density distribution.
To demonstrate the computing efficiency, we carried out computations for three cases: without chemical reaction, chemical reaction without ionization, and chemical reaction with ionization. For each case, the program runs 20,000 steps to achieve the assumed stable state before sampling, then another 40,000 steps to obtain the routine statistical sampling for the macroscopic flow variables. The CPU time is given in Table 4, which shows that the computational cost of each case differs slightly. Comparing with that in Ref. [40], where a DELL cluster with 64 processors was used and 20 h of CPU time was consumed, the present proposed model and the developed code can be said to be highly  efficient. On the other hand, the case in Ref. [40] is 2D computation, and the present computation is 3D with different CPU time, in which the total number of simulated particles in the present case is about 5 × 10 7 , while that in Ref. [40] is 10 7 .

Ionization reactions during the Stardust hypervelocity reentry
Since the Stardust has the highest manmade reentry speed, ionization is more significant. It has a nose radius of 0.22 m and a head with half cone-angle of 59.5 o . The total length is 0.5 m and the bottom radius is 0.25 m, and the configuration is illustrated in Fig. 10. Similarly, the computational conditions are shown in Table 5. To compare with Ref. The weighting factor W s of trace species is set to be 0.1. Surface temperature is set to be 1000 K, and the wall is set to be diffusive. Cell self-adaption and sub-cell are also used. The macroscopic flow variables are shown in Fig. 11, which suggests that the air is more seriously compressed as the free-stream speed is up to 12.8 km/s. However, the detached shock is formed further away from the Stardust surface than that of the RAM-C II, as the Stardust is a larger blunt vehicle with the head half cone-angle of 59.5 o and a wider area of strong disturbance. The highest temperature behind the shock reaches above 36,000 K, indicating that serious ionization reactions occur near the surface region. The negligible influence of weighting factor scheme on macroscopic flow variables is also confirmed, see Fig. 11. The major concern is still the electron number density distribution in the flowfield, especially near the front surface of the Stardust, as shown in Fig. 12a. The maximum electron density is larger than 5.0 × 10 20 /m 3 , which implies that the ionization can reach 3% for spacecraft reentering with an extremely high speed of 12.8 km/s. Figure 12b and c also show the number density distributions of N + and O + . In contrast with the case of the RAM-C II, the main ionization source is the direct ionization of N and O atoms. The ionization reaction area is larger, which almost envelops the head cone surface so that the communication blackout is formed for the Stardust reentry flow at H = 80km flying with a speed of 12.8km/s. Due to the relatively stronger ionization, the fluctuation of rare species without weighting scheme is smaller. However, the usage of rare species weighting factor scheme can still improve the smoothness of their distribution contours, especially, for the extremely rare species such as O 2 + . Figure 13 shows the comparison of electron number density distribution in the flowfield between the present 3D calculation (lower part) and the axis-symmetric simplified model (upper part) [3]. In general, the electron number density contours agree well, although small difference exists in its thickness and maximum value. At this zero angle of attack, both methods predict almost the same standoff distance for the detached bow wave. However, quantitative comparisons of the electron, N + and O + number density profiles along the stagnation line show significant difference, see   Fig. 12, different from the prediction of the axis-symmetric model [3], where a significant amount of NO + still exists. From the computation and comparison of the ionization reactions of the RAM-C II and Stardust reentry flows, the present three-dimensional DSMC algorithm can capture fine flow field better than the axissymmetric DSMC model. The weighting factor scheme of trace species is essential to describe weak flow information of the electron, N + and O + number densities. Fig. 13 Comparison of the present 3D calculation and axis-symmetric results [3] for Stardust at 80 km

Rarefied ionization analysis and forecast for the Chinese lunar capsule reentry
The Chinese lunar capsule returned on November 1, 2014, and its configuration is illustrated in Fig. 15. Since it re-entered the earth atmosphere with double jump reentry from the second cosmic velocity to the first cosmic speed, its atmospheric ionization was significant and suffered a much longer time period of the radio black-out with two access communication interruptions from high rarefied flow to near space flight environment. Due to the existence of its ear piece accessories and the flying angle of attack, flow simulations of a 3D realistic geometry are required. Here, the ionization process is numerically forecasted by the present DSMC modeling with three different stages: (1) the first reentry, between the altitudes of 90 and 97 km, (2) the first jump-out, between 80 and 90 km, (3) the second reentry, between 80 and 85 km. For all the three cases, diffusive wall reflection with fixed temperature 300 K is used, and cell self-adaption and sub-cell are also adopted.

Ionization analysis for the first reentry
According to its trajectory, typical related parameters are listed in Table 6, with an initial cell system of 240 × 250 × 250, and 60 CPUs carry out the parallel computation   Table 6, the corresponding electron number density distributions are shown in Fig. 16, with maximum values 5.216 × 10 19 /m 3 , 5.047 × 10 18 /m 3 and 1.988 × 10 18 /m 3 . Since the corresponding maximum particle number densities are 9.35 × 10 21 /m 3 , 3.19 × 10 21 /m 3 and 2.05 × 10 21 /m 3 , the ionization degrees are extremely low. For these cases, the ionization degree can be defined as the ratio between electron numbers and particle numbers in the cell, which shows the orders of 1% to 1‰ for these heights. Similar to the case of the Stardust, the main ionization source is the direct ionization due to the collisions between N/O atoms and neural atoms/molecules, and this point can be supported by the number density distributions of N + and O + in Fig. 17. Here, the sum of N + and O + number densities is almost equal to that of the electron, and the contribution of NO + is in the order of 1%, while the contributions of N 2 + and O 2 + are lower for at least one order of magnitude.
For the interest of radio communication design, the main focus is the electromagnetic shielding. And for the weak plasma formed by the rarefied atmosphere ionization, the plasma electron oscillation frequency [57] is determined by where n e is the electron number density, e is the electron charge, m e is the electron mass, and ε 0 is the dielectric constant. If n e is expressed in unit of cm − 3 , the corresponding radio frequency can be approximated as Since the capsule returns with an AoA between -20 o and -25 o , and there is a certain angle between the capsule and the ground radar station, the electron number density contour that envelopes the whole bottom part of the capsule will lead to radio blackout. The ground radar frequency is about 2.3GHz, which is plotted in Fig. 18 and reveals that the complete communication blackout will start at about between 95 and 97 km, and the flight observation value is about 97.1 km. Fig. 18 The corresponding electromagnetic wave frequency of lunar exploration type capsule at different altitudes for the first reentry (Unit: Hz)

Ionization analysis for the first jump-out
Similarly, the typical parameters are shown in Table 7, with an initial cell system of 480 × 500 × 500. Here, the weighting factor for rare species is 0.01, while other conditions are the same as in Section 5.1. The electron number densities and electromagnetic wave frequency corresponding to Figs. 16 and 18 are given in Figs. 19 and 20. Since that the reentry speed reduces to lower than the first cosmic speed, the electron number density at 90 km decreases for about one order of magnitude, and the degree of ionization reduces significantly. Thus, the demand for weighting factor scheme behaves more serious, so that the smaller weighting fact is employed for the rare species in the present computation. According to Fig. 20, the complete communication blackout will occur at the altitude of 85 km, or a little higher.
The main ionization source is still the direct ionization due to the collisions between N/O atoms and neural atoms / molecules, here the further discussion will not continue due to the length of the paper.

Ionization analysis for the second reentry
Similarly again, typical parameters are shown in Table 8. All other parameters are set as same as those in Section 5.2. The electron number densities and electromagnetic wave frequency are given in Figs. 21 and 22. The electron number density distribution in Fig. 21 differs slightly with that in Fig. 19, and the reason is that the speed difference is not significant at the same altitude, comparing with the first jump-out. Figure 22 shows that, for the S-band electromagnetic wave at the frequency of 2.3GHz, the complete communication blackout will take place at about 85 km, or a little lower. And the conclusion for main ionization source is still the same.

Comparison of the present numerical forecast with the flight observation
In summary, the comparisons of complete communication blackout between the present forecasts and the flight observations can be shown in the following Table 9. Obviously, the differences between the present numerical forecasts and the observed value from the flight test of the Chinese lunar exploration capsule launched in 2014 are on the order of 2 km, which reveals the excellent capabilities of the present numerical forecast of the ionization DSMC modeling based on the rare species weighting factor schemes.

Conclusion
Rarefied gas ionization effects on aerodynamics of spacecraft hypervelocity reentry are of the greatest practical interest. A DSMC modeling with the weighting factor approach of trace species has been developed, which enables 3D simulations of rarefied air ionization with complex configurations of spacecraft. The present algorithm has been used to compute the weak rarefied atmospheric ionization of the RAM-C II vehicle, the moderate ionization of the Chinese lunar capsule, and the strong ionization of the Stardust at the flying altitude ranging from 80 km to 97 km.
By comparing the flow field electron number density with those obtained by the flight measurement experiments at different probe locations, it is shown that the present DSMC results are in excellent agreement with the flight measured data. In comparison with the previous axisymmetric simulations for hypervelocity reentry flow at a zero angle of attack, the present DSMC algorithm can predict the distributions of the electron better. In addition, the present method is also able to simulate the flows over complex geometry of spacecraft at varying angles of attack. For that reason, it is also applied to analyze rarefied ionization for the reentry flow of the Chinese lunar capsule at the flying altitudes between 80 and 97 km. In the rarefied transitional flow regime, the ionization reaction is found to cover the windward capsule surface and leads to communication blackout, which is confirmed by the actual reentry flight test data.
The introduction of rare species weighting factor scheme can significantly improve the smoothness of rare species number density contours, especially for that of electron in weak ionization case, while it has negligible effect on the macroscopic flow variables. From the present DSMC computations of the reentry flows of the RAM-C II, the  Stardust, and the Chinese lunar capsule, we find that ionization is extremely weak for a spacecraft reentry at the first cosmic velocity, while it is significant to cause communication blackout for a spacecraft reentry with the second cosmic speed, an important designing issue to prevent communication interruption. For example, the maximum electron number density is only about 0.01% of the molecular number density for the RAM-C II with V ∞ = 7800m/s, while it can be up to 3% for the Stardust reentry at an extremely high speed of V ∞ = 12.8km/s. As this work is only the beginning of numerical forecast for reentry communication interruption of the thermochemical ionization reaction, further investigations on the hypervelocity reentry flows around recyclable spacecraft with real gas effects involving internal energy excitation, thermochemical ionization radiation non-equilibrium reactions, need to be studied in the future.