PIC-DSMC Simulation of a Hall Thruster Plume with Charge Exchange Effects Using pdFOAM

Abstract

To develop technologies for the stable operation of electric propulsion systems, the effects of charge exchange (CEX) on the exhaust plume of a Hall thruster were studied using the particle-in-cell direct simulation Monte Carlo (PIC-DSMC) method. For the numerical analysis, an OpenFOAM-based code, pdFOAM, with a simple electron fluid model was employed. In an example problem using the D55 Hall thruster exhaust plume, the results showed good agreement with experimental measurements of the plasma potential. In the results, CEX effects enhanced Xe⁺ particle scattering near the thruster exit. However, due to the increase in the plasma potential with CEX effects, fewer Xe²⁺ particles were near the thruster exit with CEX effects than without CEX effects.

1. Introduction

Electric propulsion (EP) systems were first applied in real space systems such as satellites and spacecrafts in the early 1960s. Initially, EP was used in simple functional demonstrations of some space missions. Subsequently, due to the high specific impulse of EP systems, the demand for EP gradually increased. In 2018, approximately 50% of geostationary satellites launched around the world were equipped with EP systems [1]. In particular, due to the need for high specific impulse, Hall thrusters and ion thrusters are used in most of the EP systems employed in spacecraft.

Most Hall thrusters and ion thrusters employ xenon as a propellant. Depending on the thruster type, Xe atoms are ionized to Xe⁺ or Xe²⁺ through various ionization processes. These ions are accelerated by electromagnetic forces and generate thrust. Since not all the Xe atoms are ionized within the thruster, the exhaust plume of the thruster typically consists of neutral atoms, positive ions and free electrons. Neutral atoms are not accelerated by the electromagnetic force, so they slowly exit the thruster. However, xenon ions are accelerated to high speeds by the electromagnetic force. Slow xenon atoms and fast xenon ions mix and collide in the exhaust plume. Depending on the collision conditions, atom–ion collisions may result in momentum exchange (MEX) or charge exchange (CEX) [2,3,4,5]. In the processes of CEX, fast ions accept electrons and become fast atoms, and slow atoms lose electrons and become slow ions. Thus, various species with different particle speeds are present in an exhaust plume. In particular, slow ions significantly influence the electric field distribution and spreading of the exhaust plume [6,7]. Depending on the conditions, the speed of ions in an exhaust plume can reach several km/s. In addition, widely spreading fast ions can impinge on the body of a satellite and cause sputtering and erosion in important devices of the satellite. Since the divergence angle of the plume is directly related to damage of various devices in a spacecraft, the effects of CEX on the exhaust plume should be accurately estimated to ensure sufficient lifespans of satellites and other spacecraft [7,8,9].

The particle-in-cell (PIC) method has been widely used for the analysis of exhaust plumes of Hall thrusters [10]. The PIC method is used to consider electromagnetic forces, calculate the acceleration of simulated particles in node cells, and move the particles based on acceleration by electromagnetic forces. However, the PIC method does not consider collisions between simulated particles. Therefore, to analyze regions with high number densities of particles, such as the region immediately behind the thruster exit, collision modeling should be added to the PIC method.

There are two methods to include particle collision models in the PIC method. The first is the particle-in-cell Monte Carlo collisions (PIC-MCC) method [11], and the second is the particle-in-cell direct simulation Monte Carlo (PIC-DSMC) method [12]. The PIC-MCC method considers collisions between a simulated particle and a target particle cloud near the simulated particle. In this method, the target particle cloud is not a simulated particle, so the PIC-MCC method does not trace the target particle after a collision. Therefore, using target particle clouds may result in errors in conservation of energy and momentum [13]. However, the PIC-DSMC method selects two simulated particles as candidates for a collision. After the collision, the PIC-DSMC method traces both simulated particles, and thus, energy and momentum are conserved.

Due to its important influence on satellite lifespans, substantial research has been conducted on the CEX phenomenon for decades [2,3,4,5,6,7,8,9]. However, it is hard to find numerical research on the difference between Xe⁺ and Xe²⁺ particle behaviors in the charge-exchange zone of the exhaust plume using the PIC-DSMC method.

In the present study, we used the PIC-DSMC method to investigate the characteristics of exhaust plumes of a Hall thruster with CEX effects. For the numerical analysis, pdFOAM, an OpenFOAM-based open source code, was employed [14]. Although pdFOAM was developed in 2017, it has not yet been used to analyze the plumes of electric propulsion systems. Therefore, the present study is the first attempt to use pdFOAM for Hall thruster plume analysis. In the following sections, the numerical method and the solution procedures used in this study are explained, and a sample problem for exhaust plume analysis using the D55 Hall thruster [3,4] is introduced. By comparison with experimental results, the present numerical method is validated. Finally, CEX effects on the exhaust plume of the D55 Hall thruster are addressed.

2. Numerical Method

In the present study, pdFOAM, a PIC-DSMC code developed in the OpenFOAM platform, was used for numerical analysis of the behavior of the plume of a Hall thruster [14]. Similar to other OpenFOAM-based codes, pdFOAM can be applied in numerical simulations associated with both structured and unstructured grids and used in parallel and serial computing.

Figure 1 shows the solution algorithm of pdFOAM. Starting from the upper-left box of the figure, the simulated particles are inserted into the domain using the inlet conditions or moved within the domain using the velocity and acceleration conditions of the existing particles. After updating the cell occupancies with the moved particles, collision partners are selected from each cell. If the collision partners satisfy the collision condition, then the particles collide. Depending on the conditions, the colliding particles may react. After updating the cell occupancies again, the electric charges of the ions are transformed to values at the cell nodes. Then, using the charges at the cell nodes, the electric field can be solved. Based on the solved electric field, the total electromagnetic forces on each particle are calculated. Using the boundary conditions and the particle conditions, the upper-left box can be revisited as in the first step. Further capabilities of pdFOAM are summarized in

Table 1. Model capabilities of pdFOAM [14].

Type	Models
Collision partner selection method	Transient Conglomerated Cell (TCC)
Collision model	Hard Sphere (HS), Variable Hard Sphere (VHS), Variable Soft Sphere (VSS)
Reaction model	Larsen–Borgnakke, Quantum-Kinetic (Q-K)
Charge transform method	Nearest Volume (NV), Composite Linear Volume (CLV)

, and more details can be found in the literature [14].

Before pdFOAM was applied to the numerical analysis of the exhaust plume of Hall thrusters, some modifications were made to the code. Since pdFOAM does not consider CEX and MEX, new modules to include those effects were added. Based on the measurements of Miller et al. [15], the following collision cross sections for CEX were used.

$${\mathsf{\sigma }}_{CEX}\left(X{e}^{+},Xe\right)=87.3-13.6\log \left(E_{lab}\right)$$

(1)

$${\mathsf{\sigma }}_{CEX}\left(X{e}^{2+},Xe\right)=45.7-8.9\log \left(E_{lab}\right)$$

(2)

Boyd and Dressler analyzed the differential cross section and found that the probability of CEX is 0.5 at scattering angles exceeding 2° [16]. Thus, they concluded that the same scattering data for CEX can be employed for MEX interactions. Therefore, the collision cross section for MEX was set to be the same as that for CEX [4,8].

$${\mathsf{\sigma }}_{MEX}\left(X{e}^{+},Xe\right)={\mathsf{\sigma }}_{CEX}\left(X{e}^{+},Xe\right)$$

(3)

$${\mathsf{\sigma }}_{MEX}\left(X{e}^{2+},Xe\right)={\mathsf{\sigma }}_{CEX}\left(X{e}^{2+},Xe\right)$$

(4)

Collisions between atoms were considered by using the variable hard sphere (VHS) model [17].

In axisymmetric 2D analysis, the simulated particles near the axis represent fewer particles than do those located far from the axis. Therefore, using radial weighting factors for the simulated particles could improve the efficiency of the computation [18]. Since the original pdFOAM does not employ radial weighting factors, the code was modified to utilize weighting factors for efficient axisymmetric simulations [19].

To consider the effects of free electrons, the pdFOAM code provides both hybrid kinetic (HK) and fully kinetic (FK) models [14]. The HK model considers free electrons as a fluid, while the FK model considers free electrons as particles [20]. The HK model in pdFOAM is a simple electron fluid model with the assumption of an isothermal, unmagnetized, inertialess electron fluid to utilize the Boltzmann relation (Equation (5)) [14]. The electron number density was calculated with the HK model in pdFOAM from Equation (6), and the potential distribution of the field was calculated by solving Poisson’s equation (Equation (7)) with the assumption of a quasi-neutral free stream plasma.

$$\varphi -{\varphi }_{ref}=\frac{k_B{T}_e}{q_e}\ln \left(\frac{n_e}{n_{e,ref}}\right)$$

(5)

$$n_e=n_{e,ref} exp\left[\frac{q_e\left(\varphi -{\varphi }_{ref}\right)}{k_B{T}_e}\right]$$

(6)

$${ϵ}_0{\nabla }^2\varphi -q_e{n}_{e,ref}exp\left[\frac{q_e\varphi }{k_B{T}_e}\right]=-q_i{n}_i$$

(7)

However, if many free electrons and ions from active ionization events are present, as in Hall thruster plumes, then the electron number density from Equation (6) can differ from the ion number density obtained by particle tracing, and it may break the quasi-neutrality of the exhaust plume. Therefore, in the present study, to ensure the quasi-neutrality of the plume, the electron number density was set equal to the charge number density of ions, as in Equation (8) [2]. The plasma potential was obtained from the Boltzmann relation as in Equation (5) [4,8].

$$n_e=Z_i{n}_i$$

(8)

Basically, this study employed the same rules for the Boltzmann relation and quasi-neutrality that were used by the original pdFOAM code, but their application sequence was changed.

3. Results and Discussion

3.1. Description of a Sample Problem

As a sample problem, the characteristics of the exhaust plume of the D55 Hall thruster were analyzed by using pdFOAM. This study also used experimental data for the D55 Hall thruster that have served as benchmarks in various numerical studies [3,4,21,22,23,24,25]. Figure 2 shows the configuration of the D55 Hall thruster. The thruster consists of an annular chamber with an inner diameter of 50 mm and an outer diameter of 60 mm and employs 4.76 mg/s xenon gas as a propellant; its specific impulse is approximately 1819 s [21].

Figure 3 shows the computational domain and boundary conditions used to observe the exhaust plume behaviors of the D55 Hall thruster. As shown in the figure, a 2D axisymmetric computational domain was employed. As in the study by Boyd [21], the flow through the cathode was neglected, and thus, the 2D axisymmetric condition was satisfied. In this study, the numerical test conditions were taken from the experiments of Domonkos et al. [22] and Zakharenkov et al. [23].

Table 2. Thruster exit conditions for the numerical simulation in this study.

Case #	Species	Number Density (m−3)	Temperature (K)	Velocity (m/s)
1	Xe	3.8 × 1018	750	281
	Xe+	3.6 × 1017	46,400	15,000
	Xe2+	9.0 × 1016	46,400	21,300
2	Xe	4.6 × 1018	750	281
	Xe+	2.4 × 1017	46,400	15,000
	Xe2+	6.0 × 1016	46,400	21,300

shows the exit conditions of the thruster used for the numerical simulation. In the table, Case 1 and Case 2 correspond to the conditions measured by Domonkos et al. [22] and Zakharenkov et al. [23], respectively. For the boundary conditions in Figure 3, the diffuse wall condition was applied to the wall, and the outflow condition was applied to the outlets (all particles reaching the outlets were eliminated from the simulation). The initial condition was set to a vacuum state, and the background pressure was assumed to be zero.

Several numerical studies have reproduced these experimental results, although their thruster exit conditions differed. In a previous experiment [22], the electron temperature was measured at approximately 10 eV near the thruster exit but approximately 3 eV downstream. For the electron temperature in the Boltzmann relation (Equation (5)), Boyd [21] used an upstream value of 3 eV, but Choi [3] and Wang et al. [4] used a downstream value of 10 eV. In this study, the median value of those in the studies by Choi and Wang, i.e., 6.5 eV, was used. Furthermore, Choi [3] and Wang et al. [4] employed nonuniform potential profiles that ranged from 0 V to 147 V at the thruster exit. However, in this study, the median between the maximum and minimum values, i.e., 73.5 V, was used. Other flow properties at the thruster exit were the same as those employed by Choi et al. [24]; they are summarized in Table 2.

In the present numerical simulation, the minimum cell size was 2 mm, and the time step was 5 × 10⁻⁸ s. The steady state was reached after 60,000 time steps, and the number of particles at the steady state was approximately 930,000.

3.2. Comparison and Validation

A numerical simulation of the sample problem defined in the previous section was performed using pdFOAM. Figure 4 shows the distribution of the Xe⁺ number density and the plasma potential around and downstream of the thruster. As shown in the figure, xenon ions exited the thruster and spread into space. Since the thruster was annular, the spreading ions were concentrated at the front of the axis of symmetry. In the Boltzmann relation (Equation (5)), the plasma potential is proportional to the logarithm of the electron number density. Because of quasi-neutrality, the electron number density was equivalent to the ion number density. Therefore, in the figure, the plasma potential was high and similar to the Xe⁺ number density at the front of the axis of symmetry.

Figure 5 and Figure 6 compare the distributions of the plasma potential vs. radial distance at x = 10, 50, and 500 mm and with experimental data. Due to the availability of the experimental data [22,23], the conditions of Case 1 and Case 2 were used for the results in Figure 5 and Figure 6, respectively. The numerical results of Choi [3] were added to the figures for comparison. Overall, the results of this study were in better agreement with the experimental data than those obtained when Choi used the simple electron fluid model. However, except for the case where x = 10 mm, the results obtained when Choi used the detailed electron fluid model showed better agreement with the experimental results than did the results of this study. The simple electron fluid model assumes an isothermal, unmagnetized, inertialess electron fluid. However, since the electron temperature in the exhaust plume actually greatly varies, the simple electron model may exhibit inaccurate results [26]. Depending on the location in the plume, the magnetization effects can also be severe [27]. To consider the variations in electron properties, Choi solved the continuity, momentum and energy equations for the electron fluid in the detailed electron fluid model [3]. Therefore, to obtain more accurate results, the detailed electron fluid model should be used. In addition, if the simple electron fluid model is to be used for some reason, validation through comparison with experimental results should precede its use.

Figure 7 and Figure 8 compare the distributions of the ion current density and electron number density at x = 10 mm with experimental data. As shown in Figure 7, the ion current density in this study was in good agreement with the experimental data. In Figure 8, the electron number density in this study showed some discrepancies from the experimental data. In particular, the result of this study showed a much lower density near the axis. The isothermal electron fluid model generally overestimates the local temperature and plasma potential at the axis; hence, it underestimates the electron density. However, the trend of the electron number density showing a relatively high value near x = 0.02 m was similar to the experimental results.

3.3. Effects of Charge Exchange

As explained in the Introduction, CEX may significantly influence spreading of the exhaust plume and the configuration of the electric field near the thruster exit. To observe these effects, another numerical simulation for Case 1 was conducted without CEX effects. In this case, momentum exchange between atoms and ions was modeled by using the conventional VHS model [17].

Figure 9 shows contours of the Xe⁺ number density and the plasma potential for the case without CEX effects. Even without CEX effects, contours of the number density of Xe⁺ in the downstream area were similar to those of the case with CEX effects shown in Figure 4. However, near the thruster exit (dotted circle), Xe⁺ particles scattered less without CEX effects (Figure 9) than with CEX effects (Figure 4). Since the plasma potential is a function of the logarithm of the electron number density, the pattern of the plasma potential was similar to that of the Xe⁺ number density. As shown in the figure, even without CEX effects, the contours of the plasma potential in the downstream area were almost the same as those of the case with CEX effects shown in Figure 4. However, near the thruster exit (dotted circle), there was a narrower zone of high potential in the case without CEX effects (Figure 9) than in the case with CEX effects (Figure 4).

To observe CEX effects more quantitatively for the cases with and without CEX effects, the plasma potentials along the radial distance at x = 10, 50, and 500 mm are shown in Figure 10. Close to the thruster exit (x = 10 and 50 mm), there was a clear difference between the cases with and without CEX effects in the plasma potential near the axis of symmetry. However, this tendency decreased with increasing distance from the thruster exit, and there was little difference between the cases with and without CEX effects at x = 500 mm.

Figure 11 shows contours of the Xe²⁺ number density for the cases with and without CEX effects. However, the distributions of the Xe²⁺ number density showed different patterns from those in the previous figures. In Figure 11, regardless of CEX effects, the number densities of Xe²⁺ in the downstream area showed almost the same pattern. However, near the thruster exit (dotted circle), the Xe²⁺ particles scattered less with CEX effects (top figure) than without CEX effects (bottom figure), which was the opposite of the pattern in the previous figures.

To clearly identify the differences in particle distributions, contours of the Xe⁺ number density for the cases with and without CEX effects are shown in Figure 12. Near the thruster exit (dotted circle), Xe⁺ particles were more scattered in the case with CEX effects (top figure) than in the case without CEX effects (bottom figure).

To determine the reason for these differences, the numerical results were analyzed more quantitatively. Figure 13 shows the accumulated number of particle pairs that underwent charge exchange in Case 1. In the results at 5 msec, the accumulated number of particle pairs was approximately one million for CEX between Xe⁺ and Xe, while it was only approximately 5500 for CEX between Xe²⁺ and Xe. This means that CEX between Xe⁺ and Xe was very active and produced many low-speed Xe⁺ ions and high-speed Xe atoms, while CEX between Xe²⁺ and Xe rarely occurred, which resulted in a small number of low-speed Xe²⁺ ions and high-speed Xe atoms.

To verify the results, the axial variations in the magnitude of the velocity of the Xe atoms and ions along the thruster channel center are shown in Figure 14. With CEX effects, the average velocity of the particles was 888 m/s for the Xe atoms, 16.9 km/s for the Xe⁺ ions, and 23.8 km/s for the Xe²⁺ ions. In other words, due to the strong electric forces on the doubly charged Xe²⁺ ions, the Xe²⁺ ions are much faster than the other particles. Therefore, sputtering and erosion by impingement of Xe²⁺ ions on the satellite body can be much more severe than that by other particles. As shown in Figure 14a, the Xe atoms were faster with CEX effects than without CEX effects. However, in Figure 14b, the Xe⁺ ions were slower with CEX effects than without CEX effects. This was because CEX between Xe⁺ and Xe generated a number of low-speed Xe⁺ ions and high-speed Xe atoms. However, as shown in Figure 14c, the speed of Xe²⁺ upstream (axial distance < 0.1 m) was not remarkably different for the cases with and without CEX effects. This meant that CEX between Xe²⁺ and Xe rarely occurred, and low-speed Xe²⁺ ions were rarely generated. Therefore, even with CEX effects, scattered Xe²⁺ ions were not noticeable in Figure 11.

However, even considering negligible CEX between Xe²⁺ and Xe, the reason for the opposite pattern in Figure 11 remained unclear. Basically, CEX effects should enhance particle scattering. However, in Figure 11, Xe²⁺ particles scattered less with CEX than without CEX, which meant that CEX effects reduced the scattering of Xe²⁺ particles. This opposite pattern can be explained by the distributions of the plasma potential as follows.

Figure 15 shows contours of the plasma potential for the cases with and without CEX. Near the thruster exit (dotted circle), there was a wider zone of high potential with CEX than without CEX. This was due to the highly scattered slow Xe⁺ ions after CEX near the thruster exit. As explained in Figure 13, CEX between Xe²⁺ and Xe was not remarkable in Case 1. Therefore, CEX did not enhance the scattering of Xe²⁺. However, as shown in Figure 15, the plasma potential near the thruster exit was higher with CEX than without CEX. When activating CEX, Xe⁺ ions are more scattered radially because of their frequent CEX collisions. Xe²⁺ ions, on the other hand, suffer very rare collisions, so that their radial spreading is more affected by the local electric field than by their own CEX collisions. Since, with CEX, the radial electric field at the plume periphery is lower, Xe²⁺ ions tend to spread less in the radial direction, thus justifying the observations.

4. Conclusions

In this study, the effects of CEX on the exhaust plume of a Hall thruster were investigated using the OpenFOAM-based open source code pdFOAM. For the numerical analysis, the code was modified to include CEX and MEX and improved to utilize weighting factors for efficient simulation of a thruster with axial symmetry. The simple electron fluid model was employed to model electrons.

As a sample problem, the exhaust plume of the D55 Hall thruster was solved and compared to benchmark data. The results showed good agreement with experimental data and with numerical results obtained from other studies in the literature.

To observe CEX effects on the exhaust plume of the D55 Hall thruster, a numerical simulation without CEX was performed. In the results, the effects of CEX were remarkable near the thruster exit, where Xe⁺ particles were scattered more than in the simulation without CEX. In addition, the zone of high potential was wider due to the scattered ions with CEX than without CEX.

However, the patterns in the number density contours for plumes with Xe²⁺ ions differed from those for Xe⁺ ions. Near the thruster exit, Xe²⁺ ions were scattered less in the simulation with CEX than in the simulation without CEX, which was the opposite of the pattern for Xe⁺ particles. In the analysis of the accumulated number of particle pairs that underwent charge exchange, many fewer underwent CEX between Xe²⁺ and Xe than between Xe⁺ and Xe; thus, Xe²⁺ scattering was not remarkable. However, due to the wide high-plasma-potential zone in the simulation with CEX, doubly charged Xe²⁺ ions were expelled from the front of the axisymmetric axis, which resulted in a small number of Xe²⁺ ions in that area.

Due to the strong electric forces on the doubly charged Xe²⁺ ions, the Xe²⁺ ions were much faster than the other particles. Consequently, sputtering and erosion by impingement of Xe²⁺ ions on the satellite body were much more severe than for the other particles. Since CEX effects led to different patterns in the distribution of Xe⁺ and Xe²⁺ ions, those effects on the exhaust plume should be accurately estimated to ensure sufficient satellite longevity.