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 particle-in-cell and direct simulation Monte Carlo method (PIC/DSMC method) which is integrated with plume-wall sputtering yield model. For low-energy heavy-ion 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, SPT-70’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 station-keeping 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]. Semi-empirical methods are based on experimental data, and their the main assumptions are that ion flow can be approximated as originating in ray-like 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 thin-film 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 low-energy heavy-ion 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 plume-wall interaction model for sputtering erosion calculation is integrated. The simulated SPT-70’ 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 PIC-DSMC 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. Plume-wall 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 atom-atom, atom-ion, electron-atom, electron-ion, electron-electron, and ion-ion collisions. Here, Xe, Xe+, and Xe2+ 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 cross-section for atom-atom elastic interactions is [16]
For atom-ion elastic collision, the cross-section given by Dalgarno, McDowell, and Williams [17] is adopted as
where is the relative velocity between two particles; is the elastic cross-section. The elastic cross-section 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 long-range integration with a cross-section 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 slow-moving 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 “wing-like” structure which is located at the backflow region and is the main factor of spacecraft surface pollution. “Wing-like” structure is not shown in this paper, the details can be found in reference [8].
Consider the following:
Here, the CEX cross-section, , for Xe-Xe+ 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 cross-section 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 cross-section; is relative velocity.
Accept-reject 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 self-induced electric and magnetic field. Maxwell equations and Poisson equation are adopted to solve electric-magnetic 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 ion-electron collision frequency; is the neutral-electron collision frequency, and these frequencies are evaluated for the xenon system using cross-sections 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 5-point 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. Plume-Wall Sputtering Yield Model
Three particles are tracked in the model: Xe, Xe+, and Xe2+. 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 half-Maxwellian. 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 low-energy heavy-ion and the high-energy light-ion 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 projectile-target combination. For low-energy heavy-ion sputtering where is very small, (2.17) reduces to low-energy heavy-ion case, whereas for high-energy light-ion sputtering where can be neglected, (2.17) reduces to high-energy light-ion 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 Thomas-Fermi 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 projectile-target combinations. The optimum incidence angle is given by as follows: is average lattice constant. Particularly in the case of heavy-ion sputtering, the description reveals a singular behavior at threshold energy which is overcome by an interpolation.
2.2.2. Sheath Model
Sheaths are non-neutral regions that normally form at plasma boundaries to balance electron and ion losses. They are one of the most prominent and well-known 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 SPT-70 (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 Xe3+ and Xe4+, 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 particle-wall 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 quasi-neutrality; (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 SPT-70 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 × 1016 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 × 1015. 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 × 1016 m−3, 3.58 × 1016 m−3, and 2.09 × 1016 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
Three-dimensional 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 SPT-70 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.