A Particle-in-cell Simulation for the Traveling Wave Direct Energy Converter (TWDEC) for Fusion Propulsion Andrew Chap* University of Maryland, College Park, MD 20740 USA Alfonso Tarditi† Electric Power Research Institute, Inc., Knoxville, TN 37932 USA John H. Scott‡ NASA Johnson Space Center EP3, Houston, TX 77058 USA A Particle-in-cell simulation model has been developed to study the physics of the Traveling Wave Direct Energy Converter (TWDEC) applied to the conversion of charged fusion products into electricity. In this model the availability of a beam of collimated fusion products is assumed; the simulation is focused on the conversion of the beam kinetic energy into alternating current (AC) electric power. The model is electrostatic, as the electrodynamics of the relatively slow ions can be treated in the quasistatic approximation. A two-dimensional, axisymmetric (radial-axial coordinates) geometry is considered. Ion beam particles are injected on one end and travel along the axis through ring-shaped electrodes with externally applied time-varying voltages, thus modulating the beam by forming a sinusoidal pattern in the beam density. Further downstream, the modulated beam passes through another set of ring electrodes, now electrically floating. The modulated beam induces a time alternating potential difference between adjacent electrodes. Power can be drawn from the electrodes by connecting a resistive load. As energy is dissipated in the load, a corresponding drop in beam energy is measured. The simulation encapsulates the TWDEC process by reproducing the time-dependent transfer of energy and the particle deceleration due to the electric field phase time variations. Nomenclature a Acceleration, m/s² C Capacitance, F E Electric field vector, V/m m Mass, kg P Power, W q Ion bunch charge, C Q Electrode charge, C r Radial position, m R Resistance, Ω v Velocity, m/s v Velocity vector, m/s W Weight x Position vector, m z Axial position, m α Specific mass, kg/kW Δh Spacing between grid nodes, m Δt Time-step value, s ρ Charge density, C/m³ Φ Electric potential, V Subscripts and superscripts i Radial grid node j Axial grid node n Time-step Physical Constants c Speed of light, 3 × 10⁸ m/s ε₀ Permittivity of free space, 8.85 × 10⁻¹² F/m *Graduate Research Assistant, Department of Aerospace Engineering, 3710 Martin Hall, AIAA Student Member. †Project Manager, Electric Power Research Institute, Inc., 942 Corridor Park Blvd, Knoxville, TN 37932 ‡Chief, Energy Conversion Branch, NASA Lyndon B. Johnson Space Center/EP3, Houston, TX 77058 I. Introduction To change a crewed mission to Mars from a 3-year epic event to an annual (or sub-annual) expedition, new power and propulsion systems must be developed beyond the capabilities of both chemical rockets as well as nuclear thermal rocket designs. The key parameter to characterize the mission capability of a power and propulsion system is the specific mass, α, defined as the ratio of total mass to power output of the propulsion system (typically expressed in kg/kW). Under current Mars mission architectures (including nuclear thermal propulsion, solar electric propulsion, and solid core fission nuclear electric propulsion) a conjunction class mission to Mars has a manned mission duration of two to three years.¹ Fusion concepts as well as advanced fission concepts have been proposed as the power system for electric propulsion systems.² Those concepts that release most of their energy in the form of neutrons require heavy shielding to protect the crew and systems from radiation, as well as a heat engine to convert the nuclear energy into electricity, necessitating large radiators to dispose of waste heat (figure 1 left). Aneutronic fusion and some types of fission, which release the majority of their energy in the form of charged particles, may utilize either heat energy conversion or the direct conversion of charged particles into electricity (such as in the TWDEC). Traditional heat energy conversion requires heavy radiators as well as additional power management and distribution hardware to condition the DC power output into AC suitable for the RF antenna in an ion engine. The TWDEC produces little waste heat, and can be physically tuned to match the frequency required for the ion engine RF antenna. By relaxing the mass requirements of radiation shielding, waste heat radiators, and DC to AC power management, a system using aneutronic fusion and a TWDEC may achieve a specific mass as low as 3 kg/kW (figure 1 center). The direct use of the momentum of fusion products for propulsion can further reduce specific mass. In this case the fusion products pass through the TWDEC en route to the direct plasma thruster. The TWDEC is used to decelerate the ions by only a small amount, to generate electricity for the spacecraft systems and maintain fusion core operation. The majority of the ion momentum is carried on to be used directly as thrust (figure 1 right). [Figure 1. A simplified comparison of the required component masses of three power and propulsion systems for a given jet power output. Box sizes roughly correspond to mass.] II. The Physical TWDEC The principal mechanism behind the TWDEC is analagous to a linear particle accelerator working in reverse. Rather than putting energy into accelerating a beam, the TWDEC decelerates a beam and energy is extracted. The beam path may be considered in three sequential regions: the modulator, an intermediate region, and the decelerator, as displayed in figure 2 and outlined in the following sections. A. The modulator For the beam to become suitable for energy extraction it must be modulated from a homogeneous state into a modulated or "bunched" state. In the modulator section, a series of ring electrodes with time varying voltages create a wave of electric potential that travels with the beam, accelerating ions from the crests [Figure 2. Design for the NASA Johnson Space Center TWDEC test article. The ion beam travels left to right, through the modulator (A), intermediate region (B), and decelerator (C).] towards the troughs, resulting in a roughly sinusoidal modulation of beam velocity while leaving unchanged the aggregate beam momentum. B. The intermediate region In the intermediate region, the perturbations in velocity space born in the modulator are given time to mature into density variations in the axial dimension, so that upon entering the decelerator section the ions are effectively "bunched" with a sinusoidal-like axial variation in density. C. The decelerator In the decelerator section the bunched ion beam passes through electrically floating electrode rings. The voltage of a floating electrode is determined by its surrounding environment; the proximity of a positively charged ion bunch will raise the floating electrode's voltage. A traveling wave of sinusoidal beam density variation, then, will induce a sinusoidal voltage variation in a floating electrode. By adjusting the spacing of the electrodes, alternating voltage differences between electrodes arise. Power is produced by connecting these electrodes in a resistive circuit. A detailed explanation of bunch deceleration As the ion bunches induce time-varying voltages on electrodes, the ion bunch dynamics are in turn affected by the presence of these electrodes and any resistance in a circuit that might connect them. If passing through a single electrode, an ion bunch will experience an attraction to the electrode, due to the rearrangement of surface charge on the electrode that minimizes field energy. The ion bunch will accelerate axially towards the electrode as it approaches, and will decelerate accordingly after passing through, resulting in no net change in velocity and energy.ᵃ Similarly, if two electrodes are connected without resistance, no net deceleration occurs, as the charge transfer is immediate between electrodes, similar to an elastic collision. The two electrodes here act as one equipotential conductor (figure 3, middle). However, if a resistance is introduced between the electrodes, the transfer of charge is delayed, resulting in a net deceleration of the particle bunch. For the following steps, refer to the third plot in figure 3. When approaching the first electrode, the bunch is accelerated, but to a lesser degree than before, as the resistor delays the negative charge flow into that electrode (1). After passing through the first electrode, the bunch is decelerated (2). Upon reaching the midpoint of the electrode pair, electrode 1 remains positively charged while electrode 2 is negatively charged, due to the non-instantaneous nature of their discharging. At this point the deceleration continues, almost until the pair discharges (3), and the bunch undergoes a short acceleration in approaching the second electrode, reaching a local maximum in velocity (4). Upon passing through the second electrode the bunch undergoes deceleration due to an attractive force to the negatively charged electrode 2 (5). ᵃThis simple explanation neglects the expansion of the beam due to its space charge. In practice, radial beam expansion can be limited by an axial magnetic field. [Figure 3. Top: The axisymmetric 2-D domain with an ion bunch (here modeled as a semicircle which extrudes a sphere) passing through the centers of two ring electrodes (squares) Middle: A plot of bunch velocity vs. axial distance when the electrodes are connected by a wire with no resistance Bottom: The same plot with the electrodes connected by a 500 Ω resistor, showing a net deceleration of the ion bunch. The relatively low magnitude of velocity loss is due to the pre-optimization stage of the model, as well as the lack of magnetic field to reduce bunch expansion.] III. The Numerical Model The physics of the TWDEC is modeled with a particle-in-cell method that couples the Newtonian equations of particle motion with the discrete Poisson equation to solve for the electric potential on cell nodes. The simulation exists in a 2-D axisymmetric radial-axial (r-z) domain, and it is assumed that motion and variance in the azimuthal dimension are negligible (d/dθ = 0). A. The particles The particles (ions) to be simulated are amalgamated into macroparticles to reduce computation times. The weight of each macroparticle, i.e. the number of particles each macroparticle represents, is proportional to the inverse of the macroparticle's initial distance from the axis upon generation. The macroparticle weighting is only dependent on this initial radial distance; the weight of each macroparticle will remain constant over its lifetime in the simulation. The purpose of this non-uniform weighting is to have a uniform density of macroparticles in the 2-D plane. B. Particle motion Particle motion is governed by the electrostatic equations of motion: dx/dt = v dv/dt = (q/m)E. Discretizing the time variable into discrete time steps of value Δt, first order finite differencing is applied to the derivatives, with Δt kept small enough (in accordance with the CFL condition) to assume other variables stay constant over this time increment. (xⁿ⁺¹ - xⁿ)/Δt = vⁿ⁺¹/² (vⁿ⁺¹/² - vⁿ⁻¹/²)/Δt = (q/m)Eⁿ Solving for the unknowns reveals the leapfrog method for particle transport. xⁿ⁺¹ = xⁿ + vⁿ⁺¹/²Δt vⁿ⁺¹/² = vⁿ⁻¹/² + (q/m)EⁿΔt C. The matrix equation At each time step, given the charge distribution on the cell nodes, a system of equations is set up to simultaneously solve for the potential on the cell nodes, the surface charge distribution on the electrodes, and the current flow between electrodes. The borders of the floating electrodes exist on the cell nodes. 1. Solving for the electric potential The potential is determined by the charge density in Poisson's equation for electrostatics, ∇²Φ = -ρ/ε₀. In cylindrical coordinates with azimuthal symmetry (d/dθ = 0), Poisson's equation is ∂²Φ/∂r² + (1/r)(∂Φ/∂r) + ∂²Φ/∂z² = -ρ/ε₀. A first order finite difference approximation results in (Φ_{i+1,j} - 2Φ_{i,j} + Φ_{i-1,j})/(Δh)² + (1/(iΔh))((Φ_{i+1,j} - Φ_{i-1,j})/(2Δh)) + (Φ_{i,j+1} - 2Φ_{i,j} + Φ_{i,j-1})/(Δh)² = -ρ_{i,j}/ε₀ simplifying to -4Φ_{i,j} + Φ_{i+1,j}(1 + 1/(2i)) + Φ_{i-1,j}(1 - 1/(2i)) + Φ_{i,j+1} + Φ_{i,j-1} = -(ρ_{i,j}/ε₀)(Δh)². For the special case when i = 0 (r = 0) it has been shown³ that -6Φ_{0,j} + 4Φ_{1,j} + Φ_{0,j+1} + Φ_{0,j-1} = -(ρ_{0,j}/ε₀)(Δh)². Thus for every discrete unknown Φ in the domain there is a constraining equation. These equations can be put into matrix form and solved implicitly. 2. Solving for the surface charge distribution on the decelerator electrodes The electrodes are approximated as perfect conductors, therefore the solver imposes the same potential on all nodes which reside on the electrode. Internal nodes are trivial and are neglected from the solver. The charge density on the surface node of the electrode is unknown and must be solved for, so these are added to the vector of unknowns. The conservation of charge on the surface provides an additional constraining equation, closing the system. 3. The circuit solution For electrodes connected in a circuit, two additional unknowns, the total charge on each of the two electrodes, are introduced. Two additional constraining equations, then, are needed to close the system. As the potential difference builds between electrodes, current flows as in a discharging capacitor: dQ/dt = (ΔV/R) e^(-t/(RC)). Where C is the capacitance between the electrodes. This can be trivially approximated using finite differencing. Qⁿ = Qⁿ⁻¹ + (ΔV/R) e^(-Δt/(RC)) Δt The second constraining equation is the conservation of charge between the two electrodes: Q₁ + Q₂ = 0 [Figure 4. The matrix system for the model, solving simultaneously for the electric potential at all nodes, the surface charge, and the total charge on each electrode] The system can now be solved as an implicit matrix system. Solving is accomplished in MATLAB via the black box "backslash" operator. D. Particle-cell interpolation At each time step, values must be commuted from the particle locations to the cell nodes and vice versa. Quadratic interpolation⁴ has been found to effectively avert the manifestation of non-physical furcation in beam expansion when lower order interpolation functions are used. The weighting scheme for quadratic interpolation is as follows. W(x) = 3/4 - (x/Δh)², |x| ≤ Δh/2 1/2 (3/2 - x/Δh)², |x| ≤ Δh/2 0, otherwise The weighting schemes for each interpolation function are plotted in figure 5. A visual comparison of the effect of interpolation functions on an expanding beam can be seen in figure 6. [Figure 5. A graphical comparison of the different interpolation weighting methods.] [Figure 6. A qualitative comparison of the effects of different interpolation methods on an expanding beam. A 100 × 50 grid is used, with approximately 33,000 particles. The chosen interpolation method in each case is used for both particle-to-grid as well as grid-to-particle interpolation.] IV. The Quasistatic Approximation The justification for a quasistatic treatment of ion dynamics is dependent upon proof that the electromagnetic energy radiation from a decelerating ion or ion bunch is negligible in comparison to the magnitude of direct energy conversion. For a test condition, an ion bunch losing 10% of its velocity through a pair of decelerator rings is considered.ᵇ Given a p-B¹¹ α-particle fusion product with a velocity of 10⁷ m/s and a spacing between electrodes of 5 cm, an ion will experience a deceleration on the order of 10¹¹ m/s². a = (1/2) ((Δv)² / Δx) = (1/2) ((10⁶ m/s)² / 0.05 m) = 10¹¹ m/s² The power of this wave (assuming an ion bunch charge of 10⁻⁹ C) can then be calculated using the Larmor formula P = (1/(6πε₀)) (q²a²/c³) = (1/(6πε₀)) ((10⁻⁹ C)² (10¹¹ m/s²)² / c³) ≈ 10⁻¹¹ W The deceleration will take place over a time of approximately 10⁻⁸ s, resulting in a radiated energy on the order of 10⁻¹⁹ J. Simulation results have shown that a similarly charged bunch will lose on the order of 10⁻⁸ J when undergoing deceleration between two electrodes. The radiated energy of a decelerating bunch is many magnitudes lower than the electrostatic energy transfer, therefore the quasistatic approximation is valid for the model. V. Computational Acceleration with GPU Parallelization Computation times are reduced through partial execution on a high performance graphics processing unit (GPU). Parallelization is achieved through Jacket in MATLAB. Jacket provides black box operation for running iteration-independent loops and vectorized operations in parallel on the GPU, resulting in significant speed increases over high-end CPU processing. While some operations in the PIC code benefit prodigiously from GPU parallelization, other operations are unable to be executed on the GPU under the current capability of the software. In particular, the large matrix inversion operation involved in solving the Poisson equation, as well as a portion of the charge distribution interpolation, are performed on the CPU. A breakdown of these computation times is provided in figure 7. ᵇIt is not expected for a single pair of decelerator rings to ever achieve this magnitude of deceleration [Figure 7. A comparison of the amount of time required to execute a single time step in a simulation on a 1000 × 375 grid with 1,000,000 particles. The CPU is a 12 core (2×6) 2.8 GHz Xeon X5660 and the GPU is an NVIDIA Tesla C2070.] VI. Model Capabilities The PIC model conceived of and developed by the authors is currently capable of modeling the modulation and deceleration of an ion beam. The model has been partially validated by noting the conservation of energy between the beam kinetic energy, the electric field energy, and the resistive load energy. Further validation was provided by comparison of various model sub-physics with replications in COMSOL. The model is suitable for investigating the relation between energy conversion efficiency and electrode shape and positioning. Although past computational⁵ and experimental⁶ studies on the TWDEC have investigated the achievable beam deceleration efficiency, no studies have investigated computationally or experimentally the energy conversion mechanism at play in these devices. [Figure 8. Numerical simulation of the TWDEC with heuristic parameters. Dots denote macroparticle locations, boxes denote electrodes, and lines of constant potential are plotted in color. The beam travels from left to right. Four modulator rings apply a traveling wave in potential with a wave velocity equal to the beam velocity. The downstream decelerator ring pairs experience alternating potentials.] VII. Future Work Future work will include improvements to the model design, as well as optimization routines to direct the construction of an experimental TWDEC test article at NASA's Johnson Space Center. Model improvement To improve the physics of the model, simulation of the collision of ions with background gas will be introduced, which will allow to model remain accurate for a laboratory vacuum. Because an integral part of TWDEC operation is the magnetic field to mitigate beam expansion, the motion of the gyrocenters of macroparticles in a magnetic field will be modeled. A proposed method of neutralization of the ion beam center with an electron beam will also be considered. TWDEC optimization The primary motivation for this work is the optimization of the TWDEC specific mass. Of particular interest is the way in which electrode design and spacing affects energy transfer efficiency, as well as what efficiencies are achievable for a given beam energy and density. The model will also be used for future vehicle studies for which the TWDEC may be considered. VIII. Acknowledgments This work was supported by the University of Maryland and by NASA Lyndon B. Johnson Space Center. References ¹B. Drake, Human Exploration of Mars Design Reference Architecture 5.0 (2009) ²J. Scott, J. George, A. Tarditi, Direct Conversion for Low Specific Mass In-Space Power and Propulsion, Proceedings of Nuclear and Emerging Technologies for Space (2013) ³E. Chao, S. Paul, R. Davidson, K. Fine, Direct Numerical Solution of Poisson's Equation in Cylindrical (r, z) Coordinates (1997) ⁴R. Hockney, J. Eastwood, Computer Simulation Using Particles, Taylor & Francis Group, New York, NY (1988) ⁵R. Kawana, M. Ishikawa, H. Takeno, T. Yamamoto, Y. Yasaka, Performance analysis of small-scale experimental facility of TWDEC, Energy Conversion and Management Vol. 49 Iss. 10 (2008) ⁶H. Shoyama, H. Nakashima, Y. Kanda, Two-Dimensional Analysis of Energy Conversion Efficiency for a Traveling Wave Direct Energy Convertor (1995)