Axisymmetric Hybrid Plasma Model for Hall Effect Thrusters

Abstract

Hall Effect Thrusters (HETs) are nowadays widely used for satellite applications because of their efficiency and robustness compared to other electric propulsion devices. Computational modelling of plasma in HETs is interesting for several reasons: it can be used to predict thrusters’ operative life; moreover, it provides a better understanding of the physical behaviour of this device and can be used to optimize the next generation of thrusters. In this work, the discharge within the accelerating channel and near-plume of HETs has been modelled by means of an axisymmetric hybrid approach: a set of fluid equations for electrons has been solved to get electron temperatures, plasma potential and the discharge current, whereas a Particle-In-Cell (PIC) sub-model has been developed to capture the behaviour of neutrals and ions. A two-region electron mobility model has been incorporated. It includes electron–neutral/ion collisions and uses empirical constants, that vary as a continuous function of axial coordinates, to take into account electron–wall collisions and Bohm diffusion/SEE effects. An SPT-100 thruster has been selected for the verification of the model because of the availability of reliable numerical and experimental data. The results of the presented simulations show that the code is able to describe plasma discharge reproducing, with consistency, the physics within the accelerating channel of HETs. A small discrepancy in the experimental magnitude of ions’ expansion, due probably to boundary condition effects, has been found.

1. Introduction

Hall Thrusters are spacecraft propulsion devices that utilize applied magnetic fields and closed electron Hall drift to accelerate quasi-neutral plasma [1]. The typical Hall Thruster consists of an annular chamber with one end open to the ambient environment (Figure 1). A neutral propellant is injected through the closed end, which also contains the anode. An externally located cathode produces electrons, a fraction of which enters the chamber and ionizes the propellant, and the other part neutralizes the ions in the plume. In order to increase the electron transit time, and hence improve the ionization efficiency, magnetic field (B) is applied over a section of the chamber. The magnetic field strength is selected such that electrons become magnetized and trapped in a closed azimuthal Hall drift at about the thruster centreline. The magnetic field thus plays an important secondary role. Since the field restricts motion of electrons in the normal direction, electrons redistribute along the magnetic field line according to the spatial variation in ion density. The magnetic field lines thus become lines of constant potential and an electric field (E), which develops in the direction normal to the magnetic field, accelerates the ionized propellant. The direction of acceleration thus depends on the topology of the magnetic field. In most of the HETs configurations, the magnetic lines are not perfectly perpendicular to the chamber. Thus, a beam divergence occurs that pushes ions also toward the channel. Ions with enough energy collide with the wall and remove part of material due to sputtering. The channel wall material (typically a ceramic) protects the magnetic circuit from hazardous particles. Once the magnetic circuitry is exposed to plasma, the thruster’s end of life is declared.

Despite over 40 years of flight heritage, the community still requires tools capable of modelling these devices to predict their lifetime and because, in spite of important research efforts, there is still the lack of understanding of such phenomenology that limits their full potential development, e.g., electrons transport across the magnetic field [2]. For the last purpose, various numerical models have been developed. Time-accurate solvers are used because HETs are subjected to various types of instabilities and oscillations, e.g., the ionization instability (called “breathing mode”), which is due to the periodic depletion and replenishment of neutral atoms in the ionization region and leads to current oscillations in the 10–20 kHz frequency range [2].

Simulations address either the near or far regions of the plasma plume, since the regimes in which they operate are quite different from the point of view of the “main” physical effects in the plasma. Typical approaches for the channel and near-plume regions are: kinetic, fluid and hybrid modelling.

In kinetic models (cf., e.g., [3,4,5,6,7,8]), electrons and heavy species (ions and neutrals) are treated as discrete particles and their trajectory is followed in the space phase considering the action of electric and magnetic fields. Monte Carlo Collisions (MCC) algorithms are often used to manage collision between particles. The majority of particle codes are based on the Particle-In-Cell (PIC) methodology, in which macro-particles, representatives of a number of real particles, are first sorted and “weighted” to a spatial mesh, which is utilized to calculate the electric and magnetic fields. This method provides the velocity distribution function (VDF) of species, but computational cost is very high because of the small time step required to follow electrons.

In fluid modelling (cf., e.g., [9,10,11]), all species involved in the plasma discharge are modelled as fluids. The fluid equations use the framework of the Boltzmann equation to describe the evolution of macroscopic quantities of interest. Each equation depends on moments, whose evolution is given by equations of higher order; therefore, the equations need to be truncated, or “closed”, at a certain order. In order to do this, an assumption can be made on the shape of the VDF (e.g., Maxwell distributions), which is used to derive a closed expression for the moment of the highest order to be conserved in the equations. This approach is computationally efficient, but the main limitation is due to the inability to solve for “kinetic” effects, which are those associated with how the various processes in the discharge affect the shape of VDF.

A compromise between previous models is represented by hybrid modelling (cf., e.g., [12,13,14,15,16,17,18,19]), typically choosing a kinetic approach for the heavy species populations, and the fluid approach for the electron population. Among the kinetic approaches, PIC is preferred to other combinations (such as Direct Kinetic, solving kinetic equations [20]) because it allows for describing accurately the VDF of heavy species, reducing computational efforts. The main limitation of hybrid-PIC simulations is the necessity of anomalous coefficient’s fluid electron description. It has been known since the early studies of HETs that collisions between electrons and neutral atoms (defined as the “classic” contribution mechanism to electrons’ mobility) are not sufficient to explain electron transport across the magnetic field lines, especially in the exhaust region where the magnetic field is maximum and the atom density is very small due to intense ionization inside the channel [21]. Various theories [1,16,17,22,23] have been proposed to explain this anomalous electron transport (e.g., electron–wall scattering, field fluctuations, Bohm diffusion) but a self-consistent model, able to accurately predict the operation of HETs and to help in designing new thrusters, is yet to be determined. Despite this limitation, hybrid models have been one of the preferred approaches over the past two decades in the simulation of HETs’ acceleration channel and plume, being able to model some of the high-level physical interactions that are responsible for the behaviour of these devices.

The HYbrid PIC-FLUid model (HYPICFLU) described in this paper is based on the legacy and lessons learned mainly from the HYbrid PIC-Hall (HPHall) code by Fife [13], and its development, Hybrid PIC Hall-2 (HPHall2) [16,17,24], but it takes also ideas from Hageelar et al. [18] and Kawashima et al. [15].

HPHall is one the first hybrid codes to model HETs that is able to reproduce “breathing mode” oscillations. It uses a PIC approach on an axisymmetric domain (2D radial-axial) and quasi-one-dimensional fluid equations for electrons, with magnetic field lines as spatial coordinates (thermalized electrons assumption, cf. [25]). The domain includes the channel and part of the near-plume region that is shaped with a circular sector. In order to accommodate the domain, a non-uniform spatial grid is used for the PIC method. A mixed electron mobility model (“classic”+“anomalous”) with one empirical parameter accounting for Bohm diffusion is used. In order to consider the effect of the Secondary Electron Emission (SEE) on the mobility, a total near the wall-current is modelled and included in the current conservation equation. Ions are generated as the result of electron–neutral collisions. Despite that only single-charge ions are tracked by PIC, an ionization rate for double-charge ions, a function of the single-charge ions density, can be enabled. Finally, HPHall uses an energy independent cross section for electron–neutral (Xenon) scattering of 30 × 10${}^{-20}$ m${}^2$.

HPHall-2, the evolution of HPHall, includes the following main changes: the modification of the algorithms for injection, reflection, and recombination of neutrals at the walls; the inclusion of ion-neutral collisions and plasma–wall interaction (i.e., the fulfilment by the ion flow of the Bohm condition required by the sheath solution, the inclusion of the “wall-collisionality” phenomenon, the adoption of new sheath model [24] which distinguishes between primary and secondary electrons and takes into account the space-charge saturation for high SEE). HPHall-2 uses a two region mixed mobility model for the electrons: the value of the empirical parameter for Bohm diffusion differs from inside to outside of the channel and a linear variation allows for a smooth transition at the exit of channel. Furthermore, electron–ions and electron–wall collision frequencies are included in the contribution of the mobility. The first one is modelled using the classic formula with the Coulomb logarithm [26] whereas the second one is derived by the new sheath model implemented, which is a function of axial coordinates and several plasma properties [17]. Finally, an energy dependent function of electron–neutral frequency, based on the Maxwellian average of the measured Xenon cross section of Nickel et al. [27], is considered.

HYPICFLU inherits from HPHall/HPHall-2 the hybrid approach using PIC on an axisymmetric domain (2D radial-axial) and 1D approximation to the electron fluid, but, to simplify the implementation process and several functions (e.g., the “sorting” and “weighting” function), a uniform Cartesian mesh on rectangular domain shape, including the channel and near-plume region, has been selected. The electron flux on the upper and lower boundaries outside the channel is set to be zero according to Kawashima et al. [28]. With the purpose to study in a first extent low-discharge voltage thruster (less than 300 V), the effects of multiple charged ion species on the plasma response have been neglected, as common practice in the literature [16,17,18,19]. The mobility model is similar to that used by HPHall-2 in which it differs by: the calculation of electron–neutral collision frequency which uses a new correlation for the most recent elastic electron–neutral cross section for Xenon [29]; the calculation of the electron–wall collision frequency by means of the simpler model proposed by Hageelar et al. [18], which considers wall collision to be constant within the channel throughout an empirical parameter (in $0.1÷1$ range) that multiplies the typical measured frequency of ${10}^7$ s${}^{-1}$; a continuous function that allows for the transition of Bohm empirical parameters from the value set within the channel to that outside. The location and the width of the axial position interval in which the variation occurs are considered as free parameters. Moreover, because near-wall current is not calculated, the effect of SEE on the mobility is incorporated into the internal Bohm coefficient.

Finally, the sheath model is similar to HPHall [12] and the algorithms for injections, ionization and wall recombination are inherited from HPHall-2 [17].

The aim of this paper was to describe and document the principles of the model implemented in HYPICFLU and the simulation results of SPT-100 discharge are presented in order to verify the simulation potential of the models implemented with respect to relevant experimental data and numerical models from the literature.

2. Hybrid PIC-Fluid Approach for HET Plasma

The unknown plasma variables of the HYPICFLU model are:

• $n_i$, ion number density;
• $n_n$, neutral particle density;
• $\dot{n_i}$, ionization rate;
• ${\overrightarrow{u}}_i$, ion velocity vector;
• ${\overrightarrow{u}}_n$, neutral velocity vector;
• $T_e$, electron temperature;
• $\varphi$, plasma potential.

The model is composed of two parts: the Particle-In-Cell (PIC) sub-model, for the heavy species (i.e., ions and neutral), and the fluid one, for electrons. Figure 2 shows the data exchanged between them. All variables are time dependent, except for magnetic field, $\overrightarrow{B}$. This means that the sub-models interact iteratively since the outputs of one are used as inputs of the other and vice-versa.

Particularly, fluid electrons equations require as input the number density of both ions and neutrals ($n_n,n_i$), the velocity of ions (${\overrightarrow{u}}_i$) and the ionization rate (${\dot{n}}_i$) and return the electron temperature ($T_e$) and the plasma potential ($\varphi$), from which we derive the electric field ($\overrightarrow{E}$) that accelerates ions in the channel; on the contrary, the heavy species sub-model needs as inputs electron temperatures and electric fields and gives back both the number density of ions and neutrals, the velocity of ions and the ionization rate. Neutrals velocity (${\overrightarrow{u}}_n$) is computed by PIC but it is not required by the fluid sub-model.

Bulk plasma is quasi neutral (i.e., $n_i\simeq n_e$) and a proper sheath model for plasma interaction with walls are taken into account in the fluid equations.

The magnetic field is an external input of only the fluid sub-model because only electrons are magnetized.

2.1. Heavy Particles

Heavy particles, i.e., ions and neutrals, are modelled with the PIC method. Each simulated particle represents an agglomerate of real particles, typically ${10}^8$ to ${10}^{12}$, saving computational time (a code simulating each particle would be resource intensive, indeed). The number of particles represented by one computational macro-particle is called specific weight and, in general, it varies among the particles and changes with time, e.g., due to ionization. The mass of one macro-particle is obtained by multiplying the mass of the single particle with the specific weight.

One limit of the PIC method is the introduction of numerical oscillation (“noise”) due to its inability to simulate a continuum of particles. The specific weight plays a major role in this sense: the lower it is, the less the numerical fluctuations are because more particles are simulated.

The governing equations for the PIC method are the laws of dynamics. Particularly, the second one, written for a generic charge particle subjected to both electric ($\overrightarrow{E}$) and magnetic fields ($\overrightarrow{B}$), is:

$$\overrightarrow{a}=\frac{q}{m}\left[\overrightarrow{E}+\overrightarrow{u}\times \overrightarrow{B}\right],$$

(1)

where $\overrightarrow{a}$ is the acceleration, $q/m$ is the charge to mass ratio of the macro-particle and $\overrightarrow{u}$ is its velocity.

The second term on the RHS of Equation (1) can be neglected since both neutrals and ions are not magnetized (i.e., their motion is slightly influenced by the magnetic field because their Larmour radius is larger than the typical channel length). Hence, the acceleration is:

$$\overrightarrow{a}=\frac{q}{m}\overrightarrow{E}$$

(2)

For neutrals, having no charge, it reduces to $\overrightarrow{a}=0$. Equation (2) can be integrated each time-step in order to update particles’ velocity and position.

Neutrals are injected from the inlet that is modelled as an annular ring positioned at the centre of anode wall. Velocities of new-injected particles are sampled from a semi-Maxwellian distribution characterized by anode wall temperature [30].

Neutral–neutral and ion–neutral collision are neglected within the channel because typically neutrals mean the free path is higher than the channel length. Hence, a macro-neutral changes its velocity only if a wall collision occurs. Following Ref. [16], neutrals that undergo a collision with the walls are scattered back in the domain diffusively with a velocity sampled from a Maxwellian distribution at the wall temperature. Usually, the gas–surface interaction is modelled with the Maxwell scattering model, in which a fraction $a_w$ of the colliding particles are reflected diffusively while the remaining $(1-a_w)$ are reflected specularly, where $a_w$ is the accommodation factor. However, the complete accommodation, that is the purely diffusive reflection, describes well the interaction of the neutrals with the walls of a Hall Thruster [31]. In fact, the sputtering and the reposition tends to increase the roughness of the walls inside HETs, which is the primary cause of the random distribution of the reflected particles [16].

It is assumed that every ion colliding with the wall recombines into a neutral, with a velocity that, again, is sampled from a semi-Maxwellian distribution at the wall temperature.

Only the ionization due to neutral–electron collisions is modelled and only single ionized ions are considered. Double ionized ions can produce appreciable effect on the performances, especially for thrusters operating at high power levels; however, a precise characterization of the discharge is beyond the aims of this work. The ionization rate, i.e., the rate at which ions are created per unit of volume, can be expressed as [13]:

$$\dot{n_i}=n_e{n}_n\zeta \left(T_e\right)$$

(3)

Here, $\dot{n_i}$ is the ionization rate, $\zeta \left(T_e\right)$ is the ionization coefficient while $n_e$ and $n_n$ are respectively the electron and neutral number densities (for quasi-neutrality $n_e\simeq n_i$). $\zeta \left(T_e\right)$ is a function only of the electron temperature (e.g., approximate formula for Xenon can be found in [1]). The bulk recombination has not been modelled because it is several orders of magnitude inferior to the bulk ionization rate, as shown in [13].

Finally, it is worth underlining that, by neglecting the neutral–ion collisions, this excludes the possibility of charge exchange (CEX), which produces slow ions and fast neutrals, everywhere. The effect of CEX is significant in the plume where it can lead to the attraction of slow ions toward the thruster. However, since CEX modelling requires the development of complex algorithms such as Monte Carlo Collision (MCC) that are computational demanding and, since the effect in the thruster performances is quite limited, the CEX has been neglected.

2.2. Electrons

The electron fluid equations used in the fluid sub-model are: the electrons momentum equation, current conservation and electron energy equations. The parameters that can be computed are the electron temperature ($T_e$), the discharge current ($I_a$), the plasma potential ($\varphi$) and the electron velocity normal to magnetic streamlines ($u_{e,\widehat{n}}$). As a result of quasi-neutrality in the bulk plasma, electron density ($n_e$) is known as the result of the PIC sub-model, i.e., $n_e\simeq n_i$.

2.2.1. Momentum Equation

Since the electrons inside the acceleration zone are magnetized and assuming that they reach thermal equilibrium, the electron momentum equation, along a magnetic field line, can be written as:

$$K_B{T}_e\frac{\partial n_e}{\partial \widehat{t}}=e{n}_e\frac{\partial \varphi }{\partial \widehat{t}},$$

(4)

where $K_B$ is the Boltzmann constant, e is the electron charge and $\widehat{t}$ denotes the tangent vector to the magnetic field lines.

Integrating Equation (4), it is possible to obtain:

$$\varphi -\frac{K_B{T}_e}{e}ln\left(n_e\right)={\varphi }_0-\frac{K_B{T}_{e,0}}{e}ln\left(n_{e,0}\right)={\varphi }^{*}\left(\lambda \right),$$

(5)

where the RHS term (${\varphi }^{*}$) is known in the literature as thermalized potential [25], and it is constant along a magnetic field line. The thermalized potential assumption allows to reduce the dimension of the problem to quasi-one-dimensional (Q1D); thus, the fluid domain is mono-dimensional. This means that both the potential and electron temperature are constant along the magnetic field streamlines ($\lambda$), but they vary with them, i.e., $T_e=T_e\left(\lambda \right),\varphi =\varphi \left(\lambda \right)$.

Electron diffusion across the magnetic field is assumed to obey a generalized Ohm’s Law:

$$u_{e,\widehat{n}}={\mu }_{\widehat{n}}(\frac{\partial \varphi }{\partial \widehat{n}}-\frac{K_B{T}_e}{e{n}_e}\frac{\partial n_e}{\partial \widehat{n}}-\frac{K_B}{e}\frac{\partial T_e}{\partial \widehat{n}}),$$

(6)

where $\widehat{n}$ denotes the normal vector to magnetic field lines and ${\mu }_{\widehat{n}}$ is the electron mobility across the magnetic field line.

By knowing the relationship between the derivatives along $\widehat{n}$ and the derivatives with respect to $\lambda$ (Equation (A7)) and by utilizing the definition of thermalized potential Equation (5) and Equation (6) can be finally written as:

$$u_{e,\widehat{n}}={\mu }_{\widehat{n}}[-rB\frac{\partial {\varphi }^{*}}{\partial \lambda }-rB\frac{K_B}{e}(ln\left(n_e\right)-1)\frac{\partial T_e}{\partial \lambda }],$$

(7)

where B is the magnetic field intensity. Equation (7) states that the electron velocity across the magnetic streamlines depends on electron temperature and thermalized potential; hence, it is function of $\lambda$.

The electron mobility is modelled as a combination between the classical and anomalous contributions, accounting for plasma turbulence and collisions with neutrals, ions and walls:

$${\mu }_{\widehat{n}}=\frac{m_e({\nu }_{ei}+{\nu }_{en}+{\nu }_w)}{e{B}^2}+k\frac{1}{16B}$$

(8)

Here, $m_e$ is the electron mass; ${\nu }_{ei}$ and ${\nu }_{en}$ are the electron–ion/neutral collision frequencies; ${\nu }_w$ is the electron–wall collision frequency. The coefficient k is chosen empirically; generally, it ranges from 0 to 0.1 in the channel and it assumes a value equal or higher than 1 outside [16,17,18,23]. The wall collision frequency ranges between $0.1\times {10}^7$ s${}^{-1}$ and $1\times {10}^7$ s${}^{-1}$, according to the model proposed by Haagelar et al. [18]. Since mobility depends on empirical parameters, the model is not self consistent.

2.2.2. Current Conservation

To conserve current:

$$I_a=I_i+I_e$$

(9)

where $I_a$ is the current collected by the anode while $I_e$ and $I_i$ are the electron and ion current, respectively. Near-wall current was not included and its contribution has been regarded by varying opportunely the coefficients of electron mobility mentioned above.

Equation (9) can be rewritten considering a surface obtained rotating magnetic streamline ($\lambda$) on the thruster axis:

$$I_a=-2\pi e{\int }_{\lambda }{n}_i{u}_{i,\widehat{n}}rds+2\pi e{\int }_{\lambda }{n}_e{u}_{e,\widehat{n}}rds,$$

(10)

thus, the integrals are performed along $\lambda$. Here $u_{i,\widehat{n}}$ and $u_{e,\widehat{n}}$ are the ion and electron velocity normal to magnetic streamline whereas r is the radial coordinate.

Expanding the electron velocity normal to $\lambda$ (Equation (7)) and simplifying, it is possible to obtain the derivative of thermalized potential:

$$\frac{\partial {\varphi }^{*}}{\partial \lambda }=\frac{-I_a-2\pi K_B\frac{\partial T_e}{\partial \lambda }{\int }_{\lambda }{n}_e{\mu }_{\widehat{n}}B(ln\left(n_e\right)-1)r^{2}ds-2\pi e{\int }_{\lambda }{n}_i{u}_{i,\widehat{n}}rds}{2\pi e{\int }_{\lambda }{n}_e{\mu }_{e,\widehat{n}}B{r}^{2}ds}$$

(11)

2.2.3. Electron Energy Equation

The electron energy equation can be written as:

$$\frac{\partial }{\partial t}\left(\frac{3}{2}{n}_e{K}_B{T}_e\right)+\nabla ·\left(\frac{5}{2}{n}_e{K}_B{T}_e{\overrightarrow{u}}_e+\overrightarrow{q_e}\right)=S_j-S_i,$$

(12)

where the RHS of Equation (12) takes into account the elastic and inelastic collision terms, and $\overrightarrow{q_e}$ represent the heat conduction.

Particularly,

$$\overrightarrow{q_e}=K_e\nabla T_e=\left(\frac{5n_e{K}_B^2{T}_e{\mu }_{\widehat{n}}}{2e}\right)\nabla T_e,$$

(13)

where $K_e$ is the heat conduction coefficient found with the assumption that heat and mass diffuse in the same way.

The term $S_j$ is the joule heating due to elastic collision and can be expressed as:

$$S_j=\overrightarrow{j_e}·\overrightarrow{E},$$

(14)

where $\overrightarrow{j_e}=e{n}_e{u}_{\widehat{n}}$ is the electron flux.

The inelastic collision term, $S_i$ represents the energy loss due to ionization and excitation of neutrals and ions. The energy loss due to the excitation of neutrals is significant and can not be neglected, while the excitation of ions is at least two orders of magnitude smaller and it is not considered. In addition, since double-ionized ions are not modelled, their contribution on $S_i$ is disregarded.

The excitation of neutrals can be modelled increasing the ionization energy with the energy associated with the excitation of neutrals:

$$\phi =E_i+\frac{{\sum }_j〈Q_j\left(u_e\right)u_e〉E_j}{〈Q^{+}\left(u_e\right)u_e〉}$$

(15)

Here, $E_i$ is the ionization energy, $Q_j$ is the cross section of excitation from the ground state to the energy state $E_j$ and $Q^{+}$ is the cross section of the ground state ionization. An analytic solution of Equation (15) is given by Dugan et al. [32] and it can be fitted [13] as:

$${\phi }^{\prime }=A{e}^{-\frac{B}{z}}+C,$$

(16)

where A, B and C are coefficients related to the gas considered while ${\phi }^{\prime }$ and z are respectively the normalized ion production cost and electron temperature:

$${\phi }^{\prime }=\frac{\phi }{E_i}$$

(17)

$$z=\frac{T_e}{E_i}$$

(18)

Finally, using Equations (15) and (16), is possible to express the inelastic collision term:

$$S_i={\dot{n}}_i{\phi }^{\prime }{E}_i$$

(19)

2.2.4. Boundary Conditions

Equation (11) is integrated to compute thermalized potential assuming Dirichlet Boundary Conditions (BC) at both cathode and anode line, ${\varphi }^{*}=const.$.

Equation (12) is integrated to get electron temperature assuming Neumann BC at the anode, $\frac{\partial T_e}{\partial \lambda }=0$ (indeed the relatively small magnitude of the magnetic field at the anode allows to assume an infinitely large diffusion across a magnetic field line) and a Dirichlet BC at the cathode line, $T_e=const.$

The electron energy loss to the wall is:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}exp\left(\frac{e{\varphi }_w}{K_B{T}_e}\right)\left[\left(2K_B{T}_e-e{\varphi }_w\right)-{\delta }_{eff}\left(2K_B{T}_{sec}-e{\varphi }_w\right)\right]\overline{{\Gamma }_e}\hfill & \left(a\right)\hfill \\ 2K_B\left(T_e-T_{sec}\right)\overline{{\Gamma }_e}\hfill & \left(b\right)\hfill \end{array}\right.\end{array}$$

(20)

Here, $\overline{{\Gamma }_e}$ is the electron flux in the bulk of the plasma, i.e., where the plasma is not affected by the sheath; $T_{sec}$ is the temperature of the secondary electrons emitted (∼1 eV); ${\varphi }_w$ is the wall potential, calculated as:

$${\varphi }_w=\frac{-K_B{T}_e}{e}ln\left(\left(1-{\delta }_{eff}\right)exp\left(0.5\right)\frac{{\overline{c}}_e}{4u_b}\right),$$

(21)

where $u_b$ and ${\overline{c}}_e$ are respectively the Bohm and the most probable thermal velocity of electron, whereas the Secondary Electron Emission (SEE) yield (${\delta }_{eff}$), that is the ratio between the emitted secondary electrons and the incident primary electrons, can be written using the Euler’s Gamma Function ($\Gamma$):

$${\delta }_{eff}=\Gamma [2+B]A{\left(\frac{K_B{T}_e}{e}\right)}^B$$

(22)

Here, A and B are two constants that depend on the wall material (e.g., A = 0.123 and B = 0.528 for BoroSil [17]).

When the electron temperature is lower than a certain limit, the break-up temperature (i.e., $T_e<T_{bk}$), the wall potential is negative. This mean that only the most energetic electron passes the sheath, arriving at the wall; thus, the electron flux to the wall is limited in this condition. Equation (20a) is used. When $T_e\ge T_{bk}$, the sheath collapses, the wall potential reaches values near $T_{sec}$ and this filtering action disappears; thus, the wall losses increase dramatically. Equation (20b) is used.

3. Numerical Model

Figure 3 shows the flow chart of the HYPICFLU code. The pre-processing phase includes the generation of both domain and mesh for PIC sub-model and the definition of the magnetic streamlines (Appendix A). After the initialization computation starts. In a typical iteration of the code, at first the equations of heavy particles’ motion are integrated, then the fluid equations are solved while the quantities related to the heavy particles, i.e., neutral and ion density, their velocities and the ionization rate, are kept constant. The governing equations of the heavy species (neutrals and ions) and electrons are integrated over the same time interval but different timescales: typically (cf. [12,13,19,24]), the time step of the fluid equations ($\Delta t_p\sim {10}^{-11}÷{10}^{-10}$ s) is two orders of magnitude smaller than the PIC one ($\Delta t\sim {10}^{-8}$ s).

3.1. Particle-In-Cell for Heavy Particles

The computational domain for the PIC method in the HYPICFLU code is a rectangle that, starting from the anode, includes the channel of the thruster and part of the near plume region.

The area is fully meshed using a uniform Cartesian grid because of the simplicity it offers in its construction, as well as for the sorting and weighting methods and for the computation of gradient fields, which may be required throughout the simulation. Additionally, it can help in avoiding undesired non-physical effects, which may appear in highly deformed meshes related to the calculation of volumetric forces over the macro-particles [5]. Cell size has been selected as a compromise between accuracy and computing efforts. Indeed, by reducing its dimension, the spatial accuracy of the resulting field is improved but more macro-particles are required to reduce numerical noise.

Each heavy particle is characterized by its position, velocity, specific weight, mass and charge. To update position and velocity, the equations of motion (Section 2.1) are integrated using Leap-frog scheme. The time step, $\Delta t$, is such that the distance travelled by a particle must not be greater than the typical length of a cell [24]; thus, the maximum value that can be used is simply calculated subdividing the dimension of cell by the theoretical exhaust velocity (i.e., the maximum value expected in the domain). This guarantees the satisfaction of Courant stability criterion for the Leap-frog Algorithm, because $\Delta t\le \Delta x/v_{max}$ [33].

Particles are not pushed in cylindrical coordinates but in a 3D Cartesian system; after the push, they were rotating back to the r–z plane (for details see [34]). Electric field, required to push ions (Equation (2)) is calculated by a fluid sub-model, interpolated to the grid nodes and successively gathered [35] to particle position.

The ionization process is divided into two different steps: first, an integer number of macro-ions are created in each cell according to the ionization rate defined over the nodes and then mass is removed from each macro neutral in order to ensure the mass conservation [13,24]. Ionization rate (Equation (3)) is computed using the electron temperature provided by fluid sub-model and interpolated to the nodes.

At the end of each time-step, particle velocity and mass are scattered [35] to the grid nodes. This information, along with ionization rate, is used as the input to solve fluid equations.

The recombination of ions on the walls is based on a Monte Carlo algorithm. Neutrals and ions are characterized by two different specific weights; particularly, being that neutral density is higher than ions density, neutral specific weight is higher (about 2 orders of magnitude). When an ion macro-particle impacts the wall, a random number in the range 0–1 is generated and compared to the ratio between ion and neutral specific weights. If this number is less than the ratio (i.e., about 0.01), the ion macro-particle is removed and a neutral is injected. To verify the algorithm, a series of ions with a fixed position and velocity were fired on a wall. In particular, the mass of neutrals created was compared with mass of ions destroyed. The results have shown that the error was below 5.

3.2. Fluid Electrons Equations

Equations (9) and (12) form a determined system with unknowns $T_e$ and $I_a$ that are solved iteratively until a number of pseudo time steps ($\Delta t_p$) to equal one PIC time step ($\Delta t$) has been completed (Figure 3). Then, thermalized potential (${\varphi }^{*}$) is obtained integrating Equation (11); using the definition of thermalized potential (Equation (5)), it is possible to compute the plasma potential ($\varphi$) and consequently, by deriving, the electric field ($\overrightarrow{E}$).

It is possible to verify that Equation (11) can be rewritten as a function of the derivative of electron temperature to $\lambda$:

$$\frac{\partial {\varphi }^{*}}{\partial \lambda }=h_1+h_2\frac{\partial T_e}{\partial \lambda },$$

(23)

where $h_1$ and $h_2$ are functions of such quantities provided by the PIC sub-model and the discharge current (see Appendix B for details). Knowing the electron temperature distribution, it is easy to integrate Equation (23) to get the thermalized potential. Equation (23) can also be manipulated and used to calculate the discharge current, if electron temperature and thermalized potential fields are provided.

Equation (12) is solved by means of a Finite Volume Method (FVM) taking into account that electron temperature can vary only with the magnetic streamlines (Q1D problem) as consequences of the thermalized potential assumption. Thus, Equation (12) is integrated using a volume element delimited by two consecutive magnetic streamlines ($S_1,S_2$) and the up/down boundary of the domain ($S_3,S_4$), Figure 4a, within which electron temperature can be assumed to be constant.

The discretized form of the equation is:

$$\begin{array}{c}{T}_e^{k+1}=T_e^k-\frac{\Delta t_p}{A_1}[A_2{T}_e^k+A_3{T}_{e,S_1}^k+A_4{T}_{e,S_2}^k+A_5{T}_{e,S_1}^k{\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_1}\hfill \\ \hspace{1em}\hspace{1em}+A_6{T}_{e,S_2}^k{\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_2}+A_7{T}_{e,S_1}^k{\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_1}+A_8{T}_{e,S_2}^k{\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_2}+A_9+A_{10}\frac{\partial T_e}{\partial \lambda }\hspace{1em}\hspace{1em}\\ \hfill +A_{11}{\left(\frac{\partial T_e}{\partial \lambda }\right)}^2+A_{12}{T}_e^k\frac{\partial T_e}{\partial \lambda }+{\dot{E}}_w+P_1\zeta \left(T_e^k\right){\phi }^{\prime }\left(T_e^k\right)],\end{array}$$

(24)

where k is the time iteration, $\Delta t_p$ is the time interval and:

$$T_{e,S_1}^k=\frac{T_j^k+T_{j-1}^k}{2}$$

(25)

$$T_{e,S_2}^k=\frac{T_{j+1}^k+T_j^k}{2}$$

(26)

$${\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_1}=\frac{T_j^k-T_{j-1}^k}{\Delta \lambda }$$

(27)

$${\left(\frac{\partial T_e}{\partial \lambda }\right)}_{S_2}=\frac{T_{j+1}^k-T_j^k}{\Delta \lambda }$$

(28)

$$\frac{\partial T_e}{\partial \lambda }=\frac{T_{j+1}^k-T_{j-1}^k}{2\Delta \lambda }$$

(29)

Here, j is the index of the fluid element, which corresponds to the jth magnetic streamline (MS), Figure 4c. The “fluid mesh” is thus defined by the magnetic streamlines, the number of which states the grid refinement level, and the surfaces $S_1$ and $S_2$ correspond to ${\lambda }_{j-1/2}$ and ${\lambda }_{j+1/2}$ (Figure 4c). The distance ($\Delta \lambda$) between two consecutive MSs is set to be constant. By considering for convenience $\lambda =0$ at the anode boundary, the developed algorithm searches in the domain for iso-$\lambda$ lines spaced by $\Delta \lambda$; as requirements, every MS must go with continuity from the lower to upper boundaries of the domain. The lines identified with this method are ${\lambda }_{j-1/2}$ and ${\lambda }_{j+1/2}$ (defined as “lateral”) and the jth magnetic streamline (“central”) lays in the middle. It is worth noting that, with this approach, the number of MSs is not known a-priori but it depends strictly on $\Delta \lambda$.

The factors $A_i$ in Equation (24) are a function of the quantities provided by the PIC sub-model and the discharge current (see Appendix B for details). They are obtained by integrating along magnetic field lines, which are divided into a number of discrete points for this purpose. The volume integrals of the ith fluid element are calculated over the ith streamline (Figure 4b), while the surface integrals over $S_1$ and $S_2$ are calculated over the streamlines $j-\frac{1}{2}$ and $j+\frac{1}{2}$. On all these nodes the data calculated by the solution of fluid equations are available. Thus, an interpolation to PIC mesh is due in order to be used for the next step. It is worth noting that the number of points along a magnetic streamline is selected to be close to that of the PIC mesh along the radial coordinate and the number of magnetic streamlines is about the number of axial nodes of the PIC mesh, in order to reduce interpolation errors.

The last two terms in Equation (24) refer to, respectively, the inelastic losses (Equation (19)) and the wall losses (${\dot{E}}_w$). The wall losses, modelled using Equation (20), are true only if the surfaces $S_3,S_4$ are within the channel. Outside, according to Kawashima et al. [28], it was neglected (${\dot{E}}_w=0$).

Thus, using a Forward Time Centred Space method (FTCS) for the time derivatives and solving the “space” derivatives by using a central scheme, the electron temperature at the pseudo-time step k + 1 can be computed.

It can be verified that Equation (24) can be rewritten as a linearized Burger’s equation with some additional terms [13]. Thus, the convective-diffusive criterion for the stability criterium for FCTS can not be applied directly, but it can be used for an iterative search of stable pseudo time steps. Thus, several runs for identical conditions were executed with different pseudo time steps: the selected value corresponded to the case where the solution was stable (i.e., electron temperature, discharge current and potential as function of pseudo time step converge to a stable value) and it does not vary considerably with inferior time step [13].

3.3. Code Development and Validation

Extensive efforts were undertaken to verify the results of the model. All new segments of the code were tested and compared with hand calculations before the next segment was begun. In this way, the techniques of ”top-down design” and ”step-wise refinement“ were used to make sure bugs were kept to a minimum. For this purpose, every function was tested stand-alone before being included in the model. For each case, experimental or numerical data were selected to compare the results. Disparities between computed and reference data were carefully examined as possibly due to bugs in the code. Particularly, several verification tests of PIC and fluid module stand-alone were conducted but they are not reported here for the sake of brevity.

4. Results

In this section, the code is used to rebuild suitable experimental data and to perform also comparison with some relevant numerical models presented in the literature. Particularly, the simulation of SPT-100 thruster discharge is selected as a crucial validation test case because a lot of experimental and numerical data about the plasma inside and outside the accelerating channel of the thruster are available (cf., e.g., [16,17,19,28,36]).

The following subsections will give the reader details about SPT-100 and they will provide a concise and precise description of the model results, compared to experimental and numerical data, and their interpretation, as well as the conclusions that can be drawn.

4.1. SPT-100 Thruster

The SPT-100 is a Hall-effect thruster designed and built in Russia by the Experimental Design Bureau Fakel which has over 15 years of on-orbit flight heritage on a Russian spacecraft with this thruster. SPT-100 is a 1.35 kW thruster and it has been well developed and proven for satellite applications. A variety of derivatives based on the SPT-100 [37] have been also built and flown. One of them has flown for the Earth–Moon Mission [37], too.

Figure 5 shows a picture of SPT-100 where it is visible the external magnetic coils which generate the magnetic field. Figure 6 shows the performance of SPT-100. For the present study, the design point condition has been selected (black circle in Figure 6). Finally,

Table 1. SPT-100: geometry and nominal operative parameters  (data from [37]).

Parameter	Value	Unit
Outer radius	50	mm
Inner radius	35	mm
Channel length	25	mm
Max. Radial Magnetic Intensity	0.016	T
Mass flow rate (Xe)	5	mg/s
Discharge Voltage	300	V
Discharge Current	4.5	A

reports its dimensions and nominal operative conditions [38].

4.2. SPT-100 Discharge Simulation

This section reports the results of two tests (

Table 2. SPT-100: characteristics of discharge simulation tests.

	Test 95	Test 106
PIC grid (nodes)	25 × 7	49 × 13
Cell dimensions (mm)	0.25	0.125
Fluid grid	22 ($\lambda$) × 7 (nodes)	42 ($\lambda$) × 13 (nodes)
Neutrals specific weigh	${10}^{12}$	$3\times {10}^{11}$
Ions specific weigh	${10}^{11}$	$3\times {10}^9$
Mobility-Bohm diffusion, $k_{in}$	0.035	0.035
Mobility-Bohm diffusion, $k_{in}$	5	5
Mobility elect-wall coll., $\alpha$	0.1	0.1
Anode temperature (K)	750	750
Walls temperature (K)	850	850
Cathode Potential (eV)	10	10
Cathode Electron Temperature (eV)	4	4

) performed to assess that the code is able to reproduce SPT-100 discharge. The results are compared to the available numerical and experimental data.

The tests mainly differ by the grid refinement level: uniform Cartesian grids with square cells and 25 × 7 (coarse grid) and 49 × 13 (fine grid) nodes are respectively used. Cell dimensions are respectively 0.25 and 0.125 mm. Fluid grid resolutions have been selected to be as close as possible to the PIC grid, in order to reduce interpolation errors. The spacings $\Delta \lambda$ between the magnetic streamlines are respectively $1\times {10}^{-6}$ Vm and $5\times {10}^{-7}$ Vm.

Figure 7 shows as an example the two meshes used in Test 106 calculation. Neutral particles are injected into the domain throughout the cells between radial positions 40 and 45 mm (blue line in Figure 7). Fluid grid nodes are the points used to discretize the “lateral” magnetic streamlines (cf. Section 3.2). Particularly. the last raw of nodes are used to represent the cathode.

PIC time step and the pseudo-time step, used to integrate the electron energy equation, are respectively $5\times {10}^{-8}$ s and $2\times {10}^{-10}$ s, for both calculations. The first one guarantees that one macro-ion advances no more than one cell per time step (stability requirement, cf. Section 3.1) in both tests and it is calculated subdividing cell’s dimension by the theoretical exhaust velocity of SPT-100 (∼20 km/s), calculated as ${(2V_{d}e/m_i)}^{1/2}$[1], where $m_i$ is the ion mass (i.e., Xenon). The pseudo-time step is selected to ascertain that the fluid solver converged (cf. Section 3.2) and the solution does not vary considerably with inferior values.

A trade-off analysis was needed to get the specific weights (Table 2) that allowed to smooth numerical “noise” due to the limitations of the PIC method in simulating a continuum of heavy particles. It was found, in accordance with the literature [19,24], that specific weights allowing to get an average number of macro-particles per cell higher than 30 fulfilled this task.

Bohm/SEE contributions to mobility outside and inside the channel have been empirically set. As a first guess, the values reported by Hofer et al. [17], 0.035 and 1, have been used, but a sensitivity analysis, using coarse grid, has shown that, to fit the plume trends, a higher value is needed (i.e., 5). To assess grid effect, these values are retained also in Test 106 using a finer grid. Inside the channel, the mobility model includes an electron–wall collision frequency of $0.1\times {10}^7$ s${}^{-1}$. The electron–neutral collision frequency is evaluated using the cross section formula proposed in Ref. [29].

It is important to review the boundary conditions. Neutrals are reflected on channel walls with a diffuse-reflection scheme and ions impinging on the wall are neutralized and neutrals diffuse according to the same scheme; the anode and channel walls are respectively set to 750 K and 850 K according to Ref. [17]; the cathode, as seen above, is a magnetic streamline that represents the external boundary of the fluid domain. A Dirichlet BC on electron temperature is set at the cathode (4 eV). A linear variation of $T_e$ is assumed downstream of the cathode where the electron temperature at the boundary of the PIC domain is kept fixed at 0.1 eV [13]. Cathode’s potential is set to 10 V and 0 V is applied at the end of domain. A potential drop of 300 V is established between the anode and cathode. A current of about 4.5 A is expected.

4.2.1. Magnetic Field Calculation

SPT-100 magnetic circuit behaviour was reproduced using FEMM software [41]. It has been adjusted until the trend of radial magnetic field along the channel centreline (radial position = 42.5 mm) assumed a Gaussian axial profile in accordance with the literature [17,28] and the value at the exit section was close to that reported in Table 1. Figure 8 shows the contour plot of magnetic field intensity computed by FEMM and interpolated on a PIC coarse grid (25 × 7 nodes). The radial component of the magnetic field, extracted from the solution domain reported in Figure 8, is depicted in Figure 9 and compared with the reference [17] for qualitative comparison (indeed the magnetic circuit models are different). The radial magnetic field assumes a Gaussian axial profile as desired; moreover, the peak at the exit section of the channel, resulting from the present FEMM simulation, differs only by 5% from the value reported in Table 1 and therefore it can be considered acceptable for the scope of the work.

The knowledge of magnetic fields allows for the calculation of the relative stream function (cf. Appendix A), $\lambda$ (Figure 10), and to build fluid meshes (cf. Section 3.2). The cathode’s BC (Table 2) is set to the magnetic streamline Num. 11 in Figure 10.

4.2.2. Initialization

To speed up the calculations, the discharge chamber was populated at first with neutrals by letting them flow from the anode to the outlet; null initial electron temperature and ion density fields were used to disable ionization (see Equation (3)). After about 30 k time steps, required to stabilize neutral dynamics (numerical oscillation must be expected due to PIC method, cf. Section 2.1), non null guesses of electron temperature and ion density fields were set to trigger ionization. About 30 k were required to complete calculations.

4.2.3. Stop Criteria

The system does not converge to a steady-state solution because the physics of plasma in HETs is intrinsically unsteady (cf. Section 1) and there is the numerical “noise” due to the PIC method (cf. Section 2.1). For this reason, solutions are considered complete when the fluctuations reach a regular frequency and amplitude, and they have repeated several periods. Or, if there is no discernable pattern to the fluctuations, the parameters are averaged over a long time scale, and iteration is stopped when averages reach a constant value (quasi steady-state regime).

The number of followed heavy particles (neutrals and ions) and discharge current are selected as global parameters to monitor the simulation and to decide when to stop the calculation. Figure 11 shows these parameters as a function of time step number for Test 106. The sampling frequency is 50 time steps.

Table 3. SPT-100 discharge simulation: time-average discharge current and number of simulated macro-particles (neutrals and ions). The average numbers of particles per cell have been reported.

	Test 95	Test 106
Discharge current (A)	4.8	4.66
Simulated Ion particles (No.)	10,700	33,357
Ion particle (No.) per cell	74	58
Simulated Neutral particles (No.)	6866	22,666
Neutral particle (No.) per cell	48	39

reports time-averaged values of neutral/ion macro-particle number, discharge current and average number of macro-particle per cell (obtained by dividing the number of simulated particles with the total cell numbers). The averages were performed after 20 k time steps.

The selected specific weights of neutral and ion macro-particles allow to get averaged numbers of macro-particles per cell higher than 30, as required to smooth numerical oscillation due to the PIC method. It is also possible to observe that the discharge current seems to converge to the expected value of 4.5 A (Table 1).

4.2.4. Mass Balance

Mass conservation is verified. It is a very important check because the functions that creates ions and eliminates neutrals due to ionization are different. The trends of time-average ion mass flow and neutral mass flow are reported in Figure 12: the first increases toward the exit channel while the second decreases, as expected. Within the channel, at any section, the sum of mass flows must be as close as possible to inlet mass flow. The value are computed averaging solutions spaced by 50 time steps in the quasi steady-state regime.

4.3. Results Analysis

Figure 13 and Figure 14 show the time-averaged contour plots of the variables typically used to represents HET discharge. Qualitatively, the plasma physics described by the plots matches what is expected: neutrals expand within the channel; plasma potential sharply decays at channel’s exit where electron temperature reaches its maximum value; maximum ion density is just upward the exit section of channel and ions move downstream at high speed (order of thousand of km/s) because of the electric field.

The effect of grid refinement is clearly visible:

• Neutral density contour plots are better characterized, particularly near the anode;
• Ion peak tends to increase its value and the field is overall less diffused with respect to the computation with coarse grid, e.g., the plume is better defined;
• Electron temperature decreases more sharply.

The results of Figure 14 are compared to those reported by Hofer et al. [17] who use the code HPHall2 [16] disabling double-charge ion species creation. It is possible to observe that:

• The peak of electron temperature computed by HYPICFLU is slightly higher than that predicted by HPHall2;
• Plasma potential profile in the wall sheath drops very smoothly with respect to HPHall2, due to a different sheath model; this would affect ion velocities toward the wall and consequently its erosion;
• Ion density field is more diffused in the channel and the peak is higher; refining the grid, diffusion could be further lowered; there is a small deviation of ion plume toward the axis, probably due to the magnetic field configuration.

The values along channel centreline (radial position = 42.5 mm) are extracted from the solution domain and compared to the data presented by Hofer et al. [17] (HPHall2 code) and Kawashima et al. [28] (Hyperbolic-Equation System, HES, approach for the computation of magnetized electron fluids, proposed by Komurasaki and Akarawa [14]), Figure 15. Axial ion velocities computed by HYPICFLU are compared to the experimental measurements by means of Laser-Induced Fluorescence (LIF) velocimetry [36], measured at two facility background pressures ($1.5\times {10}^{-5}$ and $5\times {10}^{-5}$ torr).

Neutral density (Figure 15b) at the inlet ($1.1\times {10}^{20}$ m${}^{-3}$) is very similar with the values computed by HPHall2. The value predicted with the HES approach is one order of magnitude less; the expansion trends are similar but the slopes of HYPICFLU curves are a little higher than HPHall2s (i.e., neutral density at exit channel is about $1\times {10}^{18}$ m${}^{-3}$ vs. $1\times {10}^{19}$ m${}^{-3}$).

The electron temperature peaks are over-predicted by ∼2 eV with respect to both reference codes (Figure 15e). Within the channel, the electron temperatures coincides, but outside, in the plume region, the decay trends differ. As seen above, grid refinement leads to sharper decrease of electron temperature. It has been found that varying Bohm constant ($k_{out}$), it is possible to change the slope of curves: lower values lead to a reduction of temperatures’ decay.

The peak of plasma density, located 10 mm downstream the anode, is ∼$1.6\times {10}^{18}$ m${}^{-3}$, matching both HPHall2 and HES predicted values (Figure 15c). After the peak, ion density decays quite rapidly with respect to HPHall2, but similarly to the HES approach; this is attributed to boundary condition effects.

The decay of potential in Test 95 (Figure 15f) well matches HPHall2 data, predicting a maximum electric field of ∼30 kV m. The main difference is in the plume where potential decreases more slowly. Grid refinement leads to higher slope at the exit section (axial electric peak in Test 106 is about ∼40 kV m) but outside the channel potential curve is closer to the HPHall2 one. It is verified that trigging Bohm constants, $k_{in/out}$, it could be possible to modify plasma potential trends as well as, as seen above, electron temperature profiles.

Axial ion velocities computed by HYPICFLU (Figure 15a) continuously increase in the channel, as expected, but values are over-predicted with respect to the experimental data. Adjusting anomalous mobility parameters, it could be possible to reduce differences as shown by Adam et al. [22].

Ionization rates are compared in Figure 15d. HYPICFLU curves are below HPHall2 ones and they have weaker peaks ($6\times {10}^{23}$ vs. $2.4\times {10}^{24}$ m${}^{-3}$ s${}^{-1}$). It is mostly due to the higher neutral expansion and partially (i.e., within the channel) due to ion expansion after the peak.

Figure 16 shows the electron mobility as a function of axial position. The main difference was outside the channel, where two different values of Bohm constant, $k_{out}$, were used (HPHall: 1; HYPICFLU: 5). It is worth noting that HYPICFLU used a simpler wall collision frequency model [18] than HPHall2 and a formula [29] that fit the most recent electron–neutral cross-section for Xenon (thus, affecting directly electron–neutral collision frequency), Figure 17. The electron–ion collision frequency instead is affected by the differences between ion density and electron temperature profiles along channel centreline.

4.3.1. Performance

Thruster performances predicted by the code, compared to numerical and experimental data, are reported in

Table 4. SPT-100 performance: HYPICFLU tests vs. Exp. [37] and Num. [17] data.

	Exp.	Num.	Test 95	Test 106
$I_d$ (A)	4.5	4.5	4.8	4.65
$P_d$ (W)	1350	1350	1440	1395
$u_{ic}$ (km/s)	15.7	17.0	19.8	19.8
${\dot{m}}_{ic}$ (μg/s)	5.29	5.00	3.93	4.19
$I_{sp}$ (s)	1600	1733	2027	2016
T (mN)	83	85	78	83
$\eta$ (%)	50	55	58	59

. They are computed according to the classical formulas [1]:

$$I_{sp}=u_{ic}/g$$

(30)

$$T=u_{ic}·{\dot{m}}_{ic}$$

(31)

$$\eta =\frac{T^2}{2{\dot{m}}_{ic}{P}_d}$$

(32)

$$P_d=V_d·I_d$$

(33)

Here, $u_{ic}$ and ${\dot{m}}_{ic}$ are the mean ion velocity and mass flow rate at the cathode line. g is the gravitational constant at Earth’s surface. Thruster efficiency, $\eta$, is defined as the thrust, T, times the exhaust velocity, $u_{ic}$, divided by the total input power, $P_d$. From the definition of specific impulse, $I_{sp}$, and thrust, it is possible to easily get Equation (32).

Ions’ exhaust velocity at the cathode line is over-predicted; thus, the specific impulse is too. Ion mass flow rate is under-predicted. Consequently, thrust level, being the product between velocity and mass flow rate, is always matched. Thrust efficiency is slightly higher than both numerical and experimental values.

4.3.2. Computing Cost

Here we report some brief considerations on computing cost. About 70% of the computing cost depends on total number of simulated macro-particles, the remaining 30% on integration of electron equations. The operations are sequentially executed on Intel(R) Xeon(R) CPU E5-2640 v4 @ 2.40 GHz. As a result that the primary goal of the study was to test the model, no care was taken in reducing computational time of the code.

Table 5. SPT-100 discharge simulation: computational time.

	Test 95	Test 106
Total Num. of macro-particles	17,566	56,023
Final time step	46,650	40,000
Time-to-final time step (days)	5	16
Time Step/hour	294	116

reports for each test the number of time step per hour, the total number of time steps and total computational time.

5. Conclusions

In this paper, we have described a newly developed code (HYPICFLU) for Hall Thruster plasma modelling. The code implements a model that simulates the temporal and spatial behaviour of Hall Effect Thrusters’ plasma, both inside the accelerating channel and in the near-plume region, by means of a hybrid approach, Particle-In-Cell (PIC) for heavy species (neutral and ions) and fluid for electrons. A two-region electron mobility model has been incorporated, the “anomalous” part of which uses empirical constants that can vary as a continuous function of axial coordinates, respectively $\alpha$ and k, for electron–wall collision frequency and Bohm diffusion/SEE effects. This model has the advantage of including recent developed relations for each relevant phenomena and therefore improving the overall reliability of modelling the plasma behaviour.

In order to verify these aspects, as a cornerstone test, the discharge of SPT-100 has been selected due to the availability of reliable numerical and experimental data. The simulation results seem to be quite in line with experimental data confirming that such a developed model catches the relevant phenomena in Hall Thrusters. The code–code comparison with HPHall2, the NASA-JPL hybrid-code, has shown that HYPICFLU shows a more diffusive behaviour in the plume region, where ions expand more rapidly (in contrast with experimental measurements, too) and that potential decay in the sheath is weaker, affecting ion-to-wall velocity. The best fit of numerical and experimental data has been obtained by using a = 0.1, k = 0.035 inside the channel and k = 5 outside; the outer value of k is higher than the value used by HPHall-2 whereas the inner one is the same set in HPHall-2 simulations. This aspect is interesting and it is worth it to be further explored. The possible causes of over-expansion could be addressed first by the instance to the effect of mobility function’s shape (i.e., by changing k values and how rapidly the curves vary from inner to outer region) and secondly to the effect of near-plume geometry (i.e., by using semicircular area, like the most of codes). These aspects will be the object of future investigations and by means of the experimental testing that will be conducted on the low-power HET developed in CIRA [42].

Other future developments of the model will include the implementation of a more detailed wall sheath model, aimed to better predict ion velocity for wall to wall erosion calculation; the use of the same grid for PIC and fluid modules, built from the magnetic streamlines, in order to avoid interpolation errors. As a result of the high computational cost, an optimization of the code, to increase the convergence speed, is also foreseen.

In the near future, a parametric study to investigate other operative conditions of SPT-100 will be conducted to asses how the mobility parameters have to be changed to match the expected performances and to further assess the code.