Numerical simulations of a magnetically shielded Hall effect thruster with a centrally mounted cathode are performed with an axisymmetric hybrid particle-in-cell/fluid code and are partially validated with experimental data. A full description of the plasma discharge inside the thruster chamber and in the near plume is presented and discussed, with the aim of highlighting those features most dependent on the magnetic configuration and the central cathode. Compared to traditional magnetic configurations, the acceleration region is mainly outside the thruster, whereas high plasma densities and low temperatures are found inside the thruster. Thus, magnetic shielding does not decrease plasma currents to the walls, but reduces significantly the energy fluxes, yielding low heat loads and practically no wall erosion. The injection of neutrals at the central cathode generates a secondary plasma plume that merges with the main one and facilitates much the drift of electrons toward the chamber. Once inside, the magnetic topology is efficient in channeling electron current away from lateral walls. Current and power balances are analyzed to assess performances in detail.
I. INTRODUCTION
Magnetic shielding (MS) of a Hall effect thruster (HET) chamber has been proven an effective technique to limit both wall erosion, due to high-energy ion bombardment, and heat loads, thus enabling the design of the next generation of HETs, featuring enhanced performances and operational lifetimes.1 Due to their recent development, just a few prototypes of MS-HETs have been experimentally tested to date. In the mid to high power range, there are the 4.5 kW BPT-4000,2 the 12.5 kW HERMeS,3 the 6 kW H6MS,4 the 9 kW H9,5 the 20 kW NASA-300MS,6 and SITAEL’s 5 kW HT5k7 and 20 kW HT20k;8 then, in the low-power range (<1 kW) there are the MaSMi-609 and the ISCT-200.10 Yet, a lower number of studies comparing experimental and simulation results have been realized. Relevant ones were performed with the multi-fluid simulation code Hall2De11 for BPT-4000,1 H6,12,13 HERMeS,14,15 H9,16 and MaSMi-60.17
The relatively low numbers of MS-HET prototypes and studies, together with the lack of predictive models of HET discharges (due to the open problems on plasma–wall interaction and electron turbulence, particularly), make uncertain the direct extrapolation of results and trends from one thruster system to another. This situation affects, for instance, the design and development of new optimized electric propulsion architectures such as the direct-drive power concept,18–20 which requires a precise characterization of the thruster performances and the cathode–anode electrical coupling over the nominal operation range for its appropriate integration with other subsystems.20 Therefore, advances on the validation of simulation tools against experimental data, capable not only of providing a full characterization of a HET plasma discharge and performances, but also of addressing thermal, electrical, and material issues related to its operation, are of central interest in HET research.
This work presents 2D numerical simulations of the 5 kW HT5k,7 and a partial validation of the numerical model with the limited experimental data existing for this prototype. The simulations are carried out with the code HYPHEN21,22 and constitute its first test with an MS topology. HYPHEN is a multi-thruster, hybrid-formulation code, which uses a particle-in-cell (PIC) model for heavy species and a magnetized drift-diffusion fluid model for electrons. Contrary to other hybrid codes, relying on a quasi-1D electron model,23–25 whose application to MS topologies is rather complicated, HYPHEN, as Hall2De, adopts a full-2D electron fluid model on a magnetic field aligned mesh (MFAM), thus allowing a complete characterization of the electron currents, which are crucial to study, for instance, plasma–wall interaction effects. The version of HYPHEN used here is also the first one incorporating a “wall” cathode instead of a “volumetric” cathode. This last one worked fine for HETs with laterally located cathodes (except for the code being axisymmetric, although electron emission becomes quickly homogeneous azimuthally26); however, it was not very suitable for the centrally mounted cathode of the HT5k.
HYPHEN simulates the slow-dynamics, axisymmetric transport of the plasma discharge. This implies that two processes, involving kinetic, non-symmetric, and high-frequency aspects, need to be modeled phenomenologically in the electron fluid equations. The first one is the estimation of the particle and energy fluxes to the thruster dielectric walls. These depend on the electron velocity distribution function (VDF), which is non-Maxwellian because of the partial depletion of the collected VDF’s high-energy tail and the secondary electron emission (SEE) by the walls. Kinetic studies27–29 and experimental evidence30,31 are, therefore, used to tune the fluxes to the walls. The second phenomenological model intends to reproduce the slow, turbulent transport of the electron fluid, resultant from averaging (on the azimuth and the high-frequency time scale) the oscillating, azimuthal forces emanating from nonlinear instabilities.32–36 Since there is not yet an established theory of this averaged turbulence force, authors have treated it as an anomalous collisional force, calibrated, when possible, with experimental data. Still, simulations differ much on the selected strength and shape of the anomalous collision frequency.15,23,25,37–40
The central goal of this work is the analysis of the 2D profiles of the discharge, the identification of the main aspects related to the MS topology and the central cathode, and the effects on thruster performances. The discussion aims at improving the understanding of the plasma physics peculiar to an MS-HET and revealing its potential advantages over a HET with a conventional magnetic topology. The document is organized as follows. Section II describes the HT5k thruster unit and the experimental data supporting this work. Section III describes the simulation model and the main settings. Subsection IV A discusses the tuning of the turbulent parameters to match the experimental data. Then, Subsections IV B and IV C analyze the 2D plasma discharge and the thruster performances. Finally, conclusions are drawn in Sec. V.
II. THRUSTER DESCRIPTION AND EXPERIMENTAL DATA
SITAEL began to develop the HT5k [Fig. 1(a)] in 2013. This 5 kW-class thruster consists of two main elements: the thruster itself and the HC20 hollow cathode. Along the last few years, SITAEL designed, manufactured, and tested three different development models (DMs) of the HT5k, as part of several national and international programmes. The different prototypes underwent several technical investigations,7,24 which permitted to demonstrate low erosion, high performance, direct-drive operations, as well as performance stability in high-vacuum conditions (pressure mbar). The design and manufacturing process of the HT5k thruster unit engineering qualification model started in 2019, when the thruster was chosen as the main propulsive unit for the orbit raising and station keeping of the Ital-GovSatCom geostationary platform.41 The HT5k DM3 prototype considered here implements a centrally mounted cathode and a non-conventional magnetic field topology, which is beneficial from the point of view of discharge channel erosion. Design and previous testing efforts were dedicated to enhance critical components and to optimize the thruster thermal behavior. The thruster tests were run in SITAEL’s IV10 facility, reaching pressures of the order of mbar (Xe) while firing at 4.4 kW of discharge power. The DM3 demonstrated competitive performance and showed stable and efficient operation in the 3–7 kW discharge power range, featuring anodic thrust efficiencies up to .
(a) HT5k thruster unit DM3 in operation. (b) HT5k thruster chamber and near plume geometry. (c) Sketch of the RLC filter unit implemented between the thruster anode wall and the cathode.
(a) HT5k thruster unit DM3 in operation. (b) HT5k thruster chamber and near plume geometry. (c) Sketch of the RLC filter unit implemented between the thruster anode wall and the cathode.
Figure 1(b) displays a sketch of the geometry of the thruster chamber and the near plume region. The geometrical parameters and correspond to the thruster chamber length and width, respectively. The plasma domain to be simulated with HYPHEN corresponds to the cylindrical axisymmetric half meridian plane, including the annular thruster chamber and the near plume region. The latter extends axially from the chamber exit plane up to and radially from the symmetry axis up to . The position of the cathode exit plane, at the thruster symmetry axis, is also indicated in Fig. 1(b). Figure 1(c) shows a scheme of the resistor-inductor-capacitor (RLC) filter unit connecting anode and cathode including the power source voltage . The discharge voltage is set between the anode and the cathode. The discharge current flowing between anode and cathode is indicated in the sketch according to the flow direction of electrons. The values of the RLC filter elements are R = 4.7 , L = 360 H, and C = 94 F.
The experimental setup is detailed in Refs. 24 and 42. The experimental data used in the simulations, with xenon as propellant, are listed in Table I. Time-averaged values of the discharge current and the thrust are available for five operation points, hereafter referred to as cases 1–5, defined by a pair (, ), where is the power source voltage and is the propellant mass flow injected through the anode to the thruster chamber. The former ranges from 300 to 400 V and the latter from 10 to 14 mg/s. For all cases, a 7.5% of is injected through the cathode. Data repeatability (i.e., standard deviation) of the measurements is 5%. The value of the time-averaged ranges from 1900s to 2100s from case to case. In addition, the operation point (, ) = (300 V, 14 mg/s) features an anodic thrust efficiency of 58.2%.
Experimental data for the discharge current Id and the thrust F for the five operation points under consideration. The background pressure is equal to 1.1 × 10−5 mbar for all cases.
Case . | Vs (V) . | (mg/s) . | Id (A) . | F (mN) . |
---|---|---|---|---|
1 | 300 | 14 | 14.6 | 269 |
2 | 400 | 14 | 14.2 | 308 |
3 | 300 | 10 | 10.3 | 184 |
4 | 350 | 10 | 10.1 | 197 |
5 | 400 | 10 | 9.6 | 208 |
Case . | Vs (V) . | (mg/s) . | Id (A) . | F (mN) . |
---|---|---|---|---|
1 | 300 | 14 | 14.6 | 269 |
2 | 400 | 14 | 14.2 | 308 |
3 | 300 | 10 | 10.3 | 184 |
4 | 350 | 10 | 10.1 | 197 |
5 | 400 | 10 | 9.6 | 208 |
III. SIMULATION MODEL AND SETTINGS
A. Simulation model
Figure 2(a) shows a schematic representation of the HYPHEN structure and simulation loop, which is briefly outlined next. HYPHEN is an axisymmetric, hybrid, OpenMP-parallelized code built modularly. The code version for HET simulations consists of three main modules: the Ion module (I-module), which follows a Lagrangian approach for simulating the dynamics of the PIC macroparticles of heavy species; the Electron module (E-module), which solves a fluid model for the magnetized electron population; and the Sheath module (S-module), which provides the proper coupling between the quasineutral plasma bulk, and the thruster walls. The E-module assures automatically plasma quasineutrality in the simulation domain. Thus, the Debye sheaths managed by the S-module are, in fact, discontinuity surfaces adjacent to the thruster walls. The three modules are coupled within a time-marching sequential loop.
(a) Simplified description of HYPHEN time-integration loop. (b) Cylindrical mesh used by the I-module. The red, green, blue, and magenta lines indicate the thruster dielectric walls, the anode, the downstream boundary, and the symmetry axis, respectively. The centrally mounted cathode is indicated by the small black box. (c) The MFAM used by the E-module. Blue and red lines are -parallel and -perpendicular lines, respectively, defining the cells. (d) 2D map of , particularized for case 2.
(a) Simplified description of HYPHEN time-integration loop. (b) Cylindrical mesh used by the I-module. The red, green, blue, and magenta lines indicate the thruster dielectric walls, the anode, the downstream boundary, and the symmetry axis, respectively. The centrally mounted cathode is indicated by the small black box. (c) The MFAM used by the E-module. Blue and red lines are -parallel and -perpendicular lines, respectively, defining the cells. (d) 2D map of , particularized for case 2.
The I-module operates on a structured mesh of the simulation domain, shown in Fig. 2(b). On the contrary, and in order to limit the numerical diffusion arising from the strong anisotropic transport on magnetized electrons, the E-module uses an unstructured MFAM,43 defined by the externally applied magnetic field and shown in Fig. 2(c). The magnetic configuration of the MS-HT5k, which features a null magnetic point inside the channel, attempts to screen well all its internal walls.
Two reference frames are considered: one is the cylindrical frame , with coordinates ; and the second one is the magnetically aligned frame , with and , and coordinates . The orthogonal magnetic coordinates and of the MFAM, Fig. 2(c), are obtained from solving and .
Let , , , and be the charge number, particle density, macroscopic velocity, and current density of the plasma species (i.e. electrons , neutrals , singly charged ions , and doubly charged ions ); be the electric field, with being the electric potential; and being the electron temperature. Every simulation step, the I-module takes as inputs , , and and performs the following tasks: (i) the propagation of macroparticles one time step forward, according to the electromagnetic fields acting on them; (ii) the injection of new macroparticles into the domain and the removal of exiting ones; (iii) the interaction of macroparticles with the thruster walls, such as neutral reflection and ion recombination; (iv) the generation of new ion macroparticles due to the ionization of neutrals; and (v) the computation, through a particle-to-mesh weighting process, of the macroscopic properties characterizing each heavy species. Further details can be found in Refs. 21, 44, and 45.
The E-module, taking these heavy-species magnitudes as inputs, solves a quasineutral, drift-diffusion fluid model for the magnetized electron population, obtaining , , , and the electron heat flux vector . The electron fluid model equations are37,46,47
Equations (1) and (2), with , correspond to plasma quasineutrality and the plasma current conservation equation, and the right hand sides are inputs from the I-module. The momentum equation (3) neglects electron inertia, assumes the pressure tensor to be isotropic, and includes the resistive force and the turbulent (or anomalous) force . The resistive force satisfies
with being the electron elementary mass, being the total momentum transfer frequency due to collisions with all heavy species, being the individual contributions for each heavy species , and
is an equivalent heavy species collisional current density.21 The turbulent force , accounting for azimuth-averaged, wave-based anomalous transport, is modeled phenomenologically as23,37,38,48
with being the turbulent collision frequency, being the electron gyrofrequency, and being the phenomenological function representing the local turbulence level.23,37 Turbulent-based contributions to the axial and radial momentum equations are negligible compared to the rest of forces there.
The main assumption of the electron drift-diffusion model is to neglect inertia in the momentum equation. This is justified as long as the electron kinetic energy is much less than their thermal energy, i.e., , a condition well satisfied in HET discharges, except in localized regions and certain operation points. The great advantage of the drift-diffusion model is that the momentum equation (3) reduces to a generalized Ohm’s law for (with an electric conductivity tensor ), which is much easier to treat numerically than the whole differential equation on .
The component of Ohm’s law along the local cross-field direction is49
where is the parallel electric conductivity, is the classical Hall parameter, and is the reduced Hall parameter when including . Therefore, the effective Hall parameter is and scales as if turbulent transport dominates. The last term on the right side of Eq. (9) represents collisions with heavy species and is negligible generally.
Equation (4) is the electron energy equation for an isotropic pressure tensor in the inertialess limit, where the second term in the left side gathers the enthalpy and heat fluxes, and the right side includes the work of the electric field and the power losses from inelastic collisions (e.g., excitation and ionization). Fourier’s law for the heat flux, Eq. (5), includes the thermal conductivity tensor , and corresponds to the drift-diffusion limit of the evolution equation for .46
The numerical treatment of Eqs. (2)–(5) in the unstructured, irregular MFAM was developed in Refs. 43, 49, and 50. As shown in Fig. 2(c), the MFAM is composed of inner and boundary cells. Inner cells are those enclosed by -parallel (blue) and -perpendicular (red) lines. Boundary cells, however, contain at least one boundary face aligned with a domain boundary, which is not a magnetic line generally. Centroids (or computational points) of both cell and faces correspond to the magnetic center or the geometric center, when the former is not available (e.g., at boundary cells). A finite volume scheme on the MFAM cells is applied to conservation Eqs. (2) and (4); and gradient reconstruction schemes, ad hoc for the unstructured MFAM, are applied to Ohm’s and Fourier’s vector laws (3) and (5). For each variable and , this generates one algebraic equation per MFAM cell. Additionally, boundary conditions specifying and (and discussed below) are imposed at each MFAM boundary face, with being the outward unit normal vector. This yields matrix equations for both and at the centroids of the cells and of the boundary faces of the MFAM. A direct solver for sparse linear systems is used for the parallelized computation of the solution.51,52
The reference for the potential is set at the cathode boundary faces, so the anode potential is . The RLC filter unit, Fig. 1(c), relates to the imposed source voltage, , and the varying discharge current, , through
This equation is integrated within the E-module over time with a first order numerical scheme.
The simulation domain extends up to the sheath edge of the quasineutral plasma, represented by the MFAM boundary faces. The solution for the (infinitely) thin Debye sheaths is detailed in the Appendix. There, Eqs. (A2) and (A7) provide the appropriate conditions at each MFAM boundary face in the form of nonlinear relations for and vs the potential jump across the sheath, .
At a dielectric wall, expressions for and are quite straightforward: a zero-collected electric current yields directly ; Eq. (A2) is then solved for , and Eq. (A7) yields . At the current-driving anode, with known potential , the determination of and , at each anode face of the MFAM, requires to compute previously in the following way: Eq. (A2) is combined with Ohm’s law (3) for yielding a non-linear implicit equation for at each anode face. These equations are linearized and introduced in the matrix system yielding all ; if needed, iterations on the implicit equations are run.
At the cathode boundary faces, the discharge current divided by the cathode area defines the electron current density . The electron energy flux at each cathode face, expressed as , is set equal to , with being the average emission energy per electron. At the (quasineutral) axis, symmetry conditions imply and ; since the I-module yields , one has and too. Then, at the plume downstream (quasineutral) boundary, a current-free condition and a Maxwellian electron heat flux are imposed.49,50
The time discretization of the electron equations follows a semi-implicit scheme,53 with a sub-time step and . This scheme allows to keep a linear system for while reducing the value of required for convergence. Finally, note that mesh interpolation of plasma variables between ion and electron modules is required.
B. Simulation settings
Figures 2(b) and 2(c) show the PIC mesh and the MFAM used to simulate the HT5k. The main characteristics of the meshes and the relevant simulation parameters are listed in Table II. A singular magnetic point () is located inside the thruster chamber. The green line at the left boundary in Fig. 2(b) represents the annular anode wall, while the small black box in Fig. 2(c) indicates the position of the central cathode boundary.
Main simulation parameters and mesh characteristics.
Simulation parameter . | Units . | Value . |
---|---|---|
PIC mesh number of cells, nodes | … | 2582, 2696 |
PIC mesh smallest grid size | mm | 1 |
MFAM number of cells, faces | … | 2848, 5830 |
MFAM average skewness60 | … | 0.058 |
Ion-moving time step, Δt | ns | 15 |
Total number of simulation steps | … | 60 000 |
Injected Xe velocity | m s−1 | 300 |
Injected Xe temperature | K | 850 |
Simulation parameter . | Units . | Value . |
---|---|---|
PIC mesh number of cells, nodes | … | 2582, 2696 |
PIC mesh smallest grid size | mm | 1 |
MFAM number of cells, faces | … | 2848, 5830 |
MFAM average skewness60 | … | 0.058 |
Ion-moving time step, Δt | ns | 15 |
Total number of simulation steps | … | 60 000 |
Injected Xe velocity | m s−1 | 300 |
Injected Xe temperature | K | 850 |
The prescribed xenon mass flow is injected from a Maxwellian reservoir through the whole annular anode featuring a flat profile with a sonic axial velocity based on its own temperature (Table II). Considering the same injection properties, a neutral mass flow is injected through the cathode boundary together with an electron current equal to . The emission energy of electrons at the cathode is set to eV;54 just as a sensitivity check, simulations run with eV have shown no observable differences in the discharge, except naturally at the cathode neighborhood.
Wall recombination of ions contributes to the neutral density. Singly and doubly charged ions are generated volumetrically by electron-neutral collisions. Single ionization rates are obtained from the BIAGI database,55 while double ionization rates follow the Drawin model,56 including the reactions and . Neutrals from ion recombination at walls are re-emitted diffusely considering complete ion energy accommodation at the wall, as suggested by several authors.57,58 Thus, the neutral emission energy is only given by the wall temperature, which is set to 850 K.59 Neutrals are reflected diffusely at the wall with zero energy accommodation; Refs. 21 and 45 provide further details on the interaction of heavy-species macroparticles with walls.
The simulations monitor independently the populations of neutrals, singly charged, and doubly charged ions. Each species population is controlled, setting a target number of 500 macroparticles per cell with a 10% of tolerance.21
The (ion) time step in Table II is set so that a typical doubly charged ion takes at least two time steps to cross the smallest PIC cell. The simulations are started by injecting neutrals through the anode and the cathode and considering a minimum background plasma density to trigger the discharge.45 Every simulation features a total of 60 000 time steps (equivalent to 900 s of simulation time) so that undergoes a sufficiently large number of low-frequency (i.e., breathing mode) oscillation cycles. Five sub-time steps per ion time step () are used to integrate electron equations.21 All the results shown in Secs. IV A–IV C are time-averaged over several cycles.
IV. SIMULATION RESULTS AND DISCUSSION
A. Fitting of turbulence parameters
In order to complete the electron model, the turbulence function in Eq. (8) must be chosen. Since for each case in Table I only two experimental parameters, and , are known, the function will be of the axial “step-out” type shown in Fig. 2(d), with two fitting parameters only, and , applying, respectively, and approximately, inside and outside the channel. These step-out profiles have provided good fittings in previous studies.39,58,61,62 Hereafter, a particular step-out profile is referred to as (, ).
For each of the cases of Table I, the pair (, ) has been tuned to reproduce with a relative error smaller than 5%, thus consistent with the repeatability of the experimental data. The results are in Table III and Fig. 3. Turbulent transport is larger in the plume by nearly one order of magnitude, in line with the existing literature. Also aligned with previous studies for MS-HETs13 and traditional HETs,62 there is a moderate change of the turbulence parameters with the operation point. The parameter features the largest variation, increasing as decreases and increases. The obtained turbulence fitting for turns out to be also good reproducing the main oscillations of , those known as the breathing mode34,48,63,64 and listed in Table III. For instance, for case 1, simulations yield an oscillation frequency kHz, close to the 22 kHz reported in experiments. For all cases, ranges between 4% and 13% of , also consistent with experiments.
Comparison of experimental (error bars) and simulation results (empty markers) for (a) the I-V curve and (b) the F-V curve. The black circle and blue square markers correspond to mg/s (cases 1 and 2) and 10 mg/s (cases 3–5), respectively. The red triangle marker corresponds to case 2* in Table III.
Comparison of experimental (error bars) and simulation results (empty markers) for (a) the I-V curve and (b) the F-V curve. The black circle and blue square markers correspond to mg/s (cases 1 and 2) and 10 mg/s (cases 3–5), respectively. The red triangle marker corresponds to case 2* in Table III.
Simulation results for the best fit of the turbulence parameters (4th column), except case 2*. Simulated results for Id and F (5th and 6th columns) are within a 5% error of the values in Table I. Frequency and relative half-amplitude of oscillation of Id are listed in 7th and 8th columns.
. | Vs . | . | (αt1, αt2) . | Id . | F . | fd . | ΔId/Id . |
---|---|---|---|---|---|---|---|
Case . | (V) . | (mg/s) . | () . | (A) . | (mN) . | (kHz) . | (%) . |
1 | 300 | 14 | (0.8, 8.0) | 15.0 | 276 | 20.1 | ±13.0 |
2 | 400 | 14 | (0.7, 5.0) | 14.8 | 319 | 23.7 | ±6.8 |
2* | 400 | 14 | (0.8, 8.0) | 16.7 | 335 | 26.2 | ± 7.2 |
3 | 300 | 10 | (0.8, 7.0) | 9.8 | 187 | 17.8 | ±7.1 |
4 | 350 | 10 | (0.7, 6.0) | 9.8 | 203 | 15.2 | ±4.6 |
5 | 400 | 10 | (0.7, 3.5) | 9.4 | 213 | 18.2 | ±7.4 |
. | Vs . | . | (αt1, αt2) . | Id . | F . | fd . | ΔId/Id . |
---|---|---|---|---|---|---|---|
Case . | (V) . | (mg/s) . | () . | (A) . | (mN) . | (kHz) . | (%) . |
1 | 300 | 14 | (0.8, 8.0) | 15.0 | 276 | 20.1 | ±13.0 |
2 | 400 | 14 | (0.7, 5.0) | 14.8 | 319 | 23.7 | ±6.8 |
2* | 400 | 14 | (0.8, 8.0) | 16.7 | 335 | 26.2 | ± 7.2 |
3 | 300 | 10 | (0.8, 7.0) | 9.8 | 187 | 17.8 | ±7.1 |
4 | 350 | 10 | (0.7, 6.0) | 9.8 | 203 | 15.2 | ±4.6 |
5 | 400 | 10 | (0.7, 3.5) | 9.4 | 213 | 18.2 | ±7.4 |
It is known that the thrust and, mainly, the discharge current are rather sensitive to the turbulence parameters. This is illustrated here by case 2* in Table III, which corresponds to the operation point of case 2 but using the turbulent fitting of case 1. Simulation errors with respect to experimental data are about 18% for and 9% for .
Nonetheless, these differences in performances do not change the main trends of the 2D discharge. Figures 4(a)–4(d) show the 1D axial profiles of , , , and , along the thruster channel midline and for cases 1, 2, and 2*. For case 1, most of the ion acceleration occurs in the near plume along 2–3 channel lengths. This agrees with existing experimental measurements for a previous HT5k prototype.65 In case 2, with higher , the total potential fall spans over a broader axial region, which would explain the slightly higher plume divergence we observe and is in line with previous numerical13 and experimental9,10 studies. The peaks of and are outside the channel and move slightly downstream for higher , in line with previous numerical studies too.13
Time-averaged axial profiles along the thruster chamber midline [(a)–(d)] for cases 1, 2, and 2* and [(e)–(h)] for cases , 2, and with , , , and , respectively, for each column.
Time-averaged axial profiles along the thruster chamber midline [(a)–(d)] for cases 1, 2, and 2* and [(e)–(h)] for cases , 2, and with , , , and , respectively, for each column.
This behavior of follows some experimental results reported in Ref. 10 but differs from those observed in Ref. 66, suggesting that the observed mild trends with depend, at least partially, on the prototype and the operational conditions.
The change of turbulent collisionality from case 2 to 2* in Fig. 4, results mainly in steeper gradients of and in the region where is transitioning from to , but plasma variables behave very similarly for cases 2 and 2*, and typical MS effects are reproduced in both cases, in line with the numerical studies of Refs. 17 and 67.
Figures 4(e)–4(h) show the effect on the discharge of moving the transition point upstream or downstream. For case 1 and at the mid-channel radius, that point was placed at (for reference, the peak of the magnetic field is at ). Cases and displace the transition point to and 1.59, respectively. The variations of are small, within the accepted range of , while the value of remains practically constant. Changes in the 2D discharge inside the chamber are barely perceptible, but mild changes on the location and values of the maximum are observed. This suggests that the location of the transition point can be used to fit the location of the peak of .
The strength of the maximum can be adjusted using a third fitting parameter around the peak of the magnetic field, as it is done in “quenched turbulence models,” in general with .68–70 One stage ahead in fitting techniques is due to Mikellides and Lopez Ortega,15 who adjust a multi-piecewise function with a larger set of experimental plasma data. Interestingly, the resultant function resembles a “step-out plus quenched model,” with .
In any case, all these are just numerical attempts to adjust time-averaged experimental data, without a firm theoretical basis (not even for extrapolating the function from one operation point to another). The fitting problem is more severe for fully 2D MS topologies than for the traditional ones, which feature quasi-radial magnetic lines. Fortunately, this weakness in modeling accurately the slow turbulent transport of electrons does not preclude the study of other central features of the plasma discharge.
B. Analysis of the 2D plasma discharge
Case 1 in Table III is chosen for this analysis. Figure 5 shows the 2D (,) maps of relevant time-averaged plasma properties inside the thruster chamber (left column) and in the whole simulation domain (right column). The neutral density map, Figs. 5(a) and 5(b), illustrates the effects of gas ionization plus the wall-born neutrals from ion recombination. The neutral density exhibits a decrease of about two orders of magnitude inside the thruster chamber, corresponding to a large propellant utilization. In addition, plot (b) shows well the injection and partial ionization of the neutral gas injected through the cathode. Then, the enhanced electron-neutral collisionality favors the coupling of cathode electrons with the main ion beam. The plasma density [Figs. 5(c) and 5(d)] inside the chamber is higher and more uniform than in conventional HETs, due likely to the acceleration region being moved outwards from the chamber. As usually, the plasma density presents its maximum at a central location inside the chamber, here close to the magnetic null point. As shown in Fig. 4, the increase in the peak at higher enhances gas ionization, augmenting slightly the peak of inside the chamber. In contrast, the lower in the near plume at higher is due to the larger ion acceleration. The right plot shows the secondary plasma plume created by ionization of cathode neutrals, which merges downstream with the main plume.
Time-averaged 2D (,) contour maps for case 1. [(a) and (b)] Neutral density , [(c) and (d)] plasma density , [(e) and (f)] electric potential , and [(g) and (h)] electron temperature . The left column plots show magnitudes inside the thruster chamber, while the right column plots correspond to the whole simulation domain. The centrally mounted cathode is indicated by the small black box.
Time-averaged 2D (,) contour maps for case 1. [(a) and (b)] Neutral density , [(c) and (d)] plasma density , [(e) and (f)] electric potential , and [(g) and (h)] electron temperature . The left column plots show magnitudes inside the thruster chamber, while the right column plots correspond to the whole simulation domain. The centrally mounted cathode is indicated by the small black box.
The electric potential [Figs. 5(e) and 5(f)] also presents a maximum close to the singular magnetic point. Potential variations are rather small inside the chamber, as expected in an MS-HET configuration. As commented already, the MS moves the acceleration region outside the chamber and equipotential lines follow approximately the magnetic lines there. Inside, the magnetic null point and the pressure gradients uncouple equipotential lines from magnetic ones. The electron temperature [Figs. 5(g) and 5(h)] peaks in the near plume, rather close to the maximum . Temperature isolines follow closely magnetic lines and temperature gradients are very pronounced when the electron flow enters into the chamber. The temperature is below 5 eV near all chamber walls, which is indeed one of the main achievements of MS topologies12,13,25 and leads to small energy losses to the walls.
Figure 6 depicts 2D vector maps for the longitudinal components of the ion, electron, and electric currents, defined as and so on. The ion current is obtained from a particle-to-mesh weighting algorithm, while the electron current comes out directly from the fluid model. The ion streamlines, Figs. 6(a) and 6(b), reflect the existence of backward, forward, and lateral ion flows. Although there is a point with , notice that the ionization source is distributed in the whole channel volume and the streamlines represent the ion macroscopic behavior. Ions, practically unmagnetized, follow the electric field and are not prevented from impacting the channel walls. Interestingly, plot (a) shows nearly wall-parallel ion current streamlines close to the channel exit, a fact also observed experimentally.5 Ion streamlines moving into the far plume present the expected divergence outside the thruster chamber. Plot (b) shows again how the secondary ion beam created around the cathode combines with the main beam downstream.
Time-averaged 2D (,) contour maps for case 1. Magnitude of the longitudinal [(a) and (b)] ion current density vector , [(c) and (d)] electron current density vector , and [(e) and (f)] electric current density vector . Blue lines with arrows depict the streamlines of [(a) and (b)] , [(c) and (d)] , and [(e) and (f)] . The left column plots show magnitudes inside the thruster chamber, while the right column plots correspond to the whole simulation domain. The centrally mounted cathode is indicated by the small black box close to the axis.
Time-averaged 2D (,) contour maps for case 1. Magnitude of the longitudinal [(a) and (b)] ion current density vector , [(c) and (d)] electron current density vector , and [(e) and (f)] electric current density vector . Blue lines with arrows depict the streamlines of [(a) and (b)] , [(c) and (d)] , and [(e) and (f)] . The left column plots show magnitudes inside the thruster chamber, while the right column plots correspond to the whole simulation domain. The centrally mounted cathode is indicated by the small black box close to the axis.
The longitudinal electron streamlines, Figs. 6(c) and 6(d), are all born from the cathode but they are split into one beam directed to the far downstream region, which neutralizes the main ion beam, and one beam moving into the interior of the thruster. This second beam follows first the magnetic lines, but then enters the chamber almost perpendicular to the magnetic field, and finally, it is again channeled by the convergent–divergent magnetic lines. Interestingly, this last configuration does not constitute at all a magnetic nozzle71 (the magnetic strength is null at the “throat,” instead of maximum) but it channels the electrons anyway. Still, some electron streamlines are driven into lateral walls to cancel the ion flows there. It has also been realized that the neutral injection through the cathode enhances much the diffusive cross-field transport of the inward beam and the electric coupling of the cathode with the main beam. Finally, Figs. 6(e) and 6(f) show the 2D streamlines of the longitudinal electric current. The zero net collected current condition imposed at the dielectric walls and along the downstream boundary (which is assumed current-free locally) yields the expected current loops connecting anode and cathode.
To complete the description of the discharge, Figs. 7(a)–7(c) plot the main plasma magnitudes along the thruster internal walls. The abscissa length runs from the inner corner at the exit of the chamber toward the anode corner, continues along the anode, and finishes along the outer, interior wall of the chamber. Figure 7(a) shows the electron temperature at the sheath edge and the sheath potential fall, , this last one computed from the S-module. The temperature near all internal walls is rather low (3–4 eV) and the ratio is between 1 and 3. This ratio is affected by the SEE and the partially depleted VDF of primary electrons, at the dielectric walls, and by the need to control the local flux of electrons, at the anode. The low makes the SEE yield small: it is 0.2 for eV.
Time-averaged simulation results for case 1 in Ref. 3. Coordinate runs along the thruster chamber walls. Profiles of (a) potential fall across the sheath edge, ; (b) ion, , and electron, , current normal to the walls; and (c) ion, , and electron, , wall-impact energy.
Time-averaged simulation results for case 1 in Ref. 3. Coordinate runs along the thruster chamber walls. Profiles of (a) potential fall across the sheath edge, ; (b) ion, , and electron, , current normal to the walls; and (c) ion, , and electron, , wall-impact energy.
Figure 7(b) plots the electron and ion currents toward the walls (which are constant across the Debye sheaths). The plotted electron current density is the net one, i.e., the difference between the currents of primary and secondary electrons at the wall. Ion and electron current densities are identical at dielectric walls. At the anode, the backward ion current density is around a 20% of the electron current, which is a percentage larger than desirable. Figure 7(c) plots the average ion and electron energy at the wall, per net particle impacting the wall, computed from Eqs. (A5) and (A8). The low values of yield ion-impact energies, , well below typical threshold values for wall erosion,13,72 which is the main advantage expected from MS topologies.
C. Current and power balances
The ion current balance at steady state can be expressed as
where is the current of ions generated by ionization in the simulation domain; , , and are the ion currents impacting the dielectric, anode, and cathode walls, respectively [and defined in Fig. 2(b)]; and is the ion beam current leaving the domain at plume boundaries, the only one contributing to thrust. All currents are defined as positive; comes out from a volumetric integration, and the other ones are computed from surface integrals at the domain boundaries.
Table IV details for cases 1–5 and how it is distributed among the different boundaries; is about one order of magnitude lower than and has not been included. The table also includes the propellant utilization, the current efficiency, and the charge efficiency, defined as
respectively. Here, and is the total ion mass flow across the plume boundaries. The trends of these partial efficiencies with source voltage and mass flow are the usual ones. Compared to a typical discharge in a conventional HET, and are rather good, and the relative amount of doubly charged ions, measured by , is similar. The relative current losses to the lateral walls, , are similar to conventional HETs, because of the compensation of two trends: a higher plasma density and a lower electron temperature. For this prototype, the relative current losses to the anode, , are high, likely due to the existence of the null point not far away from the anode. Although further analyses would be needed, the conclusion here on the current balance is that the MS topology of the HT5k does not present clear advantages over a more conventional one. However, this is not going to be the conclusion with the power balance.
Value of Iprod and fractions of Iprod corresponding to the different contributions to the current balance in Eq. (11) for cases 1–5. Values of ηu, ηcur and ηch for cases 1–5.
Case . | Vs (V) . | (mg/s) . | Iprod (A) . | Ii∞/Iprod . | IiD/Iprod . | IiA/Iprod . | ηu . | ηcur . | ηch . |
---|---|---|---|---|---|---|---|---|---|
1 | 300 | 14 | 27.6 | 0.42 | 0.39 | 0.18 | 0.94 | 0.77 | 0.90 |
2 | 400 | 14 | 33.0 | 0.36 | 0.42 | 0.21 | 0.94 | 0.78 | 0.89 |
3 | 300 | 10 | 17.4 | 0.45 | 0.37 | 0.17 | 0.91 | 0.79 | 0.92 |
4 | 350 | 10 | 18.6 | 0.42 | 0.38 | 0.19 | 0.90 | 0.79 | 0.91 |
5 | 400 | 10 | 18.1 | 0.44 | 0.37 | 0.18 | 0.92 | 0.85 | 0.85 |
Case . | Vs (V) . | (mg/s) . | Iprod (A) . | Ii∞/Iprod . | IiD/Iprod . | IiA/Iprod . | ηu . | ηcur . | ηch . |
---|---|---|---|---|---|---|---|---|---|
1 | 300 | 14 | 27.6 | 0.42 | 0.39 | 0.18 | 0.94 | 0.77 | 0.90 |
2 | 400 | 14 | 33.0 | 0.36 | 0.42 | 0.21 | 0.94 | 0.78 | 0.89 |
3 | 300 | 10 | 17.4 | 0.45 | 0.37 | 0.17 | 0.91 | 0.79 | 0.92 |
4 | 350 | 10 | 18.6 | 0.42 | 0.38 | 0.19 | 0.90 | 0.79 | 0.91 |
5 | 400 | 10 | 18.1 | 0.44 | 0.37 | 0.18 | 0.92 | 0.85 | 0.85 |
The plasma power balance for the steady state discharge is
where is the total power deposited into the plasma discharge, sum of the discharge power and the net power delivered through cathode electron emission ( amounting to 1%–2% of ); is the plasma energy flow through the plume boundaries; and are the power losses at the dielectric walls and anode, respectively; and corresponds to the power losses due to inelastic (ionization and excitation) collisions. All powers are defined as positive. is obtained from a volumetric integral, comes from a surface integral at the plume boundary, and and are computed from surface integrals at the respective walls (not at the Debye sheath edges). Equations (11) and (13) have also served to check that the numerical errors given by HYPHEN simulations are acceptable (below 2%).
Inelastic losses correspond entirely to electrons and they are roughly proportional to (notice that single and double ion creation has different ionization energies). Energy losses to walls D and A present contributions of ions and electrons, which were discussed in Sec. IV B, and they are proportional to their respective currents , , and . The downstream energy flow corresponds mainly to ions but it also includes the residual electron energy flow in the far plume.
Table V lists the main contributions to the power balance and some related magnitudes for cases 1–5. Observe that while the current losses to the lateral walls amount to about a 40% of , the energy losses to these walls amount to a mere 7%. Adding the energy losses to the anode, the total energy losses to the walls are just 9%–12%, which can be considered an important achievement of MS topologies. Inelastic losses are consistent with the ion production: the ratio yields 23 eV as effective single-ionization cost (including the contribution from excitation collisions).
Value of P and fractions of P corresponding to different contributions to the power balance in Eq. (13) for cases 1–5. Values of η, ηene, ηdiv, and ηdisp for cases 1–5.
. | Vs . | . | P . | . | . | . | . | P∞/P . | . | . |
---|---|---|---|---|---|---|---|---|---|---|
Case . | (V) . | (mg/s) . | (kW) . | η . | Pinel/P . | PD/P . | PA/P . | ( = ηene) . | ηdiv . | ηdisp . |
1 | 300 | 14 | 4.43 | 0.57 | 0.15 | 0.07 | 0.05 | 0.74 | 0.89 | 0.87 |
2 | 350 | 14 | 5.73 | 0.57 | 0.13 | 0.07 | 0.04 | 0.74 | 0.86 | 0.90 |
3 | 300 | 10 | 2.91 | 0.56 | 0.14 | 0.06 | 0.05 | 0.74 | 0.88 | 0.85 |
4 | 350 | 10 | 3.40 | 0.56 | 0.13 | 0.06 | 0.05 | 0.75 | 0.85 | 0.88 |
5 | 400 | 10 | 3.76 | 0.57 | 0.11 | 0.05 | 0.04 | 0.78 | 0.84 | 0.86 |
. | Vs . | . | P . | . | . | . | . | P∞/P . | . | . |
---|---|---|---|---|---|---|---|---|---|---|
Case . | (V) . | (mg/s) . | (kW) . | η . | Pinel/P . | PD/P . | PA/P . | ( = ηene) . | ηdiv . | ηdisp . |
1 | 300 | 14 | 4.43 | 0.57 | 0.15 | 0.07 | 0.05 | 0.74 | 0.89 | 0.87 |
2 | 350 | 14 | 5.73 | 0.57 | 0.13 | 0.07 | 0.04 | 0.74 | 0.86 | 0.90 |
3 | 300 | 10 | 2.91 | 0.56 | 0.14 | 0.06 | 0.05 | 0.74 | 0.88 | 0.85 |
4 | 350 | 10 | 3.40 | 0.56 | 0.13 | 0.06 | 0.05 | 0.75 | 0.85 | 0.88 |
5 | 400 | 10 | 3.76 | 0.57 | 0.11 | 0.05 | 0.04 | 0.78 | 0.84 | 0.86 |
The thrust efficiency is defined and then factorized as
where is the thrust measured from plasma properties at the plume boundary. The energy, divergence, and dispersion efficiencies are defined, respectively, as
with being the flow of axial plasma (ion and electron) energy across the plume boundaries. Here, quantifies the relative power in the downstream plume, assesses the plume divergence based on axial energy and total energy flows, and quantifies the level of velocity dispersion of all plasma species (which would be one for a mono-velocity gas). Plume energy flows include the residual energy of electrons, coming from their incomplete expansion in the finite simulation domain, but this is quite low: about 4% of for case 1.
Setting , the half-divergence angles in these simulations are –. As commented already, plume divergence increases slightly for higher . The corresponding decrease in is compensated by and , thus yielding mild variations of . Previous studies for a different HET prototype also reported slight changes of in the range V.13 Here, we find (or 61% if the “anodic” thrust efficiency is used), the overestimate relative to measured values being just attributed to the slight overestimate of the simulated .
V. CONCLUSIONS
The numerical simulations of the HT5k prototype with HYPHEN have allowed the testing of this code with MS topologies and a centrally located cathode and its partial validation with the limited experimental data available.
Turbulence transport has been modeled with a step-out profile of a phenomenological anomalous frequency, which has been fitted with experimental data of and . The observed trends and the moderate changes of the fitting with the operation point tend to agree with previous works in the literature. Oscillations of the discharge current agree well with experimental ones. A larger set of experimental data would allow a finer tuning of the turbulent transport, but it would probably not alter the present comprehension of the physics of the MS-HET.
A detailed picture of the plasma discharge has been shown for one operation point. The discussion has been focused on those aspects highlighting the effects introduced by the MS topology and the central cathode, over conventional thrusters with lateral cathodes. The MS topology shifts the acceleration region outwards and keeps a high-density, low-temperature plasma inside the chamber. This combination of plasma properties explains that plasma flows to walls are similar to those in conventional HETs. The low plasma temperature in the chamber implies low electric fields, small Debye sheaths around lateral walls, and low SEE, leading to a small energy deposition at the walls and small ion impact energies, well below the usual threshold for material sputtering.
In the near plume, a central cathode with both electron and neutral emission has been simulated. A secondary plasma plume is generated, which merges with the main one downstream. A fraction of the electron emission drifts toward the interior of the thruster perpendicular to the magnetic lines but, once inside, this current is again channeled by the magnetic lines until crossing the null point, except for a small part leaking with ions to the lateral walls.
Finally, detailed experimental measurements of the plasma discharge, in progress, would allow to validate more firmly the present simulation tool, especially those aspects where phenomenological approaches are still important, such as turbulent transport, plasma–wall interaction, and downstream boundary conditions.
ACKNOWLEDGMENTS
This work has been supported by the EDDA project, funded by the European Union’s Horizon 2020 Research and Innovation Program, under Grant Agreement No. 870470.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts of interest to disclose.
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: THE SHEATH MODEL
The S-module of HYPHEN solves a planar, unmagnetized, collisionless, kinetic model of the thin Debye sheaths develop around the walls. The model includes SEE and retains other non-Maxwellian features of the electron VDF.73 The S-module provides the appropriate boundary conditions for the quasineutral electron fluid equations, specifying, at each quasineutral MFAM boundary face (i.e., Debye sheath edge), the perpendicular electron current and heat flux.
Taking as reference a ceramic material with large SEE, the electric current density of primary (p) electrons from the quasineutral plasma into the wall is assumed to follow
which corresponds to a partially depleted, partially reflected Maxwellian VDF.73 Here, is the potential fall in the sheath, , is the fraction of primary electrons reaching the wall but being reflected back elastically by it, and estimates the replenishment fraction of the VDF tail corresponding to impacting electrons. If is the SEE yield, the net, local electron current density to the wall is
Kinetic studies of plasma–wall interaction27,28 show that because of the weak electron collisionality. Simulations results with between 0.05 and 0.30 have revealed small effects on thruster performances and relevant plasma profiles. This is in agreement with experimental results evidencing a reduced influence of plasma–wall interaction on the plasma discharge in magnetically shielded thrusters.30,31 Based on this, is set for all simulations here, leaving further sensitivity analyses or more accurate models to further work.
For a conducting wall, such as the anode, we just take . For the lateral ceramic walls, and are modeled according to Refs. 73–77,
with , , and being material dependent parameters and being the effective upper-bounded SEE yield, corresponding to a space-charge limited (SCL) sheath. For the boron nitride walls of our thruster, the following values were taken: , eV, eV, and .
The local net power density deposited by the whole electron population at the wall is
where is the average energy per net collected electron and is the average energy per wall-emitted electron (equal to 4 eV in the simulations here). Then, the net power density of electrons at the sheath edge is
At that edge, the kinetic electron model inside the sheath is matched with the fluid electron model of the outer quasineutral domain. The matching of is obvious, but for the energy fluxes, the “kinetic” one, of Eq. (A6), is matched to the “fluid” one, . This yields the heat flux at the sheath edge as
which is used as boundary condition on the quasineutral domain for the electron fluid.
For ions, the PIC solution in the quasineutral domain is also matched to the kinetic sheath model at the sheath edge. Dedicated particle-to-surface weighting schemes for the PIC formulation of ions21,78,79 yield directly the net ion power density at (the PIC side of) the sheath edge, . Then, the power density deposited at the wall is and the average energy per wall-impacting ion is