Performance Simulation of a 5 kW hall Thruster

Hall thruster is a kind of plasma optics device, which is used mainly in space propulsion. To simulate the discharge process of plasma and the performance of a 5 kW hall thruster, a two-dimensional PIC-MCC model in the R-Z plane is built. In the model, the anomalous diffusion of the electrons including Bohm diffusion and near-wall conduction is modeled. The Bohm diffusion is modeled by using a Brownian motion instead of the Bohm collision method and the near-wall conduction is modeled by a secondary electron emission model. In addition to the elastic, excitation, and ionization collisions between electrons and neutral atoms, the Coulomb collisions are included. The plasma discharge process including the transient oscillation and steady state oscillation is well reproduced. First, the influence of the discharge voltage and magnetic field on the steady state oscillation is simulated. The oscillation amplitude increases as the discharge voltage gets larger at first, and then decreases. While the oscillation amplitude decreases as the magnetic field gets stronger at first, and then increases. Later, the influence of the discharge voltage and mass flow rate on the performance of the thruster is simulated. When the mass flow rate is constant, the total efficiency initially increases with the discharge voltage, reaches the maximum at 600 V, and then declined. When the discharge voltage is constant, the total efficiency increases as the mass flow rate rises from 10 to 15 mg/s. Finally, a comparison between simulated and experimental performance reveals that the largest deviation is within 15%, thereby indirectly validating the accuracy of the model.


INTRODUCTION
Hall electric propulsion is one of the most mature and widely used electric propulsion technologies (Pidgeon, 2006;Mathers et al., 2009). LHT-140 ( Figure 1) is a 5 kW hall thruster developed by the Lanzhou Institute of Physics. At present, the discharge characteristics of the plasma and the performance of a hall thruster are evaluated mainly by conducting high-cost and time-consuming experiments. As research on plasma discharge mechanism becomes in depth and numerical methods improve continuously, simulation plays an increasely important role in understanding the micro physical process of hall thruster (Arkhipov and Bishaev, 2007;Hofer, 2008;Zhang et al., 2011;Zhang et al., 2014). Lentz built a one-dimensional particle-in-cell (PIC) model (Lentz, 1993) and simulated the plasma characteristics and performance of a Japanese hall thruster. Noguchi, Martinez-Sanchez and Ahedo studied axial plasma discharge oscillation by using a one-dimensional PIC model (Noguchi et al., 1999). Hirakawa built a two-dimensional PIC model (Hirakawa and Arakawa, 1996) and simulated the electron diffusion. Later, Beidler (1999), Szabo (2001), Sullivan (2004) and Adam et al. (2004) improved the simulation method continuously.
In most PIC models, Bohm collision was applied to simulate the anomalous diffusion of electrons (Koo and Boyd, 2006;Cho et al., 2013). The limitation is that not all the probable diffusion magnitudes can be simulated. In this paper, we applied the method of Vincent (2002) to simulate Bohm diffusion, to allow the simulation of any diffusion magnitude. In addition, in most present models, the Coulomb collisions between electrons and ions, as well as electrons and electrons, are excluded, we simulated the Coulomb collisions by using the method of Szabo (2001).
This paper aims to evaluate the performance and discharge characteristics of LHT-140 by using the built PIC model. First, the numerical methods applied in the building process of the model are introduced. Second, the plasma discharge characteristics and thruster performance are simulated. Finally, the simulated and experimental performance results are compared.

Simulation Region
The simulation region ( Figure 2) consists of the discharge channel and the near plume region. The length and width of the near plume region is 1.5 and 4.5 times that of the discharge channel, respectively. The boundary of the discharge channel consists of a metal anode wall, ceramics wall, symmetric axis, and free space.

Model Simplification
Lightening ions using an artificial mass ratio and increasing the vacuum permittivity using artificial permittivity are conventional methods to simplify the particle-in-cell (PIC) model (Szabo, 2001). Speeding up the ions could shorten the time of convergence. In this study, we apply an artificial mass ratio of f M real /M comp 1,000, which indicates that the ions are 1,000 times lighter than they should be. Increasing the vacuum permittivity enables us to use a coarser grid. In this study, we increase the vacuum permittivity by 2,500 times. Consequently, the size of rectangular grid applied in our model could be magnified by 50 times (Szabo, 2001). The LHT-140 Hall thruster is similar to the 5 kW P5 Hall thruster (Vincent, 2002), and the Debye length of 0.02 mm in reference (Vincent, 2002) is used in our model. Finally, the mesh size becomes 1 mm after increasing the vacuum permittivity.

Potential and Electric Field Solving
The potential was acquired by solving the Poisson equation. The boundary conditions of the potential are as follows: • At the anode surface, the potential is set to discharge voltage.
• At the symmetric axis, the normal electric field E ⊥ 0.
• At the free space boundary, the potential is set to 0 V. • At the boron nitride ceramic wall surface, the boundary condition is (Vincent, 2002): For LHT-140 hall thruster, both in and outside the discharge channel, ε dieletric E dielectric ⊥ is negligible compared with E plasma ⊥ , similar to the situation for the P5 thruster (Vincent, 2002). Thus, ε dieletric E dielectric ⊥ could be neglected, and Eq. 2 becomes: The electric field was solved by Eqs 4, 5 (Birdsall and Langdon, 1991): where E x j,k and E y j,k are the electric field on the grid point (j, k) in the x and y directions, respectively, and ϕ k+1,j−1 , ϕ k+1,j+1 , ϕ k,j−1 , ϕ k,j+1 , ϕ k−1,j−1 , ϕ k−1,j+1 , ϕ k−1,j and ϕ k+1,j are the electric potential on the grid points around (k, j), and Δx is the width of the square grid.

Moving the Particles
Both neutrals and charged particles were modeled by PIC-MCC method. After the electric filed was determined, electrons and ions moved by the electric and magnetic fields were simulated using the leapfrog method (Birdsall and Langdon, 1991). The neutrals were not affected by the electric and magnetic fields, and there's no need to update their velocity during each time step if there's no collision. Because the density of neutrals is much larger than that of charged particles, a larger weight should be used for neutrals, which was 60 in our model.

Cathode Modeling
The cathode is modeled by injecting electrons from the upper free boundary with 0.1 eV of kinetic energy. The number of electrons injected in each time step is determined by current balance.
Cathode current getting into the discharge region I cd , ion current leaving from the free boundaryI + b , and electron current leaving from the free boundary I az satisfy the following constraints: From Eq. 7, the number of electrons injected in each time step could be determined.

Anomalous Diffusion
The Bohm diffusion and near-wall conduction are the main factors that lead to the anomalous diffusion of electrons. Usually, a virtual collision called "Bohm collision" is introduced to simulate the Bohm diffusion. The Bohm collision frequency ] β , is acquired by solving the following equation: 1 16 kT e eB kT e m e 1 ω c β 1 + β 2 (7) β ω c /(] c + ] β ) is the hall parameter that indicates the magnetization of electrons. Two limitations occur. Within one time step, the direction of electron velocity can be randomized only once, thus, ] β < 1/Δt. When β 1, ] β reaches the maximum in Eq. 8.
Thus, not all the probable diffusion magnitudes can be simulated. Vincent (2002) applied a Brownian motion to simulate the Bohm diffusion. After the electron was moved by the leapfrog algorithm, it was then moved in the two-dimensional plane perpendicular to the local magnetic field where D is the electron diffusion coefficient; Δt is the time step; and R 1 and R 2 are uniform random numbers between 0 and 1. Through this approach, any Bohm diffusion magnitude can be simulated.

Coulomb Collision
According to J. J. Szabo's method (Szabo, 2001), the Coulomb collision frequency can be calculated as follows: ] ee n e Q ee |v e |.
The collision section can be calculated as follows: where ln Λ is the Coulomb logarithm.

Near-Wall Conduction
A secondary electron emission model of Richard R. Hofer (Hofer et al., 2007) was applied to simulate the near-wall conduction of electrons.
where Γ is the gamma function (a 0.123, b 0.528). The number of secondary electrons emitted from the ceramics wall after an electron collision is where δ w means rounding δ w down to an integer, and R f is a random number between 0 and 1. If R f < δ w − δ w , the primary electron is reflected back into the plasma, otherwise accumulates on the ceramics wall. The velocities of the secondary electrons after colliding against the wall are determined as follows: where v 1 and v 2 are perpendicular to the direction of effusion and v 3 is in the direction of effusion.
R 1 , R 2 , and R 3 are random numbers between 0 and 1.

Performance Model
The main performance of a hall thruster includes thrust, specific impulse and efficiency. The thrust is calculated by summing the product of mass flow rate and axial speed of ions that leave from the thruster exit.
αrepresents the charge state of ions, α 1 represents singlecharged ions, and α 2 represents double-charged ions. The anode efficiency of the thruster is where P d V d I c is the discharge power, and _ m a is the anode mass flow rate. The acceleration efficiency of ions is defined as the ratio of the ion current to the discharge current The utilization efficiency of propellant is defined as the ratio of the mass flow rate of ions leaving the discharge channel to that of propellant atoms entering the discharge channel The efficiency of electrical energy converted into kinetic energy of ions is defined as: The relationship between the total efficiency and other efficiencies is where c is the thrust loss factor.

Plasma Discharge
The plasma discharge process of LHT-140 when the discharge voltage is 300 V and the propellant mass flow rate is 10 mg/s is simulated. Figure 3 shows the change of currents with time. The process consists of transient oscillation and steady-state periodic oscillation. During the transient oscillation, the discharge and ion   currents increase quickly, reaching a maximum, which is far larger than the steady-state value at approximately 15 μs. Figure 3 shows that the transient maximum value of the discharge current is twice that in the steady state. After the transient reached the maximum value, the currents decreased and then oscillated periodically, which indicated the beginning of the so-called "breathing oscillation" (Fife et al., 1997;Boeuf and Garrigues, 1998). Figure 4 presents the change of maximum value (I max ), minimum value (I min ), and average value (I ave ) of the discharge current and the oscillation frequency (f osc ) with discharge voltage. The discharge voltage significantly affected the oscillation. When the discharge voltage increased from 300 to 700 V, the oscillation amplitude (I max −I min ) increased from 1.8 to 7.1 A. When the discharge voltage increased from 700 to 900 V, the oscillation amplitude decreases from 7.1 to 3.1 A. The oscillation frequency increased from 12 to 20 kHz. Figure 5 displays the change of I max , I min , I ave , and f osc with magnetic field intensity. We can see that the magnetic field also exhibited a great influence on the oscillation. The oscillation amplitude decreased from 4.6 to 1.8 A when the magnetic field increased from 150 to 300 Gs and increased from 1.8 to 4.3 A when the magnetic field increased from 300 to 450 Gs. The oscillation frequency decreased from 17 to 10 kHz.
Figures 6A-C show the time-averaged two-dimensional distributions of the plasma potential, density, and ionization rate during 1,000 steps (almost 1/5 of the oscillation period), respectively. As shown in Figure 6A, the ionization region is located at one to two normalized distances from the anode. Figure 6B indicates that the largest plasma density is 2 × 10 12 cm −3 and 0.6 normalized distance upstream of the thruster exit. Figure 6C illustrates that the potential drop occurs mainly in the ionization region, where the magnetic field is strong and the electron conductivity is low. The strength of the electric field must be high to compensate for the low electron conductivity, and to guarantee the continuity of the discharge current.

Thruster Performance
Performance experiments (Figure 7) with different discharge voltages and propellant mass flow rates of LHT-140 were conducted in a large vacuum chamber TS-7 with a diameter and a length of 3.8 and 8.5 m, respectively. The pumping speed of TS-7 can reach 1.8 × 10 5 L per second, providing a background pressure of around 2.0 × 10 −3 Pa during testing. The discharge voltage ranged from 300 to 900 V, and the propellant mass flow rate which was regulated by a flowmeter ranged from 10 to 15 mg/ s, while the magnetic field remained unchanged. Thrust was measured by a laser interferometer (showed in Figure 7) through the relationship between thrust and displacement. Thrust measurements were conducted in intervals of about 40 min to minimize the thermal drift. Considering the thermal drift and repeatability, the average variation in thrust for a given point was estimated to be about ±2%. Given the uncertainty of the thrust, mass flow rate, current and voltage, the uncertainty for efficiency was about ±3.1%.
The simulated results of the total, acceleration, utilization, and energy conversion efficiencies are shown in Figure 8 The total efficiency initially increased with the discharge voltage, reached the maximum at 600 V, and then declined. Blateau's simulation study (Vincent et al., 2001) showed that the magnetic field and discharge voltage require a reasonable match. After the discharge voltage increased to a certain value, the magnetic field must increase as well. However, Figure 8 shows that the magnetic field remained unchanged. After the discharge voltage increased to 600 V, the magnetic field was not large enough to constrain the electrons, which indicated that the electrons were absorbed by the anode before the ionization collisions with the propellant atoms. Thus, the total efficiency declined after the discharge voltage increased to 600 V. The acceleration efficiency is insensitive to the change in discharge voltage, and remains at 77%. A similar relationship is found in experiments (Ashkenazy et al., 1998) for a hall thruster with a diameter of 70 mm. Moreover, with the increasing discharge voltage, the ionization of propellant atoms becomes increasingly adequate, thereby increasing the utilization efficiency. Once the magnetic field is no longer large enough to constrain the electrons, the utilization efficiency starts to decline. Figure 9 shows a comparison of the experimental and simulated results of total efficiency when the discharge voltage ranged from 300 to 900 V and the propellant mass flow rate was 10 mg/s. The tested maximum value is 58%, and   the simulated maximum value is 51%. When the discharge voltage increased from 300 to 600 V, ηincreased with η u . When the discharge voltage increased from 600 to 900 V, η decreased with η u . The simulated results are smaller than the experimental results by 9-13% because the influence of background pressure in the vacuum chamber is not included. Figure 10 shows a comparison of the experimental and simulated results of thrust when the discharge voltage ranged from 300 to 900 V and the propellant mass flow rate was 10 mg/ s. As in the previous analysis, when the discharge voltage increased from 300 to 600 V, the ionization became increasingly adequate, the amount of produced ions increased, and the speed of the ejected ions became increasingly fast. Thus, the thrust increased with the discharge voltage. When the discharge voltage increased from 600 to 900 V, the ionization efficiency decreased dramatically and the number of produced ions decreased, but the speed of the ejected ions still increased. Thus, the thrust still increased with a slower speed. The simulated results are smaller than the experimental results by 8-14%. As mentioned earlier, the influence of background pressure in the vacuum chamber is not included. Figure 11 shows the relationship between thrust and efficiency. The maximum simulated efficiency is 51.1%, and the corresponding thrust is 315 mN. The maximum experimental efficiency is 58%, and the corresponding thrust is 350 mN. The trend demonstrated that efficiency increased as    thrust increased initially and then decreased as thrust increased continuously.
Only single-charged ions Xe + and double-charged ions Xe ++ are included. Figure 12 shows the ratio of Xe + current and Xe ++ current to the total ion current when the discharge voltage ranged from 300 to 900 V and the mass flow rate is 10 mg/s. The ratio of Xe + current decreased from 96 to 82%, while the ratio of Xe ++ current increased from 4 to 18%. In accordance with the variance of propellant utilization efficiency or ionization efficiency, the ion current increased when the discharge voltage increased from 300 to 600 V and then decreased when the discharge voltage increased from 600 to 900 V. The trend of beam current, Xe + and Xe ++ fractions were similar to that in reference (Vincent et al., 2001;Vincent, 2002), in which, the Xe ++ fraction increased from 20.2 to 30.3% when the discharge voltage ranged from 300 to 1200 V, and the beam current increased at first, then started to decrease.
The frequency of ionization collisions increased with the increase in the propellant mass flow rate. Consequently, the utilization or ionization efficiency increased as well, as shown in Figure 13.

CONCLUSION
The discharge process of the plasma, including transient and steady-state oscillations, is well reproduced. During transient oscillation, the discharge current increases rapidly to a value that is far larger than the steady-state value. These characteristics must be considered when designing the power processing unit. During steady-state oscillation, the discharge amplitude increased first and then decreased as the discharge voltage increased, while it decreased first and then increased as the magnetic field became stronger. The performance of LHT-140 was well predicted by the built model with a deviation of less than 15% compared with the experimental results, while the deviation was mainly caused by the background pressure of the vacuum chamber. Both discharge voltage and propellant mass flow rate have a great influence on the performance of the thruster. The trend of the total efficiency was consistent with either the propellant utilization efficiency or the ionization efficiency. As the discharge voltage increased, the thrust became larger and the specific impulse became higher, but the efficiency of the thruster decreased. The research results of this paper is helpful to determine a reasonable range of discharge voltage. This paper is helpful for the optimization of the performance and selection of ideal operating points for the LHT-140 hall thruster.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.