# Numerical Simulation of Interaction between Hall Thruster CEX Ions and SMART-1 Spacecraft

Authors: Kang Shan; Yuchuan Chu; Qingyu Li; Liang Zheng; Yong Cao
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

The interaction between the plume of Hall thruster and the surface of the SMART-1 spacecraft is investigated by developing a three-dimensional IFE-PIC-MCC code, with the emphasis on the
effect of the disturbance force and thermal loading caused by charge exchange ions (CEX) impingement on the surface of the spacecraft.  The parameters such as heat flux and forces of CEX ions which impinge on SMART-1 and solar arrays are obtained.  The disturbance force of CEX ions to the spacecraft is calculated for different divergence angles and different solar array rotation cases.  The simulation results show that the disturbance force and heat flux on spacecraft change very little as the divergence angle changes. The effect of maximum disturbance force can be neglected since it is so small comparing with the nominal value of the main thrust. Solar arrays receive the least thermal heating from the CEX ions when the beam ions flow is perpendicular to the solar array plane.


---

## Body

## 1. Introduction

The motion of the satellite is usually controlled by the ejection of the plume from the thruster into the space. Then, the interaction between the plume and the spacecraft surface may cause undesirable effects such as causing the disturbance force and thermal loads and contaminating sensitive equipment and sensors. The disturbance force can be a fraction of the total thrust, while thermal loads on the surface of spacecraft body result in the heating of the surface and affect the working status of electronic components which can only function properly in a range of temperatures. So the accurate modeling and predictions of these effects are very crucial to the design of a satellite [1–3].The interaction between the exhausted plume of thrusters and the satellite components has been studied by some researchers for both chemical thrusters and electric propulsions thrusters [2, 4]. Park et al. [2] used three-dimensional discrete simulation Monte Carlo (DSMC) to investigate the interaction of the chemical thruster (a 4.45N MRE-1 monopropellant hydrazine liquid rocket engine) plume with satellite components in KOMPSAT-II. The results showed a negligible disturbance force/torque and thermal loading compared with its nominal thrust/torque and solar heating. Xiao et al. [5] analyzed molecules adsorption and transmission on the surface of satellite by using numerical simulation and ground experiment method. Also, the motion of plume pollutants which leads to performance degradation of satellite key functional surfaces (optical systems, solar panels, thermal control object surface, etc.) is calculated. Different from the plume of the chemical thrusters, the plume of electric propulsion is plasma which consists of a large number of ions and electrons except neutral atoms. In addition, charge exchange collisions will occur between the high-speed ions and neutral atoms which result in the generation of low-speed CEX ions that have significant impact on the plume characteristics. These charged particles are affected heavily by their self-consistent electric fields. Therefore, compared with the plume of the chemical thrusters, the plume of the electric propulsion thruster is different in not only the ingredients of the plume but also the flow characteristics of the plume. A number of simulation models or numerical methods have been developed for the plume in the electric propulsion thruster to investigate the interaction with the spacecraft surface and the results from these numerical models were verified through the comparison with the experiment data in recent years [6–13]. Yan et al. [10] used particle in cell (PIC) code with DSMC techniques to model Hall thruster plume and sputtering erosion on SPT-70. Kafafy and Cao [11] investigated plume effects from indirect plume impingement on formation flying satellites using ion propulsion by developing an immersed-finite-element PIC (IFE-PIC) algorithm on parallel computers. Tajmar et al. [12] developed a hybrid PIC code with Monte Carlo collision (MCC) to study spacecraft-environment interaction. Boyd [13] studied the ion current density profile and ion energy distribution by using a detailed particle-fluid PIC-DSMC model, which was compared with the experimental measurements taken in space. However, most of them were concerned about CEX ions sputtering erosion or the accuracy of simulation model, but few of them were concerned about the disturbance force and thermal loading on the spacecraft which are caused by the impinging of the backflow CEX ions on the surface of the spacecraft.Therefore, in this paper, the study of the force and thermal loads on the spacecraft caused by backflow CEX ions is performed by using a three-dimensional IFE-PIC-MCC code. The PIC-MCC [12] code is used to simulate the generation and movement of CEX ions. The DSMC method [14] is applied to model neutral atoms. Electric field in the plume is obtained by solving Poisson’s equation which is calculated by IFE-PIC [15, 16] method which is designed to handle complex boundary conditions accurately while maintaining the computational speed of the standard PIC code. The code is then applied to the numerical simulations of the SMART-1 spacecraft which had traveled to the moon using a PPS-1350 Hall thruster with the maximum thrust of 70 mN [17].Section2 describes the interaction model between the plume of Hall thruster and the surface of the spacecraft. The numerical method is presented in Section 3. The simulation results are then shown in Section 4 and some discussions on these results are carried out. Finally, the summary and conclusions are presented at the end of this paper.


## 2. SMART-1 Spacecraft-Plume Interactions Model

The geometry and dimensions of SMART-1 spacecraft model are illustrated in Figure1. The main body of SMART-1 can be considered as the cubic shape with the dimensions of l×w×h=1100mm×1100mm×900mm. In this model, the Hall thruster is simplified as a cylinder with the diameter of 100 mm and the height of 50 mm; two thin rectangles with 5400 mm length and 1000 mm width are utilized to represent the solar arrays which can rotate around the satellite.Figure 1
The geometry and dimensions of SMART-1 (in mm).The PPS-1350 Hall thruster emits an ion beam out of a ring-shaped anode with a divergence angleθ=45°. Typical operating parameters of PPS-1350 Hall thruster are summarized in Table 1 [13]. At the exit of the thruster, the densities and temperatures are assumed to be radically uniform, and the velocity vectors vary uniformly from -45° at the lower edge of the channel exit to +45° at the upper edge. The temperatures of the electrons and neutrals are assumed to be 2 eV and 750 K, respectively. The heavy particle densities and velocities are subsequently obtained from the mass flow rate and integrated ion current. The plume in the PPS-1350 Hall thruster assumingly consists of the following components:(i)
propellant beam ions,(ii)
unionized propellant neutrals,(iii)
slow propellant ions created by CEX reaction collisions,(iv)
neutralising electrons.Table 1
Operating parameters of PPS-1350 Hall thruster.
Parameter
Value
Thrust
70 mN
VoltageUacc
350 V
Ion currentIi
3.8 A
Mass flow rate
4.2 mg/s
Specific impulse
164 s
Total efficiencyη
51%
Outer insulator diameterro
100 mm
Inner insulator diameterri
56 mm
Neutral temperature
750 K
Electron temperatureTe
2 eVThe plasma in the beam of the Hall thruster typically has a number density of more than1015 m−3. The widely used PIC approach is applied in this research. The so-called super particle representing a number of 109 real particles is introduced here. As a result, the total number of simulated entities of Xe and Xe+ is kept below 107; therefore the simulations can be carried out on standard workstations [17].X
e
+ particles are given a constant velocity at the thruster exit as
(1)vi=2eUaccmi,
where Uacc is accelerating voltage, mi is xenon ion mass, and vi denotes the ion velocity.Initially the ions are uniformly distributed across the thruster opening between the inner radiusri and the outer radius ro. The direction of the velocity vi is determined in spherical coordinates (vi;α1;α2) [17], where α1 is randomly chosen in 0<α1<2π and α2 is varied in accordance with the radial position between the inner and outer beam-spreading angles αi and αo by
(2)α2r=αo-αiro-rir-ri+αi,
where ro and ri are the outer and inner radii at the exit of the thruster, respectively.During the simulations, the number of ejectedXe+ beam ions by the Hall thruster at each time-step is obtained by the ion current:
(3)Δni=IiΔt,
where Ii denotes ion current, Δni is number of ions generated at each time step, and Δt represents the time of each step.The electrons are assumed to be collisionless, currentless, isothermal, and unmagnetized. The plasma is quasineutral and electrons satisfy the Boltzmann [12] relation as
(4)ne=ne0exp⁡eϕ-ϕ0kTe0,
where ne is electron number density, ϕ is potential, and ne0, ϕ0 represent the reference parameters.The electric field surrounding the spacecraft is solved from Poisson equation(5)-∇·ε∇ϕ=eni-ne,
where ε denotes permittivity of vacuum, ne is electron number density, ni represents ion number density, and ϕ is the potential.Charged particles in the simulation area are accelerated due to the electric field. Their movement can be determined by the integration of the equation of motion, which can be described with Newton’s second law:(6)ddtmv=F=qE,v=dxdt.NeutralXe propellant atoms are introduced with no drift velocity but the thermal velocity related to the temperature of the neutral Xe propellant atom. As neutral Xe atoms move, the DSMC method is applied to simulate collisions of neutral Xe propellant atoms since neutral propellant can be considered as a rarefied gas flow. In this method, each unit cell, containing many particles, has the dimensions on the order of a mean free path. Pairs of these particles are then randomly selected and a collision probability is evaluated which is proportional to the product of the relative velocity and cross section of the collision for each pair. The probability is compared with a random number to determine if that collision occurs.Let us consider a single simulation cell containing a number of simulation macroparticles. Bird [14] gives the probability of a collision for a macroparticle in this cell:
(7)P=FNσTcrΔtVc,
where FN is the macroparticle weight, σT is the total collision cross section, cr is the relative velocity, Δt is the simulation time step, and Vc is the cell volume.One order to check for collisions would be to iterate over allN particles and compute probability with all remaining particles. This would result in N(N-1)/2~(N2)/2 pairs. For a large number of particles, this method clearly becomes very computationally inefficient. Bird’s no time counter (NTC) method was designed to provide help with this issue. It allows us to estimate ahead of time the maximum number of pairs that need to be checked. The maximum collision probability is
(8)Pmax⁡=FNσTcrmax⁡ΔtVc,
where σTcrmax⁡ is a parameter chosen ahead of time using approximate predictions of cross section and velocity. As Bird points out in his book, the actual value is not all that important since it ends up getting cancelled out. The number of pairs to check is then
(9)12NN-FNσTcrmax⁡ΔtVc.
Here the second N- is the average particle count. The average is used to reduce the statistical time-step to time-step oscillations. For each pair, we then compute the probability as follows:
(10)P=σTcrσTcrmax⁡.
This value is compared to a random number and the collision occurs if P>R.For atom-atom elastic collisions, the variable hard sphere collision model is employed. For xenon, the collision cross section is(11)σELXe,Xe=2.12×10-18g2ωm2,
where g is the relative velocity and ω denotes the viscosity temperature exponent for xenon which has the value of 0.12 [14].Charge-exchange collisions occur between the fast beamXe+ ions and the slow neutral Xe atoms. After these collisions, the fast beam Xe+ ions turn to be neutral and keep the fast speed; and the slow neutral Xe atoms change into ions which are accelerated due to the effect of the electric field. These changes can be expressed as follows:
(12)Xeslow+Xefast+⟶Xeslow++Xefast.To model the CEX process, we use the MCC method. For eachXe+ ion we calculate the collision probability. This probability is then compared to a random number. If the probability is larger than the random number, the collision occurs and a process-specific collision handler is called. The collision probability is given by Birdsall [18] as
(13)P=1-exp⁡-νΔt=1-exp⁡-nnσgΔt.nn is the density of the target gas at the location of the particle, σ is the collision cross section, g is the relative velocity, and Δt is the time difference between collision checks. This time-step may correspond to the time-step used to move the ions in the PIC method.In this research, the transfer of only one electron between an atom and an ion is considered during the charge-exchange process. For singly charged ions, the following cross section measured by Pullins et al. and Miller et al. [19, 20] is used:
(14)σCEXXe,Xe+=-23.30log⁡10g+142.21×0.8423×10-20m2,
where g is the relative velocity between the beam ion and neutral atom.Boundary conditions for computing of the Hall thruster plume are presented as follows. Both field and particle boundary conditions are required at the outer boundary of the computational domain. The field condition is simply set as the electric field normal to the boundary edges equal to zero. The particle boundary condition is set to be that any particle crossing the boundary is simply removed and will no longer be calculated in the model. The solid surfaces of the Hall thruster and the satellite are also considered during the simulations. Along these surfaces, the potential of surfaces is set to be at a certain given value (−2 V). Any ions colliding onto the surfaces are absorbed by the surfaces.


## 3. Numerical Method

A particle simulation code using the IFE-PIC-MCC algorithm is developed to solve the plume interaction problem. The IFE-PIC-MCC code solves the generating of CEX ions and their trajectories and the electric field surrounding the spacecraft.The simulation is run in three phases: modeling beamXe+ ions, neutral Xe, and CEX ions, respectively. The reason that these different species particles can be simulated asynchronously is as follows.(1)
BeamXe+ ions have so high energy that they are hardly affected by self-consistent electric field.(2)
NeutralXe atoms do not get charged so that they are not affected by electric field.(3)
CEX ions are in small numbers in plasma plume and their generation affects the plume character very little.In the first phase, we run the PIC code to generate and trace beamXe+ ions and neutral Xe until steady-state trajectories are obtained. In this phase, we use a simplified approach to model beam Xe+ ions. Those beam Xe+ ions are accelerated in the Hall thruster channel with the velocities on the order of 10 km/s. So most high-energy beam Xe+ ions are unaffected by self-consistent electric field and follow straight trajectories. The DSMC method is used in the computation of Xe-Xe collisions. After this phase, the distribution of beam Xe+ ions and neutral Xe is gotten in the computation region, respectively.In the second phase, the IFE-PIC-MCC code is used to trace CEX ions until steady-state trajectories are obtained. CEX ions are generated from the collision of steady-state beamXe+ ions and neutral Xe by MCC at each time step. The initial velocity of CEX ions is given a uniform velocity distribution with a temperature corresponding to that of the neutral propellant. That is to say, initial CEX ions only have thermal velocity with no drift velocity which is relatively small and is affected heavily by electric field. Using IFE-PIC, the electric field is calculated which pushes the movement of CEX ions at each time step. Those processes are cycled until steady-state CEX ions trajectories are obtained.In the final phase, CEX ions are collected which impinge on SAMRT-1 and solar arrays; then all those CEX ions are used to calculate the flux density of CEX ions, the energy distribution, and pressure distribution on SMART-1 and solar arrays.


## 4. Simulations Results

Due to the symmetry of SMART-1 (as illustrated in Figure1), a half-symmetry model is employed to simulate the interaction between the plume and the surface of SMART-1, as shown in Figure 2. The position of SMART-1 in the simulation region is also shown in Figure 2. The simulation domain has a size of 4.0m×7.5m×4.0m with 81×151×81 grid points. During the simulations, approximately 8 million simulation particles are employed. The time-step is on the order of 10-6 s. The simulation reaches a steady state for CEX ions after about 800 iterations, and solutions are then averaged over further 200 iterations. The total computation time is about 24 h on an HP computer server.Figure 2
The dimensions of the simulation region.Figure3 illustrates the distribution of the potential under the conditions of θ=45° and γ=0° (here γ is the angle between the normal of the solar array plane and the center line of the beam flow). As shown in Figure 3, the backflow of the CEX ions results in the fact that the isopotential surface expands toward a direction opposite to the ion beam movement. It makes the pattern of the distribution look like a mushroom.Figure 3
Plasma potential distribution (atθ=45°, γ=0°).
(a)
(b)Figures4 and 5 present the pressure and energy distribution of CEX ions on solar arrays, respectively. As solar arrays can rotate around the axis to change their direction, the effect of the angle between the normal of solar array plane and the beam flow c enter line is investigated. As illustrated in Figures 4 and 5, three different angles are chosen as γ=0°, 45°, and 90° and their effects are presented. It is clear to see that the solar arrays receive more CEX ions and bear much more pressure at γ=0° than at γ=90°.Figure 4
Pressure distribution of CEX ions on the solar array of SMART-1 (atθ=45°) (in N/m2).
(a)
γ
=
0
°
(b)
γ
=
4
5
°
(c)
γ
=
9
0
°Figure 5
Energy distribution of CEX ions on the solar array of SMART-1 (atθ=45°) (in W/m2).
(a)
γ
=
0
°
(b)
γ
=
4
5
°
(c)
γ
=
9
0
°Figures6 and 7 present the energy and pressure distribution about CEX ions impinging on the surface of SMART-1 main body. CEX impingement on main body and CEX ions distribution on main body are almost the same for different solar array rotation angles, since different solar array positions can hardly affect the electric field distribution in the vicinity of the main body. Also, it can be observed that the less backflow CEX ions on the main body can be obtained along with the increasing of the distance to the location of the Hall thruster. The maximum energy of CEX ions reaching onto the main body of the spacecraft is 17.056 W/m2.Figure 6
CEX ions energy distribution on the surface of SMART-1 main body (atθ=45°) (in W/m2).Figure 7
CEX ions pressure distribution on the surface of SMART-1 main body (atθ=45°) (in N/m2).Tables2 and 3 list the values of the disturbance normal force and thermal loading on the SMART-1 spacecraft for different solar array rotation angle cases, which are calculated based on the backflow CEX ions energy and pressure distribution. It can be seen that the values of the average energy and disturbance force on the main body of the spacecraft are the same under different solar array rotation angle, while the values of the average energy and disturbance force on solar arrays decrease as the solar array rotation angle γ increases. The maximum value of disturbance force onto spacecraft can reach 0.7 mN when γ=0°. The tangential force acting on the main body is symmetric due to the symmetric distribution of the pressure on the main body (as shown in Figure 7). Although tangential force on the half main body can reach a value of 0.1512 mN, the total tangential force on main body is counteracted and can be neglected. As the rotation angle of the solar array increases, the tangential force in the direction of array plane length direction decreases due to less CEX impingement on solar arrays. But the total effect of the tangential force on the solar arrays is limited, since the total tangential force on both solar arrays is mostly counteracted with each other.Table 2
Parameter values on main body for different solar array rotation angle cases (atθ=45°).
Solar array  rotation angleγ (in deg)
Average energy(in W/m2)
Disturbance normal force (in mN)
0
2.46513
0.5652
45
2.46513
0.5652
90
2.46513
0.5652Table 3
Parameter values on one of solar arrays for different solar array rotation angle cases (atθ=45°).
Solar array rotation angleγ (in deg)
Average energy (in W/m2)
Disturbance normal force (in mN)
Tangential force in array width direction (in mN)
Tangential force in array length direction (in mN)
0
0.16692
0.09449
0
0.2071
45
0.12463
0.05512
0.0354
0.1562
90
0.02903
0.01105
0.00732
0.0748Figure8 presents the distribution of the plasma potential for different divergence angle. It can be seen in Figure 8 that the distributions of the plasma potential are quite different under the condition of the same rotating angle but different divergence angle. The divergence angle also affects the CEX ions impinging on the spacecraft. The effects of three different divergence angles are studied for the same mass flow rate (4.2 mg/s). The solar arrays are set in the position that solar arrays are perpendicular to ion beam (γ=0°), and the parameters of CEX ions impingement are recorded for the three divergence angles: θ=15°, 30°, and 45°.Figure 8
Distribution of plasma potential for different divergence angles.
(a)
Atθ=15°, γ=0°
(b)
Atθ=45°, γ=0°As listed in Table4, the average energy on both the main body and solar arrays does not change much when the divergence angle changes. As divergence angle increases, the disturbance force on the main body increases very little. The disturbance force to solar arrays is relatively small comparing to the disturbance force on the main body, since solar arrays are far from the flux beam flow and receive less CEX ions. The maximum disturbance force is less than 0.1% of the nominal value of the main thrust. Thus, the disturbance force has little influence on the dynamic system of the satellite, as well as the propellant budget.Table 4
Parameter values for different divergence angle cases (atγ=0°).
Divergence angleθ (in deg)
Average energy on main body (in W/m2)
Average energy on solar arrays (in W/m2)
Disturbance force to main body (in mN)
Disturbance force to solar arrays (in mN)
15
2.43187
0.18829
0.54366
0.17317
30
2.42073
0.19279
0.54532
0.19483
45
2.46513
0.16692
0.56523
0.18897


## 5. Summary and Conclusion

In summary, the interaction between the plume of the Hall thruster and the SMART-1 spacecraft has been investigated using a three-dimensional IFE-PIC-MCC code. In this research, the backflow CEX ions impingement on the SMART-1 is investigated, and the effects of the disturbance force and the thermal loading due to the impingement are studied. The effects of different divergence angles and solar array rotation angles on the disturbance force are considered. Simulation results revealed that the maximum disturbance normal force can only reach a value of 0.7 mN, which is relatively small comparing with the main thrust (70 mN). As the solar array rotation angle increases, both the disturbance force on the spacecraft and the average energy on solar arrays decrease substantially. With the increment of the divergence angle, both the disturbance force and the average energy on the surface of spacecraft are slightly changed under the same solar array rotation angle conditions, whereas the disturbance force in the main body has a slight rise. But the maximum disturbance force is less than 0.1% of the nominal value of the main thrust, which means the effect of the disturbance force can be ignored.



---

*Source: 418493-2015-06-23.xml*