On the calculation of the electron temperature flowfield in the DSMC studies of ionized re-entry flows

Currently available procedures of electron temperature calculations in studying ionized flows around reentry spacecraft by the direct simulation Monte Carlo (DSMC) method are analyzed. It is shown that the heat conduction of electrons is not taken into account in these procedures. The contributions of various effects to the electron energy balance are calculated by an example of the RAM-C II capsule, and a numerical solution of the electron energy conservation equation is obtained, which refines the electron temperature distribution used in the DSMC computations. A method of coupled calculation of the electron temperature within the framework of the continuum approach and modelling of ionized gas flow by the DSMC method is proposed.


Introduction
The direct simulation Monte Carlo (DSMC) method [1] is traditionally used for studying high-altitude aerothermodynamics of space vehicles. At altitudes above 80 km, where the mean free path of molecules becomes comparable with the spacecraft size, the level of thermochemical nonequilibrium of the flow becomes significant. Conventional methods of computational fluid dynamics based on the Navier-Stokes equations become inapplicable, and the DMSC method turns out to be the main tool for the numerical analysis of gas flows around space vehicles.
A shock layer is formed around a space vehicle entering the atmosphere with a high velocity. The gas temperature in the shock layer is determined by the total enthalpy of the free-stream flow and can reach several tens of thousands degrees. One of the phenomena inherent in such high-enthalpy flows is ionization. The ionized shock layer reflects radio waves and leads to blackout, which prevents radio communication with the space vehicle. At high velocities of atmospheric reentry (more than 10 km/s), when the ionization degree becomes rather high, ionization affects both the gas flow around the space vehicle and the convective heat flux onto the vehicle surface. The heat flux value has to be taken into account in thermal protection system design. Therefore, the study of the plasma surrounding the reentry vehicle is of significant interest, and effective models that can be used in DSMC computations have to be developed for investigating high-altitude aerothermodynamics of prospective space vehicles [2].
The main difficulty in the DSMC simulation of unbound electrons in the plasma is caused by the small mass of electrons; for this reason, it is necessary to consider processes of different scales for electrons and heavy particles. Therefore, the thermal motion of electrons is not modeled in the DSMC applications to the reentry flows, and the electron concentration is determined from the ion concentration [3,4]. Various approximations and procedures that take into account the energy exchange between electrons and heavy particles in different ways are used for calculating the electron temperature distribution. Criteria for validation of these procedures are rather uncertain.
In the present work, modeling of the electrons in the DSMC modeling of the hightemperature ionized shock layer around the reentry vehicle is considered. The goal of the paper is to analyze available simplified models of electrons as applied to the calculation of plasma near the reentering orbital spacecraft and to formulate new algorithms for obtaining physically justified data on plasma parameters.
The available electron models in the DSMC method are analyzed by using appropriate conservation equations. Comparisons of the electron energy balance in these models and conservation equations written in the most general form are performed.
For the conservation equation written in the general form, the influence of various effects on the electron temperature distribution is evaluated by an example of a typical ionized flow around an orbital reentry vehicle. This evaluation shows which effects have to be taken into account for obtaining a correct electron temperature distribution and which effects can be neglected. The flow around the RAM-C II capsule at an altitude of 73 km is simulated as a typical example [5]. This flow was chosen because of a small effect of charged particles on neutral particles. For this reason, the distributions of macroparameters of the neutral species are not affected by changes in plasma parameters, which is convenient for the analysis of DSMC results.
The measurements of the number density of electrons along the RAM-C II capsule surface [6] were widely used for verifying the possibility of an adequate description of high-enthalpy hypersonic air flows. As an example, we can mention [1] and also [7], which contains a detailed numerical and theoretical analysis of this flow. Advanced capabilities of taking into account nonequilibrium effects in solving Navier-Stokes equations for the flow around the capsule at descent trajectory altitudes of 61 and 71 km were demonstrated in [8].
Based on this analysis, a combined procedure of DSMC computation of ionized flows is proposed, which consists of solving the electron energy conservation equation within the framework of the continuum approach and calculating heavy particles by a standard DSMC approach. The main requirement in the development of this new procedure is the consistency between the predicted electron temperature distribution and the electron energy conservation equation written in the general form.

Conservation equations for electrons
The equations of conservation of mass, momentum, and energy in a steady highenthalpy flow can be written in the following form [9]: div n e u e ð Þ¼R ð1Þ div ρ e u ei u e Here n e , u e = (u e1 , u e2 , u e3 ), and p e are the number density (m -3 ), velocity, and pressure of electrons; m e and e are the electron mass and elementary charge. The mass density is ρ e = m e n e , and E is the electric field strength.
In Eq. (1), R is the number of electrons generated and lost in ionization and recombination reactions per unit volume per unit time. In Eq. (2), τ ij is the stress tensor and F c is the volume force corresponding to electron momentum relaxation in elastic collisions with heavy particles. In Eq. (3), W = 3n e k B T e /2 is the volume density of thermal energy of electrons, S e is the source term that describes energy production and relaxation processes, and q e is the heat flux vector. The equation of state for electrons is p e = n e k B T e .
Five relations (1)-(3) are written in the most general form and involve all processes that can affect the electron temperature distribution in the shock layer. These relations include 18 unknowns (n e , u e , T e , τ ij , q e , F c , S e ) and the electric field E; therefore, they are not enough for finding the solution. At the same time, if distributions of all flow parameters are obtained with some approximate model, then substitution of these parameters into the conservation equations for the electron component allows one to estimate how accurately these parameters' distributions satisfy the general conservation relations.
Estimation of the terms in the energy conservation Eq. (3) can be used to determine the most influential effects on the electron temperature distribution. The equation can be applied to analyzing the energy balance regardless of the way the flow fields were obtained. Examples of such an analysis for orbital and superorbital reentry are presented in [4].
Various procedures of the DSMC simulations of ionized re-entry flows take into account only those processes of electron energy exchange that are considered as the most important ones by the authors of procedures. It is usually assumed that a simplified conservation equation that involves only the processes taken into account in the procedure ensures a sufficiently accurate description of energy exchange. In this case, the electron temperature distribution obtained in the approximation used satisfies the corresponding approximate conservation equation and must satisfy the general Eq. (3) with a sufficient degree of accuracy. Thus, in order to analyze the accuracy of a DSMC electron model, it is necessary to estimate the correspondence between the numerical solution obtained with this model and Eq. (3).
In the present paper, we consider various procedures of the DSMC modeling of electrons. Accuracy of the approximate relations for the electron energy corresponding to these procedures is analyzed.
For the procedures implemented by the authors of the present paper, the distribution T e obtained in DSMC computations is refined by solving the equation of electron energy conservation (3) numerically for a stationary case. The DSMC macroparameter flowfields are used to calculate the terms of the equation. As a result, accuracy of the model is analyzed and the processes responsible for the differences with the general energy conservation law are determined.

Simulation of electrons in the DSMC method
The most relevant processes for the electron energy balance in hypersonic air flows are energy gain and loss during the formation and recombination of the ion-electron pair S i , relaxation with gas translational temperature due to elastic collisions S eT (collisions with neutral particles and ions), electron-vibrational exchange with diatomic molecules S eV (S e = S i + S eT + S eV ), and also energy transfer due to convection and heat conduction [9].
Electrons are much lighter than gas molecules; therefore, the velocity of thermal motion of electrons is much greater than the velocity of heavy particles. At the same time, the motion of electrons is determined by the electric field, which constrains the diffusion of electrons and does not allow the charges to be separated to distances greater than the Debye length (the so-called ambipolar diffusion).
DSMC modeling of the electron component together with ions and neutral particles, combined with finding the electric field, allow ionized flows to be fully studied within the DSMC approach (PIC-DSMC, [3]). However, it makes the numerical algorithm prohibitively computationally expensive (though the DSMC method itself requires extensive computational resources). Hence, only ions are modeled by the DSMC method, while the plasma neutrality condition and the ambipolar diffusion approximation are used for electrons. This approximation implies that the electron pressure gradient is compensated by the electric field: If the electron density at a distance L changes by Δn e , then the charge separation is is the Debye length. The charge separation can be neglected at λ D ≪ L [10]. In most cases, the Debye length is of the order of fractions of a millimeter. The electron density gradient in the shock layer around the reentry vehicle is directed across the shock layer. If the shock layer thickness L is about 0.01-1 m, and Δn e is about n e , then the neutrality condition is satisfied fairly well. While the charge neutrality approximation is used in all DSMC models of electrons, the approximation used for energy exchange of electrons are different. Some of the models are discussed below.

Model of Bird, 1987 [11]
The plasma neutrality assumption is used to determine the electron density in the simplest model of ionized air for the DSMC method [11]. In this model, an electron is assumed to be connected to the ion with which it was generated, and dissociative recombination reactions, which depend on the electron energy, are ignored. Therefore, the electron velocity coincides with the ion velocity, and Eq. (1) coincides with the continuity equation for ions, the electron momentum is assumed to be equal to zero, and the ion energy is not considered in this model.

Model of Farbar and Boyd 2010 [3]
The model proposed in [12] was further developed in [3]. In this model, each electron has a velocity used for modeling electron interaction with heavy particles, but moves with the mean velocity of ions in the cell, which is determined by sampling over a certain time interval. If the instantaneous ion velocity is denoted by u + and its averaged value by u þ avg , then the energy conservation equation for this model takes the form where S clmb and S cld are the terms that describe electron collisions with ions and neutral particles. Comparing Eqs. (3) and (5), we can conclude that the model [3] ignores the heat conduction of electrons.

Models of Shevyrin et al., 2016 [4, 5]
In our previous studies, we used two different models [4,5]. The first model (for weakly ionized air plasma) implies intense resonant e-V relaxation on nitrogen molecules N 2 : T e = T v (N 2 ). Moreover, it is assumed that Eq. (3) contains only the relaxation term. The second model was formulated under the assumption that T e = const, which implies intense heat conduction. These models are the limiting cases of the conservation equation div q e ð Þ ¼ S eV ð6Þ and correspond to the limits as τ eV → 0 and τ eV → ∞, respectively. Neither model can reproduce the intermediate situation.

Model of Kusov, 2016 [13]
The recently developed DSMC code [13] simulates electron collisions with ions and neutral particles. At each step, mixing of electrons between ions to which they are attached is simulated. In addition, collisions of each electron with all electrons and all ions in a cell are performed (obviously, providing e-T relaxation and Maxwellization of the electron distribution function). Thus, the model [13] satisfies the conservation equation The magnitude of the heat flux induced by mixing of electrons between ions in a cell depends on the cell size and time step, which may fail to be a fairly good approximation of heat conduction of electrons. The author did not provide the electron temperature distribution, but it follows from his paper that the electron temperature should be close to the translational temperature of ions.
None of the simplified schemes considered above describes such an important process as heat conduction of electrons correctly; therefore, the predicted electron temperature distribution may be appreciably different from the real distribution satisfying Eq. (3).

Analysis of energy exchange of free electrons
Let us first simplify Eqs. (1)-(3). The plasma is assumed to be neutral: n e = n + . The plasma satisfying the ambipolar diffusion approximation (4) is considered. It is also assumed that viscous friction in the electron component is absent [14], therefore, τ ij = 0; and the electron bulk velocity rapidly relaxes to the gas velocity owing to collisions: u e = u. The energy, associated with electron bulk velocity, is much smaller than the thermal energy, and the term ρ e u 2 /2 in Eq. (3) can be neglected.
Using Eq. (4), the combination associated with the work of pressure and electric field div(u e p e ) + en e (E,u e ) can be written as p e div(u e ). Thus, Eq.
Heat conduction of electrons depends on the velocity distribution function of electrons. If the number of electrons in the domain with sizes equal to the Debye length is large, intense interaction of electrons occurs. Then the distribution function is weakly nonequilibrium, and the Fourier approximation for heat conduction is valid: where κ is the thermal conductivity of electrons. In the present work, we use the value taken from [4] κ(T e ) = 3.67·10 − 11 T 5=2 e W/(m K) based on [10]. The relaxation terms S eT and S eV are calculated on the basis of the relaxation expression S = 3/2 n e k B (T * -T e )/τ, where T * is the equilibrium temperature and τ is the relaxation time. The values of τ eT and τ eV are taken from [15,16], respectively.

Parameters of DSMC simulations
By example of the flow around the RAM-C II capsule at an altitude of 73 km, the relation between the simplified model [5] T e = T v (N 2 ) and Eq. (8) was considered. The capsule is a cone 1.295 m long with an apex angle of 9°and spherical bluntness with a radius R = 15.24 cm. The flow around the capsule with a velocity of 7650 m/s is calculated in the axisymmetric formulation with the free-stream parameters corresponding to the altitude of 73 km: the number density of the gas is 9.7 × 10 20 m − 3 , and the temperature is 210 K. The free-stream molar fractions of O 2 and N 2 are assumed to be 21.7% and 78.3%, respectively. We consider a nine-species air mixture (N, O, N 2 , O 2 , NO, N 2 + , O 2 + , NO + , e) with the set of reactions used in [5] (15 dissociation reactions, 4 exchange reactions, and 6 reactions of associative ionization/dissociative recombination). The chemical reactions were modeled with the use of the TCE model and chemical reaction rate constants from [17]. In contrast to [5], where the number of reactions of dissociative recombination (AB + + e = A + B, where A and B are atoms) was calculated macroscopically, with the use of a modified Arrhenius dependence for the reaction rate constant, here we use the method [18], where the probability of recombination of a randomly chosen electron and ion is taken to be P(E) = CE B , where E = E e + E int is the sum of the electron energy and internal energy of the ion. The diffuse reflection with complete accommodation of energy and neutralization of the ion charge was assumed on the capsule surface; the surface temperature was 1000 K.
To suppress the statistical error, we use statistical weights for charged components, which increases the number of simulation molecules for ions and electrons by a factor of 64.

Solution procedure
For the considered flow, we obtain the distribution of macroparameters on a rectangular Cartesian grid in the domain shown in Fig. 1 To speed up the solution convergence, the macroparameters are first interpolated to a coarse grid (in the present case, 224 × 96). When the convergence is reached, the resultant distribution of T e is transferred to the initial grid, and the procedure of solving the problem was continued. The calculations were carried out using a single-core program on the Intel i7-7700 3.6GHz processor. To obtain the result presented in this paper, we first performed 50,000 iterations on a 224 × 96 grid, which took 60 s, and then 90,000 iterations on a 2240 × 960 grid, which took 3 h. The calculation of the macroparameter fields by the DSMC method was carried out using the MPI technology in the multiprocessor code, and 98 CPU days were spent in total.

Results of calculations
The resultant distribution of T e is shown in Fig. 4. A comparison of the distributions of T e and T v (N 2 ) is of major interest because the latter is used as T e in the DSMC simulation (model T e = T v (N 2 )). The main change in T v (N 2 ) (Fig. 3) is directed across the shock layer away from the capsule surface. However, as it is seen from Fig. 4, the change in the electron temperature across the shock layer is rather small, and the temperature mainly changes along the capsule surface or along the streamlines. The value of the peak temperature located near the stagnation point is lower by 680 K than the nitrogen vibrational temperature: the maximum values of T v (N 2 ) and T e are 9400 and 8720 K, respectively. The distribution of T e and source terms corresponding to the processes of heat conduction div(q e ) and e-V relaxation S eV along the stagnation line l is shown in Fig. 5. It should be noted that the l axis is directed away from the body, i.e., opposite to the x axis: l = 0 corresponds to the wall and l > 2 cm actually corresponds to the free stream. Figure 5а also shows the profile of T v (N 2 ). A comparison of different terms in the equation shows that the main energy exchange is provided by heat conduction and e-V relaxation. The terms that describe other processes are smaller by two or more orders of  At l > 1.6 cm, the value of T e does not decrease below 8100 K, forming a plateau in the curve T e (l). The main reason for the formation of this plateau is the small density of N 2 in the upstream region and the small rate of e-V exchange. The electron density in this region decreases by orders of magnitude as compared to the peak values; correspondingly, the volume density of electrons becomes smaller. Therefore, the heat flux generated in the process of heat conduction of electrons leads to intense heating of electrons in this region.
The distribution of the source terms along the stagnation line in Fig. 5b testifies to intense heating of electrons both at the distance l < 2.4 from the wall and at l > 1 cm.
The heat flux q e toward the wall along the stagnation line reaches 100 kW/m 2 (at l = 2.7 mm). It should be noted that this is a rather small value under the conditions considered in the study (the convective heat flux at the stagnation point is 1.4 MW/m 2 ), and the flow pattern is not much affected by this heat flux.
The temperature gradient near the side surface of the capsule is significantly smaller than that near the stagnation point; therefore, the conditions of energy exchange of  electrons may be significantly different. Let us consider the distribution of parameters along the normal to the capsule surface at x = 0.232 m (the line is shown by the arrow in Fig. 4). As the plot in Fig. 6а shows, the temperature of electrons along this line varies within narrow limits (from 4250 to 4990 K). Therefore, the values of the source terms in Eq. (8) are smaller almost by two orders of magnitude than those on the stagnation line (Figs. 5b and 6b, respectively). However, the difference in the peak values of T e and T v (N 2 ) reaches 1000 K, i.e., it is greater than the difference on the stagnation line.
Let us consider the validity of assumptions made in the beginning of the present section. The correspondence of Eqs. (8) and (3) is provided if both conditions are valid: 1. Approximation of ambipolar diffusion (4), which is valid if the Debye length λ D is small; 2. Approximation of the heat conduction relation (9), which is valid if the local distribution of electron velocities is close to equilibrium. If the number of electrons n e λ D 3 in the volume with dimensions equal to the Debye length is large, then the electrons intensely interact with each other, and the local distribution may be expected to be close to the equilibrium form.
The field of the local values of λ D is shown in Fig. 7 and confirms that the ambipolar diffusion approximation is indeed valid almost in the entire computational domain. The field of n e λ D 3 (Fig. 8) shows that the condition n e λ D 3 > > 1 is also satisfied every- where. The insets of Figs. 7 and 8 show the distributions of the corresponding variables along the stagnation line. Thus, the validity of approximations (8) and (9), as well as the approximation of local neutrality of the plasma in DSMC procedures for this flow is estimated. This estimation shows that Eq. (8) approximates the general Eq. (3) quite well. As a whole, the results obtained in the present study show that the distribution of T e , which is the solution of the general energy conservation Eq. (3) is significantly different Fig. 6 Electron temperature and vibrational temperature (а); profile of S eV and heat flux divergence (b) along the line indicated in Fig. 4 from the approximation T e = T v (N 2 ) used in the present DSMC simulations. Therefore, for computing the electron temperature distribution around the RAM-C II capsule, both e-V relaxation and heat conduction of electrons have to be taken into account simultaneously.

Conclusions
A method of a detailed analysis of energy balance for the electron component in weakly ionized high-enthalpy air flows around space vehicles entering the Earth atmosphere has been proposed. The method is based on the electron energy conservation equation, which is written in the most general form. Accurate distributions of the electron density and temperature must satisfy this equation.
The presented analysis of plasma models for the DSMC method consists of two parts. First, the electron energy conservation equation satisfied for each particular DSMC plasma models was derived. A comparison of this equation with the general equation demonstrates the accuracy and applicability of the considered numerical model. Second, after calculating the flowfields of macroparameters in the shock layer, one can calculate the terms in the general conservation equation and the new (refined) distribution of the electron temperature. A comparison of the new distribution and that used in the DSMC simulations allows one to draw conclusions on the accuracy of the DSMC solutions and on the necessity of taking into account these or those processes of electron energy exchange in simulation procedures for a particular flow.
The analysis of the currently used simplified models of electrons in the DSMC method (proposed both by the authors of the present work and by other researchers) shows that  none of the models is capable of relevant reproduction of heat transfer by free electrons. The only approach completely satisfying the conservation equation in the general form (fully kinetic PIC-DSMC formulation) is too computationally expensive and is not suitable for solving complicated 3D problems of spacecraft aerothermodynamics. It is demonstrated by an example of the flow around the RAM C-II capsule that the electron temperature distribution significantly differs from the approximation previously used by the authors: T e = T v (N2) [5,19]. It is shown that the difference is determined by heat conduction by electrons, which has to be taken into account in addition to e-V relaxation. Evaluation of other energy exchange processes (Coulomb and elastic collisions, convection, and changes of energy in reactions) shows that these processes can be neglected under the conditions considered in the present study.
The main novelty of the present work is the analysis of electron energy exchange in the full-scale formulation (where the dimension of the energy conservation equation coincides with the spatial dimension of the problem considered). This analysis provided a possibility of the maximum possible accuracy of the description of the processes near the side surface of the capsule, which cannot be described if the flow is considered in the one-dimensional formulation.
The numerical data obtained in the present simulations and the results of the analysis of the plasma models for the DSMC method show that the currently available simplified models of the plasma have to be improved. A possible aspect for improvement is consideration of electrons within the framework of the continuum approach, i.e., solving the electron energy equation by finite difference methods together with DSMC procedures.
The following algorithm for obtaining the electron temperature distribution in the DSMC method is proposed: 1. Macroparameter flowfields (including electron density) are computed by the DSMC method in one of simplified approximations; 2. The electron energy conservation equation is solved, and a new distribution of T e is obtained; 3. The DSMC computation is continued with the new distribution of T e until the steady state is reached. 4. Steps 2 and 3 are repeated until the convergence is reached.
The calculation of T e performed in the present work for the flow around the RAM-C II capsule is the result of the first iteration in this algorithm. The main advantage of the proposed method is the possibility of sufficiently accurate simulation of the electron component by the continuum approach which requires much less computational resources than kinetic particle methods.
While only typical orbital entry of blunt body RAM C-II is considered in the present analysis, the proposed method seems to be promising for other application cases. For example, for a spacecraft reentering the Earth's atmosphere with speeds of more than 10 km/s, the energy exchange of electrons can be substantially different from orbital reentry. The most interesting question then could be a comparative analysis of the orbital and superorbital entries, showing a qualitative difference between these cases. Such a study is planned in the future works.