Journal of Applied Mathematics

Volume 2012 (2012), Article ID 327021, 16 pages

doi:10.1155/2012/327021

Research Article

Numerical Study of Hall Thruster Plume and Sputtering Erosion

Li Yan,¹ Ping-Yang Wang,¹ Yang-Hua Ou,¹ and Xiao-Lu Kang²

¹School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai 200240, China

²Department of Electric Propulsion, Shanghai Institute of Space Propulsion, Shanghai 200233, China

Received 14 August 2012; Accepted 28 October 2012

Academic Editor: Yansheng Liu

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.

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.

Figure 1: One time step of the PIC-DSMC algorithm.

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 Xe²⁺ 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 Xe²⁺. 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.

Figure 2: Spacecraft geometries.

Figure 3: Simulated geometry.

The species, such as Xe³⁺ and Xe⁴⁺, 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.

Table 1: Exit parameters of plasma plume.

Species
Xe
Xe⁺
Xe²⁺

Velocity, m/s
325
15000
21000

Temperature, K
1000
35000
35000

Number density, 10¹⁷ m⁻³
12.40
2.72
0.30

Table 2: Computational parameters.

Parameter
Value

Electron temperature, eV
2.0

Reference electron number density, m⁻³
6.5 × 10¹⁵

Reference potential, V
0

Time step, s
1.0 × 10⁻⁷

Table 3: Sputtering parameters.

Parameter
Value

Sublimation energy, eV
4.28

Screening radius
0.0976

Absolute threshold energy, eV
46.9

Pseudo threshold energy, eV
70

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

Figure 4: Comparisons between computational results and experimental data.

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 × 10¹⁶ 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¹⁵. 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.

Figure 5: Contours of number density for slice through surface 1.

Figure 6: Contours of number density for slice through surface 2.

Figure 7: Sputtering depth for surface 1 (μm/200 h).

Figure 8: Sputtering depth for surface 2 (μm/200 h).

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¹⁶ m⁻³, 3.58 × 10¹⁶ m⁻³, and 2.09 × 10¹⁶ m⁻³, 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.

Figure 9: Electron density at = 25 cm for different mass flow rates.

Figure 10: Current density at = 25 cm for different mass flow rates.

Figure 11: Max sputtering depth on surfaces.

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.

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.

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.

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.

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.

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.

R. Subrata and B. P. Pandey, “Modeling the effect of plasma-wall interaction in a hall thruster,” in Proceedings of the 41st Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, 2003.

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.

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.

Y. Yamamura and H. Tawara, “Angular dependence of sputtering yelds of monatomic solids,” Institute of Plasma Physics IPPJ-AM-26, Nagoya University, 1996.

N. Matsunami, Y. Yamamura, Y. Itikawa et al., “Energy dependence of the ion-induced sputtering yields of monatomic solids,” Atomic Data and Nuclear Data Tables, vol. 31, no. 1, pp. 1–80, 1984.

Y. Yamamura, Y. Itikawa, and N. Itoh, “Angular dependence of sputtering yields of monatomic solids,” Institute of Plasma Physics IPPJ-AM-26, Nagoya University, 1983.

D. Y. Oh, Computational Modeling of Expanding Plasma Plumes in Space Using a PIC-DSMC 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 thruster-100 thruster plumes and implications for satellite design,” Journal of Propulsion and Power, vol. 15, no. 2, pp. 345–357, 1999.

Y. Choi, M. Keidar, and I. D. Boyd, “Particle simulation of plume flows from an anode-layer hall thruster,” Journal of Propulsion and Power, vol. 24, no. 3, pp. 554–561, 2008.

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.

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.

S. H. Pullins, Y. Chiu, D. J. Levandier, and R. A. Dressler, “Ion dynamics in Hall effect and ion thrusters-Xe⁺ + Xe symmetric charge transfer,” in Proceedings of the 38th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nev, USA, 2000.

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.

M. Mitchner and C. H. Kurger, Partially Ionized Gases, Wiley, New York, NY, USA, 1973.

W. M. Ruyten, “Density-conserving shape factors for particle simulations in cylindrical and spherical coordinates,” Journal of Computational Physics, vol. 105, no. 2, pp. 224–232, 1993.

P. Sigmund, “Theory of sputtering. I. Sputtering yield of amorphous and polycrystalline targets,” Physical Review Online Archive, vol. 184, pp. 383–416, 1969.

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. 1-2, pp. 65–86, 1983.

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.

J. Lindhard and M. Scharff, “Energy dissipation by ions in the kev region,” Physical Review, vol. 124, no. 1, pp. 128–130, 1961.

C. Shannon, Computational modeling of hall thruster plasma plume in a vacuum tank [M.S. thesis], Massachusetts Institute of Technology, Boston, Mass, USA, 2003.