Research Article  Open Access
Numerical Study of Hall Thruster Plume and Sputtering Erosion
Abstract
Potential sputtering erosion caused by the interactions between spacecraft and plasma plume of Hall thrusters is a concern for electric propulsion. In this study, calculation model of Hall thruster’s plume and sputtering erosion is presented. The model is based on three dimensional hybrid particleincell and direct simulation Monte Carlo method (PIC/DSMC method) which is integrated with plumewall sputtering yield model. For lowenergy heavyion sputtering in Hall thruster plume, the Matsunami formula for the normal incidence sputtering yield and the Yamamura angular dependence of sputtering yield are used. The validation of the simulation model is realized through comparing plume results with the measured data. Then, SPT70’s sputtering erosion on satellite surfaces is assessed and effect of mass flow rate on sputtering erosion is analyzed.
1. Introduction
Hall thrusters have become a tempting alternative to traditional chemical propulsion systems due to the great mass saving that they provide through high specific impulses. Large specific impulses show that Hall thrusters are well suited for missions such as stationkeeping and orbit transfers. However, a major stumbling block to their widespread integration is the uncertainty about the thruster plume’s interaction with spacecraft components. For Hall thrusters, the divergence angle of particles is relatively large, leading to the possibility of direct impingements of high energy propellant ions onto the spacecraft surfaces which result in sputtering and degradation of material properties. In order to effectively assess the contamination of plasma plume and improve the performance of electric thrusters, it is crucially important to develop a computational sputtering erosion method.
Hall thruster plume is a complex rarefied flow with several species: atoms, positively charged ions, and electrons. There are two main ways for calculation of plume: semiempirical methods [1–3] and PIC/DSMC method [4–8]. Semiempirical methods are based on experimental data, and their the main assumptions are that ion flow can be approximated as originating in raylike fashion from one or two effective centers and that ion energy is purely a function of axial and radial distance. But, it is difficult to obtain the full information of the plasma plume and there is a need for abundant precise experimental data. PIC/DSMC method which is a general way to simulate plasma plume tracks every real physical process, such as particle moving, particle collision, and electric field calculation. The simulation successfully captures detailed plume structures and plume interactions.
Sputtering is commonly used for thinfilm deposition, etching, and analytical techniques. There are a lot of experiments [9] to study the sputtering yield. And lots of semiexperimental sputtering yield formulations between inlet ion and project material which based on Sigmund theory are built. For lowenergy heavyion sputtering in Hall thruster plume, the Matsunami formula [10] for the normal incidence sputtering yield and the Yamamura angular dependence of sputtering yield [11] are used in this study.
A 3D PIC/DSMC model is built for plume and sputtering yield calculation, in which plumewall interaction model for sputtering erosion calculation is integrated. The simulated SPT70’ propellant is xenon, typically with a power of 660 W, a thruster of 40 mN, and a specific impulse of 1500 s. Plume sputtering erosion on satellite’s surface is calculated. An experiment performed in a vacuum chamber [8] is used to verify simulation results by comparing with experimental data. Additionally, propellant mass flow rate effect on sputtering erosion is assessed.
2. Mathematics Models
A plasma plume is a complex rarefied flow with several species: atoms, positively charged ions, and electrons. It is firstly suggested by Oh et al. [12, 13] to combine DSMC and PIC to solve a rarefied plasma ultrasonic jet flow problem. The direct simulation Monte Carlo (DSMC) method simulates the collisions of heavy particles (ions and atoms) and the motions of particles, while PIC method is employed to calculate the electric flied and accelerate the charged particles. Electrons are treated using a fluid description, because electrons which have significantly lighter mass, can adjust their velocities more quickly than ions or atoms [14]. The algorithm of PICDSMC hybrid method is shown in Figure 1. The simulation successfully captures detailed plume structures and plume interactions. Ions which impact walls deposit their energy near the surface and a collision cascade develops which lead to targeting atoms ejected from the surface. Plumewall sputtering model is integrated in the PIC/DSMC model, in which the Matsunami formula for the normal incidence sputtering yield and the Yamamura formula for angular dependence of sputtering yield are used.
2.1. PIC/DSMC Method
2.1.1. Collision Dynamics
Collision dynamics of dilute plasma plume is involved with atomatom, atomion, electronatom, electronion, electronelectron, and ionion collisions. Here, Xe, Xe^{+}, and Xe^{2+} are taken into consideration. Collisions between these particles can be divided into two types: momentum exchange collision, which is treated as elastic collision, and charge exchange (CEX) collision. Because the mass of the electron is far less than that of ion, directly simulating its collision is lacking efficiency. Here, electrons are treated as fluid with momentum but without mass. Also collisions between ions are not modeled for calculation because of their too high collision frequency.
(1) Momentum Collision
Momentum collisions are treated as elastic collisions with variable hard sphere (VHS) [15] model. The crosssection for atomatom elastic interactions is [16]
For atomion elastic collision, the crosssection given by Dalgarno, McDowell, and Williams [17] is adopted as
where is the relative velocity between two particles; is the elastic crosssection. The elastic crosssection for interaction between an atom and a doubly charged ion is twice of that of an atom and a singly charged ion.
(2) Charge Exchange Collision
CEX collision is very important for studying Hall thruster plume and is a longrange integration with a crosssection relatively large compared to that of momentum collision. A fast ion from the main beam undergoes a CEX collision with a slow neutral particle results in a slowmoving ion and a fast neutral particle, as shown in (2.3), where “” indicates fast and “” indicates slow. Such collisions are significant for studies of spacecraft contamination. Ions with relatively high velocities tend to follow straight trajectories, while ions with slow velocities are quickly turned towards the edge of the plume. In ion thrusters, the ions in the main beam are moving in a small range of velocities which is much faster than ions created by CEX, so the CEX ions exit to the sides of the plume and form a “winglike” structure which is located at the backflow region and is the main factor of spacecraft surface pollution. “Winglike” structure is not shown in this paper, the details can be found in reference [8].
Consider the following:
Here, the CEX crosssection, , for XeXe^{+} collision measured by Pullins et al. [18] and Scott Miller et al. [19] is used as
It is assumed that there is no momentum transfer accompanying the transfer of electron(s) in CEX collisions. The crosssection for interactions between atoms and doubly charged ions is half of that of atoms and singly charged ions.
(3) No Time Counter (NTC) Method
NTC is used to determine whether a collision is happening or not [15]. In every cell, collision pairs between particles and particles could be determined by (2.5) as follows:
where is number of collision pairs; and are the numbers of and particles, respectively; is the weight of calculation zone; is time step; is volume of cell; is crosssection; is relative velocity.
Acceptreject method: collision probability is firstly calculated, and then a random number, between 0 and 1, is generated. If , the collision is accepted, that is, a collision is supposed to be happening. No collision occurs if vice versa.
2.1.2. Calculation of Electric Field
PIC method is used to trace the charged particles movements in selfinduced electric and magnetic field. Maxwell equations and Poisson equation are adopted to solve electricmagnetic field and static electric filed, respectively. Although magnetic field is significant in the accelerating channel, experimental results show that magnetic field in the plume exit decreases fast, and magnetic field in most regions is very weak. Therefore, magnetic field throughout the plume region is neglected. Based on a quasineutral assumption, the electron momentum equation (2.6) is written as (2.7), while assuming that the electrons are massless, unmagnetized, and collisionless. Further assuming that the electrons are isothermal, and that their pressure obeys ideal gas law, the Boltzmann relation, as shown in (2.8), is used to solve for the plasma potential in simulations [8, 12] as follows: where is the mass of electrons; is the number density of electrons; is the time; is the charge of an electron; is the velocity of electrons; is the electric field intensity; is the magnetic induction; is the electron pressure; is the electron temperate; is the friction term defined as (2.9); is the potential and the superscript; “*” indicates a reference value for the corresponding parameters.
Consider the following: where is the current density; is the electrical conductivity; ; is the ionelectron collision frequency; is the neutralelectron collision frequency, and these frequencies are evaluated for the xenon system using crosssections provided in [20].
In this paper, the friction term is included in the simulation though it is a weak term, which is different from the accomplished work using Boltzmann relation in [8, 12]. This treatment induces the solving process of electric field to be more complex, but it is more accurate than the Boltzmann relation. Introducing the plasma potential , a generalized Ohm’s law is obtained as follows: and with the charge continuity condition as the discretized form of (2.12) at node is where is Boltzmann constant. Using the 5point formula and assuming that the electrons are isothermal, the plasma potential can be solved by using an alternating direction implicit (ADI) scheme. Here the parameters, such as the number density of ions, accumulated onto each node of the computational domain are obtained by Ruyten’s method [21].
2.2. PlumeWall Sputtering Yield Model
Three particles are tracked in the model: Xe, Xe^{+}, and Xe^{2+}. For neutral, it is assumed that particles that impact the channel walls thermalize at the wall temperature and are diffusely scattered back into the channel. Similar to anode injection, the velocity distribution is assumed to be halfMaxwellian. Particles which exit the computational domain in the plume region are no longer tracked by the simulation. Different to the neutral treatment, ions which impact walls recombine to form neutrals and deposit their energy near the surface and a collision cascade develops which leads to target atoms ejected from the surface. Sputtering yield, the number of target atoms removed per incident particle, is dependent upon incident particle energy, incident angle, and target surface composition.
2.2.1. Sputtering Yield Model
Early in 1969, Sigmund [22] derived the formula for the sputtering yield by solving the linearized Boltzmann equation based on the assumption that the collision cascade is well developed in the infinite medium and the heat of sublimation is the surface binding energy.
Consider the following: Here, where the numerical coefficient is in units of , is the deposited energy at the surface, and is the inlet energy. is the atomic number density, is the surface binding energy of the target solid. is the coefficient of the nuclear stopping cross section.
At extreme, in cases of the lowenergy heavyion and the highenergy lightion sputtering, the original Sigmund formula should be modified. Therefore, Yamamura et al. [23] written the sputtering yield formula as an interpolation between two cases, that is, as where is the reduced energy where and are nucleus charge numbers; and are mass numbers of projectile and target; and are the reduced nuclear stopping power and the reduced electronic stopping power, respectively; is the nuclear stopping cross section; is the threshold energy; , , and are constants of the projectiletarget combination. For lowenergy heavyion sputtering where is very small, (2.17) reduces to lowenergy heavyion case, whereas for highenergy lightion sputtering where can be neglected, (2.17) reduces to highenergy lightion case.
In order to propose an empirical formula as simple as possible, Matsunami et al. [10] used in (2.17) the power approximation for [24] and the Lindhard electronic stopping power for [25]. Since of (2.17) has a dependence, it is factorized as , the factor corresponded to . The nuclear stopping cross section in the present formula is calculated by is ThomasFermi potential [10] as Then, empirical formula is finally given as is the parameter about target atoms as where and are empirical constants. The Matsunami formula is calculated for the normal incidence sputtering yield. For the dependence of the sputtering yield on the angle of incidence of the bombarding particles, Yamamura et al. [11] proposed a procedure which is based on the assumption that the angular dependence can be described by a factor to the yield at normal incidence as Here, is the incidence angle at which maximum sputtering occurs, and is a fit parameter. Yamamura gave empirical formulae for the parameters and, basing on a large number of experimental sputter data of various projectiletarget combinations. The optimum incidence angle is given by as follows: is average lattice constant. Particularly in the case of heavyion sputtering, the description reveals a singular behavior at threshold energy which is overcome by an interpolation.
2.2.2. Sheath Model
Sheaths are nonneutral regions that normally form at plasma boundaries to balance electron and ion losses. They are one of the most prominent and wellknown features of confined plasma. A crude surface model predicts the interaction between particles and surfaces in the original plume simulation models, in which the quasineutrality is assumed and resolution of the nonquasineutral sheath boundary is not accomplished. Thus, a sheath interaction model is needed to account for this sheath region [26]. The potential drop across the sheath is calculated with where is ion flux, is the electron temperature, and and are the electron charge and Boltzmann constant, respectively. Energy due to acceleration through this sheath potential is added to the particle impacting the object surface. This energy is then used to calculate the sputtering yield.
2.3. Computational Domain and Parameters
The simulation engine is a SPT70 (7 cm diameter) which works at nominal flow and current. As shown in Figure 2, thruster is placed at the middle of the surface, which has a potential influences on surface 1 and surface 2. From symmetry considerations, it is adequate to solve 1/4th sector of the satellite for high efficient calculation, as shown in the yellow part in Figure 2. The simplified computational model is shown in Figure 3, in which the calculation domain has a direction of 1.1 m, and and direction of 0.8 m, then from −1.0 m to 0.07 m in direction is the backflow region. Here the thruster inner radius, , and the outer radius, , at the thruster exit are 20.8 mm and 35 mm, respectively. The outer radius of the thruster, , is 50 mm.
The species, such as Xe^{3+} and Xe^{4+}, are not counted due to their tiny number [10]. The boundary treatments are as follows. (1) The temperature of every solid wall is assumed to be 500 K; (2) the diffusive reflection model is assumed for simulation particlewall interaction; (3) once a charged particle strikes the wall, it is neutralized and reflected back with a speed characterized by the wall temperature. Meanwhile, one electron reaches the wall in order to maintain quasineutrality; (4) the wall potential is set to 0 V (wall is metallic).
The plume parameters used for computation are listed in Tables 1 and 2. Part of the parameters is measured by a Longmuir probe. But an accurate measurement of the SPT exit plane parameters is almost impossible to achieve, so some assumptions are necessary for simulation, such as assuming the electrons are isothermal (the electron temperature is fixed at 2.0 eV). Satellite surfaces are modeled as iron. Table 3 summarizes sputtering yield parameters of iron by incident xenon particles.



3. Results and Discussion
3.1. Plume Parameter Validation
The angular profiles of the electron temperature and density of plume flow field are measured by the double probe when the thruster is working under a rated propellant mass flow rate of 2.68 mg/s and a discharge voltage of 300 V [8]. Here, is defined as the distance from the center of thruster exit plane.
Figure 4 shows the simulation results and experimental results. It is shown that the variations of both experimental and computational results are of similar shape, and simulation results overestimate both number density and current density at the center of plume while underestimate these parameters at the wings. It is evident that the simulation results are reasonable and referential at most of the measuring points though at some places the errors between the measured data and the computational results are relatively large. The errors between the experimental and the computational results are observable at some locations due to the following reasons: (1) experimental error consists of at least four parts in measuring: misalignment of probe, contamination of probe (surface), effect of sheath, and effect of chamber wall, which results in a definite measuring error; (2) simplifications and assumptions for calculation, such as the assumption of isothermal electrons is not accurate, especially near the exit, the simplified distribution radial velocities of injected particles at the thruster exit, and the treatment of backpressure. One thing that should be noticed is that the effect of chamber wall is seen at large angles since the experimental tank is not large enough.
(a) , 50 cm
(b) , 100 cm
3.2. Results
The sensitive surfaces that should be paid attention to are surfaces 1 and surface 2. They are main functional surface which are exposed to the SPT70 plasma plume. Figures 5 and 6 are the total particles number density contour of  slice through the surface’s center line when the thruster is working under a rated propellant mass flow rate of 2.68 mg/s and a discharge voltage of 300 V. It is seen that surface 1’s maximum number density is 1 × 10^{16} which is at the top of the surface and that the number density value reduces gradually. For surface 2, number density value also reduces, but its maximum is 1.5 × 10^{15}. Sputtering erosion is the cumulative effect of each incident particle, so surface 1 may have a high influence. The calculation result of this two surfaces’ sputtering depth is given in Figures 7 and 8. These results show that surface 1 and surface 2’s max sputtered depth is 1.8 μm and 0.5 μm per 200 hour operation time, respectively. As expected, surface 1 has a higher sputtering depth than surface 2. In the condition of fixed ion and surface material, the differences between the two faces are due to the following reasons: (1) the distance between target surface and the engine plume: ion number density and ion energy are decaying with the distance increasing; (2) incident angle: it is the angle between surface and engine plume. Usually, when incident angle is between 50 and 70, sputtering yield reaches a maximum.
Fix the discharge voltage at 300 V, the effect of mass flow rate is studied. Here, three mass flow rates are 3.12 mg/s, 2.68 mg/s, and 1.56 mg/s. The distance between sensitive surface and outlet of the thruster is about 25 cm. As shown in Figures 9 and 10, both electron number density and ion current density at 25 cm are presented. Electron number density and ion current density are increasing with the increase of flow rate and their variations are of the same shape, though the ratios of corresponding parameters are not strictly equivalent to those of flow rates. In fact, the increment or the decrement rate of both number and current density seems to be less than those of flow rates. For example, the electron number densities at the centerline according to three different mass flow rate levels, 3.12 mg/s, 2.68 mg/s, and 1.56 mg/s, are 4.17 × 10^{16} m^{−3}, 3.58 × 10^{16} m^{−3}, and 2.09 × 10^{16} m^{−3}, respectively. It is also shown that electron and current density’s decrement rate is larger than divergence angle rate. Max sputtering depth on surface 1 and surface 2 is shown in Figure 11. It can be seen that sputtering depth increased with flow rate, but the rate of change is different in the two faces. Since the two faces are in different plume divergence angle range in which surface 2’s divergence angle is greater than surface 1, increment rate of max sputtering depth on surface 1 is obviously larger than that on surface 2.
4. Conclusions
Threedimensional PIC/DSMC computational modeling is built for Hall thruster plume and sputtering erosion calculation. The Matsunami formula for the normal incidence sputtering yield, the method of Yamamura for angular dependence of sputtering yield and sheath model are integrated to PIC/DSMC program for sputtering calculation. Evaluation of SPT70 Plume and sputtering erosion on satellite surfaces is calculated. Computational plume results agree well with the experimental data. The maximum of sputtering depth on satellite surfaces is 1.8 μm/200 h. In addition, effect of mass flow rate on plume parameter and sputtering erosion is considered. It is shown that the increment or the decrement rate of both number and current density seems to be less than those of flow rates and sputtering depth increasing rate with mass flow rate on surface 1 is obviously larger than that on surface 2.
Acknowledgments
This paper is supported by the Shanghai Nature Science Fund (no. 12ZR1414700) and Innovation Fund of Shanghai Jiao Tong University.
References
 P. M. Burgasov, G. V. Grigorian, A. A. Izmaelov et al., “Some results of complex work concerning problems of electric rocket thruster integration with a spacecraft and its subsystems,” in Proceedings of the 2nd European Spacecraft Propulsion Conference, Noordwijk, The Netherlands, 1997. View at: Google Scholar
 D. Borie, V. Perrin, S. Khartov, and A. Nadiradze, “The I.S.P. software: calculation of the SPT jet influence,” in Proceedings of the 2nd European Spacecraft Propulsion Conference, Noordwijk, The Netherlands, 1997. View at: Google Scholar
 A. B. Nadiradze and S. A. Kharov, “A 3D model calculating sputtering and depositing processes under electric propulsion thruster testing in a vacuum chamber,” in Proceedings of the 29th International Electric Propulsion Conference, Princeton University, 2005. View at: Google Scholar
 E. Sommier, M. K. Allis, N. Gascon, and M. A. Cappelli, “Wall erosion in 2D Hall Thruster simulations,” in Proceedings of the 42nd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, pp. 3335–3344, Sacramento, Calif, USA, July 2006. View at: Google Scholar
 S. Kay, M.S. Manuel, and B. Oleg, “Integration of a sputtering model into a full PIC hall thruster simulation,” in Proceedings of the 39th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, Huntsville, Ala, USA, 2003. View at: Google Scholar
 R. Subrata and B. P. Pandey, “Modeling the effect of plasmawall interaction in a hall thruster,” in Proceedings of the 41st Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, 2003. View at: Google Scholar
 A. Bernard and M.S. Manuel, “Computation modeling of plasma plume in a vacuum tank,” in Proceedings of the 26th International Electric Propulsion Conference, Kitakyushu, Japan, 1999. View at: Google Scholar
 Q. Zhong, W. Pingyang, D. Zhaohui, and K. Xiaolu, “Study of plume characteristics of a stationary plasma thruster,” Plasma Science and Technology, vol. 10, no. 5, pp. 612–618, 2008. View at: Publisher Site  Google Scholar
 Y. Yamamura and H. Tawara, “Angular dependence of sputtering yelds of monatomic solids,” Institute of Plasma Physics IPPJAM26, Nagoya University, 1996. View at: Google Scholar
 N. Matsunami, Y. Yamamura, Y. Itikawa et al., “Energy dependence of the ioninduced sputtering yields of monatomic solids,” Atomic Data and Nuclear Data Tables, vol. 31, no. 1, pp. 1–80, 1984. View at: Google Scholar
 Y. Yamamura, Y. Itikawa, and N. Itoh, “Angular dependence of sputtering yields of monatomic solids,” Institute of Plasma Physics IPPJAM26, Nagoya University, 1983. View at: Google Scholar
 D. Y. Oh, Computational Modeling of Expanding Plasma Plumes in Space Using a PICDSMC Algorithm, Massachusetts Institute of Technology, Cambridge, Mass, USA, 1997.
 D. Y. Oh, D. E. Hastings, C. M. Marrese, J. M. Haas, and A. D. Gallimore, “Modeling of stationary plasma thruster100 thruster plumes and implications for satellite design,” Journal of Propulsion and Power, vol. 15, no. 2, pp. 345–357, 1999. View at: Google Scholar
 Y. Choi, M. Keidar, and I. D. Boyd, “Particle simulation of plume flows from an anodelayer hall thruster,” Journal of Propulsion and Power, vol. 24, no. 3, pp. 554–561, 2008. View at: Publisher Site  Google Scholar
 G. A. Bird, Molecular Gas Dynamics and the Direct Simulation of Gas Flows, Oxford University Press, New York, NY, USA, 1994.
 K. Koura and H. Matsumoto, “Variable soft sphere molecular model for air species,” Physics of Fluids A, vol. 4, no. 5, pp. 1083–1085, 1992. View at: Google Scholar
 A. Dalgarno, R. C. McDowell, and A. Williams, “The mobilities of ions in unlike gases,” Philosophical Transactions of the Royal Society A, vol. 250, pp. 411–425, 1958. View at: Publisher Site  Google Scholar
 S. H. Pullins, Y. Chiu, D. J. Levandier, and R. A. Dressler, “Ion dynamics in Hall effect and ion thrustersXe^{+} + Xe symmetric charge transfer,” in Proceedings of the 38th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, 2000. View at: Google Scholar
 J. Scott Miller, S. H. Pullins, D. J. Levandier, Y. H. Chiu, and R. A. Dressler, “Xenon charge exchange cross sections for electrostatic thruster models,” Journal of Applied Physics, vol. 91, no. 3, p. 984, 2002. View at: Publisher Site  Google Scholar
 M. Mitchner and C. H. Kurger, Partially Ionized Gases, Wiley, New York, NY, USA, 1973.
 W. M. Ruyten, “Densityconserving shape factors for particle simulations in cylindrical and spherical coordinates,” Journal of Computational Physics, vol. 105, no. 2, pp. 224–232, 1993. View at: Publisher Site  Google Scholar
 P. Sigmund, “Theory of sputtering. I. Sputtering yield of amorphous and polycrystalline targets,” Physical Review Online Archive, vol. 184, pp. 383–416, 1969. View at: Publisher Site  Google Scholar
 Y. Yamamura, N. Matsunami, and N. Itoh, “Theoretical studies on an empirical formula for sputtering yield at normal incidence,” Radiation Effects Letters, vol. 71, no. 12, pp. 65–86, 1983. View at: Google Scholar
 J. Lindhard, V. Nielson, M. Scharff, and P. V. Thomsen, “Integral equations governing radiation effects,” Kongelige Danske Videnskabernes Selskab Biologiske Skrifter, vol. 33, no. 10, pp. 1–42, 1963. View at: Google Scholar
 J. Lindhard and M. Scharff, “Energy dissipation by ions in the kev region,” Physical Review, vol. 124, no. 1, pp. 128–130, 1961. View at: Publisher Site  Google Scholar
 C. Shannon, Computational modeling of hall thruster plasma plume in a vacuum tank [M.S. thesis], Massachusetts Institute of Technology, Boston, Mass, USA, 2003.
Copyright
Copyright © 2012 Li Yan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.