FIG. 3.1. Simplified flowchart for the preliminary computational model.

FIG. 3.2. Illustration of Wirz's predictor/corrector method. The magnetic field value calculated at the midpoint is used to execute the corrected gyromotion.

FIG. 3.3. Magnetic field configuration used in the Discharge Experiment. The field is computed with the analytical equations for permanent magnets.^5,6

FIG. 3.4. The adaptive Cartesian mesh used in the simulation. Grid refinement is provided near the downstream cylindrical magnet.

FIG. 3.5. Error in electric field using the standard finite difference equations for a rectilinear mesh. The error is computed for the cells inside the yellow box. Relatively higher error is seen at the coarse/fine grid interfaces. The convergence of the solution is clearly seen from the error in three different grid levels.

FIG. 3.6. Error in electric field using the least squares quadratic fitting method. The error is computed for the cells inside the yellow box. The error at the coarse/fine grid interfaces is mitigated compared to the case in Fig. 3.3.

FIG. 3.7. Electric potential (V) calculated by the computational model

FIG. 3.8. Primary electron density (m^-3).

FIG. 3.9. Normalized current density profile along the downstream plate radius. The data are normalized to the maximum values for each species.

FIG. 3.10. Examples of ion trajectories very near the cusp.

Table 1. Characteristic radii (mm) for the Cusp Confinement Discharge Experiment.

Recent Accomplishments

The preliminary model has been developed and used to simulate the Cusp Confinement Discharge Experiment for a sparse plasma created by the electron gun. Although the model has not been quantitatively validated with the Discharge Experiment, the model has shown interesting features of near-surface cusp confinement for the plasma.

Overview

The computational model employs a particle-in-cell Monte Carlo collision (PIC-MCC) method, treating high energy (primary) electrons, ions, and secondary electrons as particles. The PIC model is required due to a large range of collisionality for the densities and energies that are used in this effort and to resolve species interactions and dynamics in the cusp region. In the model, particles are moved in electric and magnetic fields while they are weighted to the computational grid at every time-step. The densities are then computed and used to calculate the electric potential and field on the grid. These calculations are repeated until the electric potential converges to a solution, as illustrated in Fig. 3.1. The key components of the model are:

• Particle tracking using a modified Boris method^1,2
• Particle weighting using generalized weighting scheme for axisymmetric domain³
• Anisotropic elastic scattering for electron-atom collisions⁴
• Magnetic field calculation using analytical equations for permanent magnets^5,6
• Adaptive Cartesian mesh generation
• Second order electric potential⁷ and field solver

Computational Methods

Particle Tracking and Collision Calculations

Each of the three charged species (primary electrons, secondary electrons, and ions) is modeled as macro-particles that are tracked separately, while neutrals are assumed uniform. The interaction between the three species is incorporated by the electrostatic forces acting on all the species. For the cases examined, the mean free path for electron collisions in the cusp region is long enough that the motion of the particles between elastic collisions can be described with an implicit particle-pushing algorithm. For this model, the Lorentz forces on the particles are decomposed into electric and magnetic forces using a modified version of a Boris-type particle pushing technique.^1,2,8 The electrostatic force is assumed to act half at the beginning, and half at the end of a given time step. This allows the circular motion of the particle due to the magnetic field to be treated in the absence of the electric field. The simple circular gyromotion of the particle is solved using a predictor-corrector algorithm as illustrated in Fig. 3.2. The predictor-corrector algorithm determines the particle trajectory for each time step by first predicting the midpoint of the trajectory assuming the magnetic field at the starting position. The magnetic field at the predicted midpoint is calculated and then used to determine a corrected final position and velocity orientation. This method has several advantages, when compared with other techniques, in both accuracy and speed.^1,2

For all species, only the collisions with background neutrals are implemented in the preliminary model since the neutral density is several orders of magnitude larger than the ion or electron densities in the experiment. The elastic collisions are approximated by the anisotropic scattering model given by Okhrimovskyy et al.⁴ for electrons and a variable hard sphere model⁹ for ions. The inelastic collisions (ionization and excitation) between the primary electron and background xenon atoms are approximated by simply depleting the primary electron current at every time-step, assuming that the primary electrons join the secondary electron population after a single inelastic collision event. These depleted currents as well as particles of different species are distributed to the grid using the first order area weighting. The particle densities and generation rate densities are then computed using corrected cell volume approach proposed by Verboncoeur.³ This approach removes the systematic errors associated with particle weighting in cylindrical coordinates.

Magnetic Field Calculation

In the computational model, analytical equations for magnetic fields induced by permanent block magnets⁵ and a cylindrical magnet⁶ are used to accurately calculate the magnetic fields. These equations assume a constant magnetization in a magnet, which is a reasonable assumption for the samarium cobalt magnets used in the experiment. In Ref. 6, an analytical equation for an axially magnetized ring magnet is provided, and from the equation, magnetic fields induced by a cylindrical magnet can be computed in the limit that the inner radius approaches zero. Both the magnetic field calculation and the particle tracking scheme determine the accuracy in computing important parameters such as the Larmor radius and the reflection point, which is critical in examining the behavior of different species at the cusp. The magnetic field configurations used in the simulations are shown in Fig. 3.3.

Adaptive Cartesian Mesh Generation

To provide adequate resolution in the near-cusp region and reasonable run times, an adaptive Cartesian mesh, with cell sizes comparable to the electron Larmor radius, is used. The mesh is generated by first creating a uniform Cartesian mesh throughout the domain and then subdividing the cells into four sub-cells depending on the magnetic field strength. For this effort, grid refinement was provided near the downstream cylindrical magnet. Fig. 3.4 shows the axisymmetric computational grid used to produce the model results discussed below.

Second Order Electric potential and field calculation

For a rectilinear mesh, the standard electric potential and field calculations using three point finite difference equations provide second order accurate solutions. However, for an adaptive Cartesian mesh, the accuracy drops by one order at the coarse/fine grid interface. In order to keep the second order accuracy, more sophisticated methods have to be implemented. The method implemented in the computational model approximates the local variation in electric potential with 2D quadratic spline equation.⁷ In computing the electric potential, the integral form of Gauss equation is applied to each cell by employing the finite volume method. Here, the electric flux across a cell face is approximated more accurately from the spline approximation to the electric potential using six cells neighboring the cell face.

In computing the electric field, the local variation in electric potential is again expressed with 2D quadratic equation given. For a given interior cell, the number of cells neighboring the current cell range from 7 to 13, which results in an over-determined linear system of equations when solving for the spline coefficients. Therefore, the least squares approach is used to obtain a unique solution for the coefficients. Once the coefficients for the quadratic equation are computed, then the electric field in the current cell can be calculated by taking the gradient of the quadratic equation. The electric field calculation method is tested for a case of a cylindrical domain with end plate biased at fixed potential. Figures 3.5 and 3.6 show the relative error in the electric field calculation for an adaptive grid using standard method and least squares quadratic fitting method, respectively. It can be readily seen that the standard method produces larger error, particularly at the grid interface, whereas the least squares quadratic fitting method produces significantly less error throughout the domain and most notable at the grid interface.

Preliminary Results

In the simulation of Discharge Experiment, the 25 eV primary electrons are injected into the discharge domain filled with 300 K xenon neutrals of uniform density 5×10¹⁹ m^-3, thus experiencing elastic and inelastic collisions with neutrals and producing ions and secondaries. The operating conditions for this experiment result in sparse plasma conditions with very low ionization levels.

Figure 3.10 shows the electric potential computed by the model. The highest potential value is seen close to the center of the domain while the potential drops rapidly toward the grounded boundaries. This potential structure pushes ions toward the walls. Figure 3.11 shows the contour plot of primary electrons. The density is relatively higher at the single magnetic cusp. By design, the stronger magnetic field strength at the ring-cusps reflect the electrons, while a larger number of electrons reach the downstream plate because of the relatively weaker field at that surface. The contours for ion and secondary generation rate density are similar to the primary density contour; thus, the vast majority of ions and secondaries are created very near the cusp for plasma condition created by the electron gun.

Figure 3.12 shows current density profiles at the downstream plate; these data are normalized to the maximum values for each species to show the proportionality of the current density. By taking full width at half maximum (FWHM) of these profiles, leak radii (summarized in Table 3.1) for individual species are obtained. In estimating leak radii from the experiment, the average FWHM values across the density peaks for radial lines from the center of the magnet are used. The computed primary loss radius is on the order of its Larmor radius, which is much smaller than the value obtained from the experiment. This disagreement is somewhat expected since, in the model, the electron gun is perfectly aligned with the cylindrical magnet, which results in primary electron velocities nominally directed largely along the axis. The primary electron tracking model has shown that the misalignment of the electron gun in the experiment reduces the electron population near the centerline field, resulting in relatively lower current density peaks shown in the data. As described above, ions and secondaries are mostly created close to the cusp. Because of the axial electrostatic force in this region, the ions are pushed toward the downstream wall as shown in Fig. 3.13, resulting in high ion current density near the cusp. This behavior results in ion loss radii much smaller than the ion gyroradius and comparable to the hybrid gyroradius evaluated for secondary electrons and ions. The secondary loss radius on the order of hybrid loss radius is explained by the magnetic field structure near the permanent magnet cusp. Unlike a spindle cusp, the permanent magnet creates a weakly convergent cusp field as shown in Fig. 3.3, and cross-section of secondary electron volume remains large even very near the cusp. By these observations, we see that the sparse plasma conditions used in this experiment do not exhibit the same mechanisms that are used to explain the hybrid gyro behavior for more highly ionized plasma.

Future Work

Once the analysis is extended to the partially ionized plasma, the model will employ a hybrid approach which will treat ions and slow electrons using a two-fluid approximation. The fluid equation will be computed numerically using finite volume method.

References

1. Wirz R., “Discharge Plasma Processes of Ring-Cusp Ion Thrusters,” Ph.D. Dissertation, Aeronautics Dept., Caltech, Pasadena, CA, 92005, 2005, http://resolver.caltech.edu/CaltechETD:etd-05232005-162628
2. Mao H.-S., Wirz R., “Comparison of Charged Particle Tracking Methods for Non-Uniform Magnetic Fields,” 42nd AIAA Plasmadynamics and Lasers Conference, Honolulu, HI, 2011, AIAA 2011-3739
3. Verboncoeur J., “Symmetric Spline Weighting for Charge and Current Density in Particle Simulation,” Journal of Computational Physics, Vol. 174, 2001, pp. 421–427
4. Okhrimovskyy A., Bogaerts A., and Gijbels R., “Electron Anisotropic Scattering in Gases: A Formula for Monte Carlo Simulations,” Physical Review E, Vol. 65, No. 037402, 2002
5. Engel-Herbert R. and Hesjedal T, “Calculation of Magnetic Stray Field of a Uniaxial Magnetic Domain,” Journal of Applied Physics, Vol. 97, No. 074504, 2005
6. Babic S. I., Akyel C., “Improvement in the Analytical Calculation of the Magnetic Field Produced by Permanent Magnet Rings,” Progress In Electromagnetics Research C, Vol. 5, 71-82, 2008
7. J. M., “Advances in Fully-Kinetic PIC Simulations of a Near Vacuum Hall Thruster and Other Plasma Systems,” Ph.D. Dissertation, Aeronautics and Astronautics Dept., Massachusetts Institute of Technology, Cambridge, MA, 02139, 2007
8. Boris J.P., “Relativistic plasma simulations-Optimization of a hybrid code,” Proceedings of the 4th Conference on the Numerical Simulation of Plasmas, NRL Washington, Washington DC, 1970
9. Bird G. A., “Molecular Gas Dynamics and the Direct Simulation of Gas Flows,” Clarendon Press, Oxford, 1994