Whole device modeling of the fuze sheared-flow-stabilized Z pinch

The FuZE sheared-flow-stabilized Z pinch at Zap Energy is simulated using whole-device modeling employing an axisymmetric resistive magnetohydrodynamic formulation implemented within the discontinuous Galerkin WARPXM framework. Simulations show formation of Z pinches with densities of approximately 1022 m−3 and total DD fusion neutron rate of 107 per µs for approximately 2 µs. Simulation-derived synthetic diagnostics show peak currents and voltages within 10% and total yield within approximately 30% of experiment for similar plasma mass. The simulations provide insight into the plasma dynamics in the experiment and enable a predictive capability for exploring design changes on devices built at Zap Energy.


INTRODUCTION
Whole device modeling of the FuZE shear-flow-stabilized (SFS) Z pinch [1,2] at Zap Energy is illuminating key physics of the device and guiding progress towards a breakeven fusion concept with a high-gain powerproducing power plant.Two-dimensional axisymmetric resistive magnetohydrodynamic (MHD) simulations using the discontinuous Galerkin finite element WARPXM framework [3] are compared to experimental measurements from the FuZE device using a suite of synthetic diagnostics.The modeling unveils the plasma dynamics inside FuZE and enables in-depth investigation of various areas including plasma compression adiabaticity and electrode erosion.Validated modeling provides a predictive tool for next-step SFS Z pinch experimental design.This paper summarizes simulation results using synthetic diagnostics to compare with experimental analogues.Section 2 describes the MHD model used while Sec. 3 briefly describes the numerical method.Section 4 describes the problem setup, Sec. 5 discusses results, and Sec.6 gives conclusions.

MODEL
For FuZE plasmas, expected temperatures are on the order of 1 keV for a pinch size of about 1 cm with enclosed currents of up to 500 kA [1,4].In this regime, the ion Larmor radius and ion inertial length are on the millimeter scale.For the MHD model to be valid, these length scales should be small compared to the system [5,6].For a 1 keV pinch with 500 kA current, a radius of 1 cm would be on the order of 10 ion Larmor radii, so we proceed with a resistive MHD model given by ∂ρ ∂t Equations ( 1)-( 9) are written in non-dimensional form with normalized cyclotron frequency (ω c τ ) and collision frequency (ν p τ ) defined as described in Ref. [7].The equations are solved, as described in the next section, in two-dimensional (r-z plane) axisymmetric form where ∂ ∂θ = 0.
A blackbody radiation model [8] is added in Eq. ( 3).Viscosity and thermal conduction terms are added as well, as seen in Eqs. ( 2), ( 3), (8), and (9).These terms are required to match the simulation dynamics with experiment, as described further in Sec. 5.

NUMERICAL METHOD
The MHD model described in Eqs. ( 1)-( 9) can be combined and expressed compactly as The equations are solved using the discontinuous Galerkin method within the WARPXM framework developed at the University of Washington and Zap Energy [3].An overview of the method and implementation is found in Refs.[3], [5], and [7].The discontinuous Galerkin method is well suited for solving the predominantly hyperbolic MHD equations on unstructured meshes using simplex finite elements, allowing for simulation of the FuZE geometry.In this work, a cylindrical formulation of the MHD equations described in Sec. 2 is used, the numerical handling of which is described in Ref. [5].Artificial dissipation is also added to stabilize solutions [9].
A circuit model is implemented to capture the interaction between the driving current bank and the plasma in FuZE in a self-consistent manner.The circuit is modeled by including the plasma dynamics as a time-dependent impedance and associated load voltage as a circuit element.The current is driven by the circuit with prescribed capacitance, inductance, resistance, charge voltage, and trigger time.

PROBLEM SETUP
The simulations performed aim to match conditions for a set of experimental operating parameters on FuZE as shown in Fig. 1.Triangular finite elements are used with a quadratic solution basis.A slug of plasma of 0.5 mg in the acceleration region is initialized where the valves inject gas.A volumetric injection of plasma is also input at a rate of 3 mg / ms at the valve location, mimicking breakdown of residual neutral gas and ongoing valve injection.
Simulated initial mass is chosen by following the trajectory of the magnetic field data for an experimental plasma pulse (#230518072) and calculating the corresponding plasma mass.Valve mass injection rates are also estimated, as described in Sec.4.1.A series RLC circuit is used to simulate the driving capacitor bank for FuZE, with circuit values C = 222 µF, L = 286 nH, and R = 1.5 mΩ.At t = 0, the capacitor charged to 25 kV is discharged.The plasma load is calculated as a voltage drop measured along the insulator of the device between the cathode (inner electrode) and anode (outer electrode).The gap voltage is measured upstream of the insulator ahead of a 70 nH region of inductance (see Fig. 1 in Ref. [10]), so to compare with experiment, it is calculated as The circuit current, converted to an azimuthal magnetic field, imposes a boundary condition on the plasma at the insulator, where the axial velocity is allowed to float if into the domain and set to zero otherwise.At r = 0, axisymmetric boundary conditions are employed, as described in Ref. [5].All other boundaries are treated as free-slip conducting walls that are thermally insulating by imposing zero-normal density and pressure gradients.
FIG. 1. Two-dimensional r-z plane slice of FuZE used as the simulation domain.Gas is puffed into the coaxial acceleration region where it is ionized by an electrical power supply between the cathode (inner electrode) and anode (outer electrode) and pushed into the assembly region to form a pinch.The simulation assumes the puffed gas is fully ionized as an initial plasma slug at t = 0. Note that the radial dimension is stretched roughly by 3-to-1.

Mass determination from slug and valve injection models
The mass of the plasma in experiment is estimated using a combination of magnetic probes along the outer electrode and gas valve injection formulas.Using the magnetic probes, the motion of an assumed slug of plasma injected by the gas puff valves is tracked to estimate its mass.Using this slug model, the plasma mass is determined using the force [11,12], Dividing Eq. ( 11) by m and integrating twice to time, t d , leads to an equation for distance traveled by the plasma, In the experiment, the times taken for the plasma slug to move from the gas puff valves to the magnetic probes along the outer electrode are known by measuring when each probe records a jump in magnetic field.A space-time plot of this magnetic probe data as shown in Fig. 2 indicates where and when along the outer electrode the jumps occur.With d and t d known for various magnetic probe locations, m is calculated using Eq. ( 12).The dashed orange line in Fig. 2 corresponds to m = 0.5 mg for the particular plasma pulse.
A combination of gas valve injection and a deflagration process ionizing remaining neutral gas sustains plasma flow to the pinch [11].To account for that effect, a mass injection model is needed.A simple model assumes gas flows out of the valve at the speed of sound such that Normalizing by non-dimensional parameters, ρ = ρ/ρ 0 and t = t/τ v where τ v = V /c s0 A is the timescale of the venting valve and assuming adiabatic expansion such that p = p/p 0 = (ρ/ρ 0 ) γ = ργ leads to d ρ The solution for ρ (0 1−γ , which describes the normalized density in the valve.The normalized density entering the device is thus ρ = 1 − γ−1 2 t + 1 2 1−γ .Substitution of the dimensional values leads to an expression in terms of initial density, ρ 0 , which is recast in terms of valve pressure, p, using the ideal gas law for deuterium (assuming V = 1 cm 3 and T = 300 K).A reference mass calculated for a pressure of 100 PSIA is M 0,ref = ρ 0,ref V = 1.11 mg.For different pressures, this reference mass is scaled linearly, leading to mass as a function valve pressure and opening time as where ∆ is a timing delay between the valve trigger and opening, N v is the number of valves used, and p is measured in PSIA.For the experimental plasma pulse in Fig. 2, N v = 4, p = 149.1 PSIA, and t = 0.26 ms.For deuterium (γ = 7/5, T = 300 K), the sound speed is 934 m/s leading to τ v = 1.070 ms (assuming A = 1 mm 3 ).Equation ( 14) is solved for ∆ using m = 0.5 mg, yielding 0.175 ms.Equation ( 14) is then differentiated with respect to time, yielding ṁ = 3 mg / ms.Such an injection rate may be considered a low estimate, since deflagration is expected to make a significant contribution as an upstream plasma source.Modeling with neutrals to capture deflagration is the topic of ongoing research.
FIG. 2. Enclosed current in experiment as measured by magnetic probes along the outer electrode as a function of time.A few probes in the assembly region did not register data for this pulse.
A simulation (#s230523033) is thus run to 30 µs matching experimental pulse settings with 0.5 mg slug mass and 3 mg / ms injection rate.The initial slug of plasma sits in the acceleration region where the valves are located, as shown in Fig. 1.

RESULTS
Traces of current and voltage are shown in Figs.   Figure 5 shows the evolution of plasma density along with lines of enclosed current.The plasma slug is pushed from the acceleration region to the assembly region through the J × B force imposed by the circuit current.A true slug motion would be indicated by the lines of radial current through the plasma as it moves through the acceleration region.At t = 8 µs, deviation from the 1D slug model is apparent, as axial components of the radial contours are seen near the outer electrode.Such behavior represents a blowby phenomenon where the imposed current does not fully couple to the plasma [13].The viscous (ν = µ/ρ) and heat (κ = kA i /ρ) diffusivities in this simulation are set to 200 m 2 / s which lead to the good agreement with experimental values shown in Figs. 3 and  4. Simulations with increased viscous and heat diffusivities show reduced blowby and increased peak current and voltage.After the slug reaches the assembly region and the flux rebounds after hitting the endwall, a pinch then forms on axis, as seen in the high density on-axis regions in the 14.7 µs and 16 µs plots.
Figure 6 shows the volumetric neutron yield rate calculated using Bosch-Hale thermonuclear fusion formulas [14] at the time of maximum total yield as well as plasma temperature.The thermonuclear yield is a function of the density and temperature.The combination of high temperature on axis above 3 keV combined with on-axis densities above 10 22 m -3 leads to the on-axis yield rate, indicating the existence of a Z pinch in the simulation.
FIG. 5. Plasma density at 8 µs, 12 µs, 14.7 µs (maximum total domain-integrated yield rate time), and 16 µs, overlaid with lines of enclosed current.The plots show current sheet canting and blowby as the magnetic flux pushes the plasma slug into the assembly region.Plasma thermal energy and magnetic flux accumulate at the endwall, causing a rebound behavior that completes pinch formation.The 16-µs plot shows the pinch formed on axis with noticeable m = 0 activity associated with reduced sheared flow.
Figure 7 shows the total integrated yield rate in the simulation versus experimental data.The experimental data is of fast plastic scintillator measurements scaled to instantaneous neutron yield rates using total yield measurements by Lanthanum Bromide detectors.The simulation data is an analogous synthetic diagnostic calculating the fulldomain integral of the simulated volumetric yield rate (the top plot in Fig. 6).The dashed lines are integrations of the total yield rates over time.The difference in simulation and experimental yield indicates an overprediction in simulation by about 300%.An adjustment to the simulation yield rate calculation is made by removing yield where the difference between velocity at the pinch edge, the radius a where p(a) = p max /4 [5], and the axis is below 10% of the Alfvén speed at a.The correction removes yield for low sheared flows assumed to be m = 1 FIG. 6. Neutron yield rate calculated from Bosch-Hale formulas for thermonuclear fusion [14] (top) and plasma temperature (bottom) at 14.7 µs, the time of maximum total domain-integrated yield rate. Lines of enclosed current are also shown.The plasma reaches above 3 keV near the nosecone on axis in the assembly region.This high temperature combined with collimated density on axis lead to the region of thermonuclear fusion seen in the neutron yield rate plot.
unstable [15,5].The total yield calculated from the correction is about 70% of the experimental yield.
In simulations without viscous and heat diffusivities (not shown), the yield rate increases early compared to the experiment.Increased blowby in such cases lead to earlier flux rebound from the endwall and thus earlier neutron production.This again suggests some viscous dissipation is needed to match experimental results.Additionally, simulations with no radiation have higher neutron production at times after the peak, unlike in the experiment where neutron production drops significantly after about 16 µs.The addition of the blackbody radiation term reduces this late time yield in the simulation, similar to experiment.The physical value for the blackbody radiation constant used in simulations is C bb = 4.559 × 10 52 J −3 s −1 .
FIG. 7. Total yield comparison between experiment and simulation.The experimental yield rate (solid blue) is found using plastic scintillators calibrated such that the time-integrated yield agrees with Lanthanum Bromide detectors.Total volume yield rate from simulation is calculated by integrating the volumetric yield as shown in Fig. 6 (solid orange).The solid green line also shows the simulation volume yield rate after removing yield in regions a pinch would be assumed unstable.Dashed lines integrate the yield rate to determine the total yield.
Further work to simulate the plasma dynamics with more realistic viscosity and radiation models are ongoing.This includes the consideration of Braginskii [16] terms as well as bremsstrahlung and line-radiation models [17].
Capturing the true nature of the endwall flux rebound is also being considered.Current modeling does not include the possibility of outflow through the endwall which has a spoked transparent geometry [18] or dissipation due to ion recycling and re-ioniziation [19] that may weaken the flux rebound.However, it is unlikely that there is full outflow [18] so, in future work, boundary conditions and related physics (including 3D effects) will be developed to capture realistic endwall behavior.
Figures 8 and 9 show simulation synthetic analogues of the space-time enclosed current plot shown in Fig. 2. Figure 8 shows the enclosed current at the outer electrode with the 0.5 mg slug trajectory superimposed.The lack of agreement with Fig. 2 is understood by observing Fig. 5. Between 12 µs and 16 µs, a clear detachment of enclosed current is seen along the outer electrode due to the formation of eddy currents, which reduce the total enclosed current measured.In contrast, the experiment shows strong magnetic fields at the outer electrode in the assembly region as seen in Fig. 2. For the simulation, a similar plot is constructed at a smaller radius inside of the eddy currents.Choosing 7 cm as shown in Fig. 9 indicates a magnetic space-time structure more in line with the experiment in Fig. 2, with the 0.5 mg trajectory close to the magnetic field wave front.

CONCLUSION
Comparisons of WARPXM simulations using two-dimensional axisymmetric resistive MHD with the FuZE experiment at Zap Energy show a remarkable level of agreement.Comparisons using a suite of synthetic diagnostics illuminate these agreements while suggesting areas that modeling does not fully capture.However, these results show that the 2D axisymmetric MHD model can be used as a predictive tool for exploring alternative geometries, electrode configurations, driving circuit parameters, and gas fill profiles which can accelerate progress toward breakeven fusion gain.Improvements to modeling are possible and will improve prediction accuracy.Extension of the MHD model to 3D can capture non-axisymmetric effects leading to m = 1 instability.The inclusion of multi-fluid physics [20,3,5] can help to accurately capture blowby, sheath effects, and instability growth rates.Multi-species continuum kinetic modeling can allow for further accuracy of sheath modeling and pinch dynamics at small pinch radius where the magnetic field goes to zero [21], with hybrid approaches helping to speed up such simulations [22,7].Other advancements include the addition of a plasma-neutral model [23] to capture plasma deflagration [11] and more accurate boundary conditions that can capture heat and particle fluxes to walls, and corresponding erosion and impurity generation [17].These advancements are actively being developed to help capture the plasma dynamics of FuZE and other devices being built at Zap Energy [4].

3 and 4
with comparisons to experimental discharge #230518072.In all experimental plots, timings are shifted by 4.2 µs to when the current rise starts.

Figure 3 FIG. 3 .
Figure3compares experimental capacitor bank current output with that of the simulation coupling the equivalent RLC series circuit described in Sec. 4 with the MHD plasma.The currents reach similar peaks before dropping, with slightly lower currents later in the experimental pulse compared to simulation.

Figure 4
Figure4shows a comparison of V gap between experiment and simulation.Early in time, the simulation exhibits jumps in voltage not present in experiment, an artifact of the fully ionized MHD assumptions of the initial slug.The incoming plasma flux is reflected by the slug, causing the jumps in V gap .Nonetheless, the gap voltage values in experiment and simulation are similar and have comparable temporal profiles with the drop in voltage at around 15 µs.The voltage drop is associated with pileup of plasma and flux at the endwall that drives a wave of magnetic flux back into the acceleration region where it interacts with the circuit.

FIG. 8 .FIG. 9 .
FIG. 8. Enclosed current in simulation at outer electrode radius.Large eddy currents in simulation reduce these values in the assembly region.FIG. 9. Enclosed current in simulation at r = 7 cm.At this radius, some of the eddy currents are removed.The 0.5 mg trajectory line matches better than in Fig. 8.