Particle simulation of grid system for krypton ion thrusters

Molin CHEN,Anng SUN,Chong CHEN,Gungqing XIA

aScience and Technology on Combustion,Internal Flow and Thermo-Structure Laboratory,Northwestern Polytechnical University,Xi’an 710072,China

bState Key Laboratory of Electrical Insulation and Power Equipment,Xi’an Jiaotong University,Xi’an 710049,China

cState Key Laboratory of Structural Analysis for Industrial Equipment,Dalian University of Technology,Dalian 116024,China

1.Introduction

In the last decade,ion thrusters have been widely used for various space missions,1–3in which xenon(Xe)is frequently used as the propellant.However,Xe is typically expensive and a cost-effective propellant is highly demanded,from the economy point of review.Recently,krypton(Kr)attracts more attentions as an alternative propellant for different electric thrusters,for example,Hall thrusters4and ion thrusters.5

The performance of an ion thruster mainly depends on the acceleration process of the propellant in the grid system.Besides experimental diagnostics,numerical simulation provides a supplementary tool to investigate the transport process of the plasma in the grid system,and to estimate the performance and lifetime of the thruster.In the past,many simulations on the transport process of Xe plasmas in grid systems have been performed to explore the following aspects:the electric field distribution in a grid system,the motion trajectory of ion beam,the erosion on the grids,back-streaming,cut-off current and the ion extraction process.For instance,Wang et al.used a three-dimensional Particle In Cell(PIC)method to simulate the transport process of the Xe ion-beam in a single grid gate aperture.6Using particle methods,Peng et al.7and Sun et al.8examined the grid erosion and the effect of the deceleration grid for two-grid and three-grid systems,respectively.Wheelock et al.investigated the neutralization process in ion beam using PIC method.9Liu et al.investigated the characteristics of Charge EXchange(CEX)ions in ion thruster optical system with a two-dimensional axisymmetric numerical model.10,11Cao et al.conducted an in-depth study on the transport process of Xe ions in the grid system by combining PIC with the Immersed Finite Element(IFE)method.12Jia et al.analyzed the grid system performance of LIPS Xe ion thrusters.13However,most of the aforementioned simulations use Xe as the propellant(hereinafter referred to as Xe ion thrusters or Xe thrusters).In this paper,Kr is chosen as a novel propellant for ion thrusters(hereinafter referred to as Kr ion thrusters or Kr thrusters).The transport processes of Kr ions in the grid system at different acceleration voltages were simulated with PIC method.By analyzing the variations of the screen grid transparency,the accelerator grid current ratio and the divergence loss,the effects of acceleration voltage on the performance of Kr ion thruster were identified to guide the design of the grid system prototype of Kr ion thrusters.The results were also compared with Xe ion thrusters for the analyses of the advantage and disadvantage of different propellant choices.

The paper is organized as follows.In Section 2,the physical model and simulation method are introduced.Simulation results including the screen grid transparency,accelerator grid current ratio,and divergence loss are described and discussed in Section 3.Finally,we summarize our results and draw main conclusions.

2.Calculation model

2.1.Physical model and simulation domain

To avoid short-circuit induced by the grid system’s thermal deformation,the grid system in the classical NSTAR-30 ion thruster was used for the preliminary design and evaluations.

Table.1 lists the specific structural parameters of this kind of grid structure.6

Due to the symmetry of the grid structure,only a quarter of grid aperture was selected as the computational domain.14Asshown in Fig.1,(a)presents 3D structure of a complete grid aperture,(b)is the left view of grid aperture,(c)is the left view of the simulation domain,and(d)is the top view of the grid aperture.

Table 1 Parameters of grid system.

An equidistant grid was used on the computational region.Considering that the density of the plasma in the discharge chamber ranges from 1016m-3to 1017m-3,the spatial grid size was set as 5×10-5m and the grid number in calculation was set as 23×39×115,while the time step was set as 1.0×10-10s.

2.2.PIC/MCC model

PIC method is a kinetic method of simulating low temperature plasmas with the aim of tracking the motion of particles and the self-consistent electric field in a coupled way.15,16Monte Carlo Collision(MCC)method has been widely used for treating collisions between charged ions and neutral gas.17,18

In general,the PIC code consists of a cycle in every time step as follows15:(A)weighting the charge of the ions and electrons to the mesh nodes;(B)calculating the electric potential and the electric field of the calculation domain;(C)weighting the electrostatic field back to the ions;(D)moving the ions according to the second Newton law method with the electric field forces obtained above.The flowchart of PIC simulation can be seen in Ref.11.

In the model,the velocity and position of ions are calculated by the Newton-Lorentz law according to

Fig.1 Illustration of computational domain.

in whichmidenotes the mass of ion;edenotes the unit charge(since the ion considered in the model is monovalent,the carried charge is an element charge);v and x denote the velocity vector and position vector of ion respectively;E and B denote the electric field intensity and magnetic induction intensity respectively;tdenotes time.The magnetic field was neglected in calculations since it is very weak in grid system.Only the electric potential in the grid system was updated by solving the Poisson’s equation:

where ε0is the permittivity of vacuum,and φ,niandnedenote the electric potential,ion number density and electron number density respectively.To accelerate the convergence of the calculation,the successive over-relaxation method is used to solve the Poisson’s equation.

The number of ions injected into the simulation domain for each time step is decided by the ion number densitynat the inlet boundary.For the pre-sheath is set as the inlet boundary,the ion number densitynat the inlet boundary can calculated with equation:n=0.61n0,wheren0denotes the plasma number density in the discharge chamber and is varied with the beam current needs of different operation modes.

The initial axial velocity of ions injected into the simulation domain is set as the Bohm velocity9:wherekis the Boltzmann constant,andTeis the electron temperature and set as 5.0 eV.6

In the model,electrons are regarded as a fluid and their number density follows the Boltzmann distribution9:

wherene，ref,φrefandTe.refdenote the number density,potential and electron temperature of the plasma at the reference point respectively.In the model,the reference point was set as the discharge chamber in the upstream of accelerator grid,and as the neutral plane of plume downstream in the downstream of accelerator grid.9

The collisions between ions and background particles(Kr or Xe atoms)are described by MCC method.The background particles are assumed as a uniform distribution,with a temperature close to the discharge chamber wall(≈500 K),6and a number density equal to the neutral density upstream of the screen grid(≈1.5 × 1018m3).6Within a time step of Δt,the collision probability between the target particle and the background particle can be expressed as

in whichntdenotes the number density of the background gas,vincdenotes the velocity of the target particle,σT（εinc）denotes the collision cross section between particles,and εincdenotes the energy of the target particle.The collision cross section data of Xe are taken from Ref.9,and the collision cross section data of Kr are taken from Ref.19.

3.Simulation results and discussion

3.1.Screen grid transparency

The screen grid transparency ηsc=Ib/Ii,defined as the ratio of beam current,Ib,to the total ion current,Ii,from the discharge chamber that approaches the screen grid,1is an important parameter in assessing the efficiency of the grid system.Using the present model,the screen grid transparency corresponding to various beam currents was calculated under different acceleration voltagesVT.1Acceleration voltage means total voltage across accelerator gap,which is defined asVT=Vsc+|Vacc|and named total voltage too.Here,VscandVaccdenote the voltage of screen grid and accelerator grid respectively.In this paper,Vscwas set to different values,such as 600 V,800 V,1074 V,1400 V and 1800 V,butVaccwas set to a constant value of-180 V.So,different acceleration voltagesVTmean different screen grid voltagesVsc.

The variation of the screen grid transparency as a function of beam current is presented in Fig.2.For comparison,simulations with Xe as a propellant were also performed,while the distribution of ion number density and electric potential for both Kr ion thrusters and Xe ion thrusters are presented.Fig.3 presents the distribution of ion number density under the conditions ofJb=0.16 mA(Jb=nivaxialAapertureis the single-aperture beam current, whereAaperturedenotes the single-aperture area of screen grid,andvaxialdenotes the average axial velocity of ions)with different acceleration voltages,and Fig.4 presents the distribution of electric potential under the same operation conditions.

Fig.2 Screen grid transparency with beam current of Kr ion thrusters and Xe ion thrusters.

Fig.3 Distribution of ion number densities under different acceleration voltages.

In Fig.2,for both Kr ion thrusters and Xe ion thrusters,the curve of the screen grid transparency can divided into a smooth section and a fast descent section.The smooth section becomes larger at a higher acceleration voltage,and the screen grid transparency is higher too.The knee point of the curve represents the current threshold and corresponds to the extreme transparency when the grid structure and acceleration voltage are fixed.When the single-aperture beam currentJbis less than the current threshold,the screen grid transparency remains almost unchanged.WhenJbis greater than the threshold,the increasing plasma density compresses and weakens the sheath range and focusing effect,and some ions cannot be focused to enter into the screen grid aperture,leading to the decrease of the axial current intensity.

On the other hand,at the same acceleration voltage,the screen grid transparency ηscof Kr ion thrusters is slightly greater than Xe ion thrusters.This is due to the fact that Kr+is smaller in mass.The smaller mass means Kr+has a larger axial velocityvaxialthan Xe+with the same acceleration voltage.When the single-aperture beam currentJbis fixed,the ion number densityniof Kr ion thrusters is obviously smaller than Xe ion thrusters in the sheath on the upper of screen grid,which can be seen in Fig.3.The ion number density at the inlet is equal to that of Kr ion thrusters and Xe ion thrusters,but near the screen grid,Kr+number density is smaller than Xe+number density.In other words,Kr ion thrusters have a greater variation in ion number density from pre-sheath to screen grid.This different ion number density distribution leads to different electric potential distribution and results in different sheath structure.At the same acceleration voltage,the greater ion number density results in larger sheath area in Kr ion thrusters,as shown in Fig.4.

3.2.Analysis of accelerator grid current ratio

During the ion thruster operating process,some ions deviate from the main beam and impact into the surface of the accelerator grid and form an electric current.This current is then referred to as the accelerator grid current,which can reflect the lifetime of the grid system.The greater the accelerator grid current is,the more erosion on the grid occurs and the shorter lifetime the grid system has.

The accelerator grid currentJaand the beam current in the upstream of the accelerator gridJwere calculated.The ratio between these two factors is defined as the accelerator grid current ratio ηa,i.e.,ηa=Ja/J.Fig.5 presents the variation of ηawith differentJbat different acceleration voltages.

Fig.4 Distribution of electric potential under different acceleration voltages.

Fig.5 Variation of accelerator grid current ratio ηawith single-aperture beam current Jbof Kr ion thrusters and Xe ion thrusters.

In Fig.5,for both Kr ion thrusters and Xe ion thrusters,as the single-aperture beam current increases,the accelerator grid current ratio first changes slowly and then rises rapidly.The abscissa value of the curve’s infection point is the cutoff current threshold of the grid system,which has limited the operating beam current range of ion thrusters.The acceleration grid also affects the operating beam current range.By comparing the accelerator grid current ratios ηaat different acceleration voltages,one can observe that,at a high acceleration voltage,the curve of the accelerator grid current ratio has a long smooth segment,i.e.,the cutoff current threshold of the grid system increases and ion thrusters can operate within a larger beam current range.For example,for Kr ion thrusters,the maximum beam current of the single grid aperture increases from 0.28 mA to 0.65 mA,if the screen grid voltage increases from 1074 V to 1800 V.

On the other hand,the accelerator grid current ratio curve of Kr ion thrusters is significantly lower than that of Xe ion thrusters,and has a long smooth segment,at the same acceleration voltage.This is because the CEX collision cross section of Kr+smaller than Xe+,i.e.,fewer CEX ions were produced when Kr was used as the propellant,and low-energy CEX ions were the primary source of the accelerator grid current.The CEX collision rate,also called the CEX ion generating rate,was calculated in the model,and the distribution of CEX collision rate and current density on accelerator grid,under the conditions ofJb=0.16 mA andVsc=1074 V,are presented in Fig.6.In Fig.6,the CEX collision ratekCEXand current density on accelerator gridIerosionof Kr ion thruster are obviously smaller than Xe ion thruster.

3.3.Divergence loss

The momentum-weighted average plume divergence is defined in Eq.(6)as the ratio of the measured thrust component directed along the centerline of ion thrusters to the theoretical thrust achieved when all ions are traveling parallel to the centerline of ion thrusters.20in which θ denotes the divergence angle,˙mdenotes the ion mass flow rate,¯vdenotes the average velocity of ions in plume,Iaxialdenotes the axial current,andIbeamdenotes the total current of plume.

Fig.6 Distribution of CEX collision rate and current density on accelerator grid by CEX ion impact.

So,the momentum losses associated with plume divergence,called divergence loss,may be calculated with knowledge of the input mass flow,measured thrust,and the mass-weighted average velocity,which is defined in Eq.(7).

Fig.7 presents the variation of the thruster’s divergence loss withJbat different acceleration voltages.At different acceleration voltages,the beam’s focusing varies due to the variation of the ions’motion trajectories.

Compared with the curves of the Xe ion thrusters,the curves of the divergence loss of the Kr ion thrusters have a right shift.This indicates that,when the same grid structures and acceleration voltages are given,greater operating beam currents are required to achieve the grid system’s focusing for Kr ion thrusters.The mass difference between Kr+and Xe+is the reason of such shifting.Kr+respond easily to the radial electric field due to its smaller mass,aggregate towards the axis and make the focus point shift towards the left side of the screen grid,giving rise to the appearance of over-focusing.To ensure that the focusing position is reasonable,greater operating beam current and ion number density are required to enhance the potential along the axis of grid aperture,and then to reduce the deviation of the ions in the downstream of accelerator grid.

By combing the analyses of screen grid transparency,the accelerator grid current ratio and the divergence loss of Kr ion thrusters and Xe ion thrusters under different acceleration voltages,we can conclude that Kr ion thrusters have a larger operating beam current range and Kr is a suitable propellant for ion thrusters.

Fig.7 Variation of divergence losses with Jbat different acceleration voltages of Kr ion thrusters and Xe ion thrusters.

4.Conclusions

Using a 3D PIC method,the plasma transport processes in the grid system of a Kr ion thruster at different acceleration voltages were investigated,and the results,including screen grid transparency,accelerator grid current ratio and divergence loss,were compared with those of Xe ion thrusters.Main conclusions are addressed as follows:

(1)As the result of Xe ion thrusters,the screen grid transparency of Kr ion thrusters also decreases with the increase of acceleration voltage.But the screen grid transparency of Kr ion thrusters is slightly higher than that of Xe ion thrusters,which means a better propellant efficiency for ion thrusters.

(2)At the same acceleration voltage,the accelerator grid current ratio curve of Kr ion thrusters is far below the curve of Xe ion thrusters.This means that fewer ions will impact the accelerator grid and the grid system will have a long lifetime.On the other hand,the cutoff current threshold of Kr ion thrusters is larger than that of Xe ion thrusters,which provides the advantage of a large operating current or a high operating power for ion thrusters.

(3)Kr ion thrusters have a better divergence loss characteristic than Xe ion thrusters for the divergence loss curve shift,which means that Kr ion thruster can operate with a large beam current left at the same acceleration voltage.

(4)Kr ion thrusters have a larger operating beam current range and Kr is a suitable propellant for the design of high power ion thrusters.

Acknowledgements

This study was co-supported by the National Natural Science Foundation of China(No.11675040)and the Fundamental Research Funds for the Central Universities of China(Nos.3102014KYJD005 and 1191329723).

1.Goebel DM,Katz I.Fundamentals of electric propulsion:ion and hall thrusters.New Jersey:John Wiley&Sons Inc;2008.

2.Patterson MJ,Sovey JS.History of electric propulsion at NASA Glenn Research Center:1956 to present.Aerospace Eng2013;26(2):300–16.

3.Yamamoto N,Tomita K,Yamasaki N,Tsuru T,Ezaki T,Kotni Y,et al.Measurements of electron density and temperature in a miniature microwave discharge ion thruster using laser Thomson scattering technique.PlasmaSourcSciTechnol2010;19(4):45009–15.

4.Hause ML,Prince BD,Bemish RJ.Krypton charge exchange cross sections for Hall effect thruster models.Appl Phys2013;133(16):113–40.

5.Rawlin VK,Williams J,Pin˜ero LR,Roman RF.Status of ion engine development for high power,high specific impulse missions.Pasadena,USA;2001.Report No.:IEPC-01-096.

6.Wang J,Polk J,Brophy J.Three-dimensional particle simulations of ion-optics plasma flow and grid erosion.Reston:AIAA;2002.Report No.:AIAA-2002-2193.

7.Peng X,Ruyten WM,Friedly VJ,Keefer D.Particle simulation of ion optics and grid erosion for two-grid and three-grid systems.Rev Sci Instrum1994;65(6):1770–3.

8.Sun AB,Mao GW,Yang J,Xia GQ,Chen ML,Huo C.Particle simulation of three-grid ECR ion thruster optics and erosion prediction.Plasma Sci Technol2010;12(2):240–7.

9.Wheelock A,Cooke DL,Gatsonis N.Investigation of ion beam neutralization processes with 2D and 3D PIC simulations.Comput Phys Commun2004;164(1):336–43.

10.Liu C,Tang HB,Gu Z,Jiang HC.Particle simulation of ion thruster optical systems using single linked lists.High Power Laser Particle Beams2006;18(7):1193–8.

11.Zhong LW,Liu Y,Li J,Gu Z,Jiang HC,Wang HX,et al.Numerical simulation of characteristics of CEX ions in ion thrust eroptical system.ChinJAeronaut2010;13(1):15–21.

12.Wang EM,Chu YC,Cao Y,Juan L.Numerical simulation of erosion mechanism for ion thruster accelerator grid aperture walls based on IFE-PIC and MCC methods.High Voltage Eng2013;39(7):1763–71[Chinese].

13.Jia YH,Li ZM,Zhang TP,Li J.Effect analysis of plasma density on ion beam extracted by grid system.Chin Space Sci Technol2012;32(3):72–7.

14.Chen ML,Xia GQ,Mao GW.Three-dimensional particle in cell simulation of multi-mode ion thruster optics system.Acta Phys Sin2014;63:182901[Chinese].

15.Birdsall CK,Landon AB.Plasma physics via computer simulation.New York:McGraw-Hill;1985.p.11–3.

16.Miyake Y,Usui H.Particle-in-cell modeling of spacecraft-plasma interaction effects on double-probe electric field measurements.Radio Sci2017;51(12):1905–22.

17.Birdsall CK.Particle-in-cell charged-particle simulations,plus Monte Carlo collisions with neutral atoms,PIC-MCC.Plasma Sci IEEE Trans Plasma Sci1991;19(2):65–85.

18.Sarikaya CK,Rafatov I,Kudryavtsev AA.PIC/MCC analysis of a photoresonance plasma sustained in a sodium vapor.Phys Plasmas2017;24(8):535–54.

19.Johnson RE.Charge transger and f i ne structure transitions in Kr+(2Pj)+Kr and Xr+(2Pj)+Xe collision.Phys Soc Jpn1972;32(6):1612–4.

20.Brown DL,Larson CW,Beal BE,Gallimore AD.Methodology and historical perspective of a hall thruster efficiency analysis.J Propul Power2009;25(6):1163–77.