Computational Modeling of a Lithium Magnetoplasmadynamic Thruster for Nuclear Electric Propulsion

Program: NASA Space Technology Graduate Research (NSTGRO25)
Principal Investigator: Prof. Mark Cappelli (Stanford University)
NASA Research Collaborator: Dr. Kurt Polzin (NASA Marshall Space Flight Center)
Role: Graduate Research Fellow

Project Overview

High-power electric propulsion, including Nuclear Electric Propulsion (NEP), is a promising technology for future missions to Mars. For high-power spacecraft, the self-field Magnetoplasmadynamic (MPD) configuration is especially attractive because it delivers strong performance.

A major barrier to operating MPDs at high current, which enables high specific impulse, is onset: a transition associated with large voltage fluctuations and severe anode erosion when the current exceeds a threshold. Onset is influenced by propellant choice, mass flow rate, and thruster geometry. However, its underlying causes remain insufficiently understood, and conventional computational tools such as MHD do not capture the physics that appears central to onset.

Goal: develop a predictive, physics-based simulation framework capable of predicting onset, to guide design and operating conditions that extend thruster lifetime for MPD thrusters relevant to Mars missions.

Technical contribution

By retaining non-equilibrium thermodynamics and ionization dynamics, the formulation captures electrothermal behavior that single-fluid MHD cannot.

Linear stability model

I implemented a two-fluid, 0D linear stability model for the near-anode region, where instabilities are most critical. In this framework, electrons and heavy species are treated with separate energy balances, enabling non-equilibrium between electron temperature and gas temperature. The model explicitly includes ionization and recombination kinetics, ambipolar diffusion, Ohmic heating with Spitzer conductivity, viscous dissipation, and electron–heavy species energy exchange.

The governing equations are nondimensionalized, yielding a cubic dispersion relation whose roots provide the local instability growth rate as a function of the normalized axial wavenumber.

One-dimensional two-fluid solver and coupling strategy

To evaluate stability from computed plasma states, I developed a 1D two-fluid solver that computes the axial evolution of plasma properties along the channel, including:

mass and momentum,
electron density,
electron temperature and heavy-species temperature,
self-induced magnetic field,
consistent chemistry terms: ionization–recombination, ambipolar diffusion, Ohmic heating, electron thermal conduction, and wall losses.

I then coupled the 1D solver to the 0D stability model: the solver provides axial profiles, a near-anode state is extracted, and the stability module computes the growth rate and evaluates an onset criterion for a magnetic field B that increases up to onset. For a given geometry, this workflow can be used to sweep operating conditions and record onset thresholds when the instability criterion is exceeded. This threshold is defined using the onset parameter, given by the square of the discharge current divided by the propellant mass flow rate, and its critical value marks the boundary between stable operation and the first observable onset symptoms. This coupled framework was exercised across representative geometries and operating conditions and was used to generate onset predictions that can be directly compared with published experimental trends.

My results

This coupled system was validated by reproducing key trends reported in the experimental literature:

Stability maps (Figure 1): Growth rate versus nondimensional wavenumber exhibits alternating stable and unstable bands.

Figure 1 — Growth rate vs. wavenumber at J = 2.0 × 10⁶ A/m²
Ionization-fraction threshold (Figure 2): Strong sensitivity to the near-anode ionization fraction. Above a critical value, modest current increases can trigger rapid stability loss.
$Growth rate vs. ionization fraction alpha$ Figure 2 — Growth rate vs. ionization fraction α
Interelectrode separation sensitivity: For a fixed mass flow rate, decreasing interelectrode separation increases stability (Figure 3). The onset parameter drops as interelectrode separation increases; the predicted inverse power-law captures the observed experimental trend (Figure 4).

Figure 3 — Growth rate vs. interelectrode separation

Figure 4 — Onset parameter vs. interelectrode separation (experiment vs. prediction)
Electrode-length sensitivity (Figure 5): The onset parameter increases with the length of the shortest electrode, as reported in the experimental literature.
Figure 5 — Onset parameter vs. length of the shortest electrode (experiment vs. prediction)
Mass-flow sensitivity (Figure 6): The onset parameter decreases with increasing mass flow rate; the model captures the experimentally observed dependence on mass flow rate.
Figure 6 — Onset parameter vs. mass flow rate (experiment vs. prediction)