Validation of extended magnetohydrodynamic simulations of the HIT-SI3 experiment using the NIMROD code

A spheromak¹ is a simply connected magnetized plasma configuration that is typically formed and sustained through the injection of magnetic helicity. After being injected, the magnetic structure experiences a relaxation event where an equilibrium that has minimized free energy forms.^2,3 While the method of injecting magnetic helicity varies, the standard method is to balance the resistive decay of helicity in the spheromak with injection. Examples of these helicity injection schemes include inductive drive in a tokamak, coaxial helicity injection (CHI), which has been used to form both spheromaks^4,5 and spherical tokamaks,^6,7 localized helicity injection (LHI) used on the Pegasus device,⁸ oscillating-field current drive,^9,10 in reversed field pinches (RFP), and steady inductive helicity injection (SIHI) on the Helicity Injected Torus with a Steady Inductive helicity injection (HIT-SI)¹¹ and HIT-SI3 devices.¹² 

The HIT-SI3 experiment forms and sustains spheromaks through the use of SIHI.¹¹ These injectors consist of half-toroids with two sets of coils to inductively drive an EMF along magnetic flux producing helical fields oscillating at a given frequency (⁠$finj=5−70$ kHz). Since both sets of coils are oscillating in phase, the handedness of the injected helicity does not change. HIT-SI3 uses the bow-tie flux conserver geometry originally used by the HIT-SI device, but differs in the number and location of the helicity injectors. While the HIT-SI device consisted of two injectors, one each on the top and bottom of the device, HIT-SI3 consists of three, all located on the “top” of the machine. These geometries can be seen in Figs. 1 and 2.

Operation of the device consists of three stages. The first stage is the build-up of a non-axisymmetric injector-driven state and concludes with a large scale relaxation event over an Alfven time where the energy is converted into the axisymmetric spheromak state. Following this, the second stage consists of the slow growth of the spheromak state over the L/R time. Finally, the third stage is the steady sustainment of the equilibrium until the injectors are turned off. HIT-SI equilibria have been observed to be stable to n = 1 kinks and pressure-driven activity¹³ and are well described by the Imposed Dynamo Current Drive (IDCD) model.^14–16 

Of interest in the experimental operation of HIT-SI3 is the mode structure of the injector perturbation. While HIT-SI consisted of two injectors with primarily n = 1 toroidal symmetry, HIT-SI3’s three injectors (named A, B, and C) can be phased in ways to create different toroidal spectra of the perturbation. This allows the testing of how different non-axisymmetric helicity injection schemes interact with the spheromak. The three primary injector phasings that will be discussed in this paper are described in Table I. Helicity injection $K̇inj=2(VAΨA+VBΨB+VCΨC)$ can be kept constant or oscillated with these injector phasings.

Extended MagnetoHydroDynamic (eMHD) models have seen success at modeling the formation of plasma equilibria via helicity injection. These include CHI on the SSPX,¹⁷ HIT-II,¹⁸ and NSTX^19,20 devices and LHI on the Pegasus Toroidal Experiment.²¹ In the past, numerical models implemented by Izzo,²² Akcay,²³ and Hansen^24,25 have succeeded at simulating a subset of HIT-SI’s parameter space, specifically the low f_inj, low-β space. This study attempts to model the operation of HIT-SI3 at the same injector frequency as previous HIT-SI studies to match parameters for comparisons with a series of experimental discharges taken in the summer of 2016.

The motivation for validation of these models is twofold. First, numerical simulations are able to provide far more detail in the dynamics of the device than typical experimental measurements, particularly in plasmas with strong 3D perturbations such as HIT-SI3 with $(δBB)2∼6%$⁠. Second, validation at a variety of parameter spaces provides predictive capability for understanding deviations away from the experimental setup, which helps guide models used for predictive applications.

This paper will be organized as follows: Section II introduces the model used to simulate the HIT-SI3 device along with a description of a typical discharge. Section III contains a set of comparisons between experimental and simulation measurements, with an emphasis on the global structure of the plasma. Section IV is a discussion of results learned from the study of HIT-SI3. The paper concludes with a summary of these results.

HIT-SI3 is modelled using the NIMROD code²⁶ which solves the equations of extended magnetohydrodynamics,

where n is the particle density, $V$ is the center-of-mass velocity, T is the temperature, $B$ is the magnetic field, and $E$ is the electric field. Additionally, the following are the relations between terms in the model:

Two forms of this model have been implemented for understanding HIT-SI3, referred to as the zero-β and finite-β models. The zero-β model, which has been used in past studies by Izzo²² and Akcay,²³ uses the approximations $∇n=0$ and $∇T=0$ and only evolves the magnetic and momentum equations, setting right-hand side of Eqs. (1) and (3) to zero. In the zero-β case, two models for the viscous stress tensor $(Π)$ are compared, involving isotropic and anisotropic viscosity. The finite-β model relaxes this approximation and allows the full evolution of the fields and uses the anisotropic form of $Π$⁠.

The inclusion of the Hall term in Ohm’s law was guided by previous modeling of the HIT-SI device²³ which found that resistive MHD (rMHD) models fail to accurately describe the evolution of I_tor when compared with experimental results with HIT-SI. HIT-SI simulations found that without the Hall-term in the generalized Ohm’s law, the formation of the spheromak is delayed and the current amplification (⁠$Itor/Iinj$⁠) is lower than experimentally measured. In the magnetic field advance, the electron-inertia term (⁠$mene2∂J∂t$⁠) is used for numerical reasons to ease the solution of the advance involving the Hall term. An enhanced electron mass (⁠$me,simulation=100me,physical$⁠) is used to strengthen this term. Previous studies^23,24 have found little effect on the final solution from the addition of this term for HIT-SI simulation.

The diffusion coefficients are determined from several sources. The particle diffusivity D = 1000 m²/s is used purely for numerical stability. It is found that when operating with a lower diffusivity two problems occur: the first is a pumping out of the density from the volume where the total particle inventory is lost to the Dirichlet condition of the wall, while the second is large density oscillations on short length scales. The resistivity $η=η0T−3/2$⁠, with $η0=461 μΩ$ m, is determined from Spitzer assuming a constant log $(Λ)$⁠. In order to match the current gain more closely with the experiment, the zero-β models use $η=10 μΩ$m, corresponding with T = 12.8 eV. The viscosity values used are obtained from Braginskii²⁷ assuming plasma parameters late in the discharge. In the isotropic case, $ν∥=ν⊥=300$ m²/s, while in the anisotropic and finite-β cases $ν∥=11ν⊥=3300$ m²/s. The thermal conductivity coefficients are determined from Braginskii²⁷ in the same method used by O’Bryan for similar modeling of the Pegasus experiment.²¹ The cross-field heat flux is assumed to be dominated by the ion contribution while the parallel heat flux is dominated by the electron contribution. $k∥$ is computed using the enhanced electron mass mentioned above, which has the result of damping $k∥$ by a factor of $me,physicalme,simulation$⁠, which is $110$ in this case. Limited simulations using $me,physical$ have been performed but were insufficient for two reasons. First, the model becomes very stiff numerically because of the large non-axisymmetric fields present in the device (⁠$δBB∼25%$⁠), significantly increasing the computational time required. Second, due to chaotic field paths, a larger $k∥$ transports excessive thermal energy to the edges of the chamber which leaves the volume. This leads to a colder (T < 6 eV) plasma, exhibiting a resistivity that is too large for a spheromak to form, disagreeing with the experimental results.

NIMROD discretizes the equations in cylindrical coordinates, with the R-Z plane being composed of finite elements and the $ϕ$ component being a finite Fourier series. Using convergence properties determined by Akcay,²³ a grid with 28 × 28 finite elements of polynomial degree 4 is used with 22 Fourier components. The geometry neglects fine experimental features, such as openings in the flux conserver for diagnostics but captures the overall shape well. The outermost layer of finite element cells is a thin region (δ = 1 mm) kept at a large resistivity $ηwall∼105ηplasma$ to simulate the insulating layer on the experiment.

Due to the 2-D grid nature of NIMROD, the injectors cannot be directly modeled in the simulation. The injectors are instead implemented as boundary conditions on $B⊥$ and $E∥$ of the surface to create the ψ_inj and I_inj seen experimentally. Figure 2 shows both the grid and the locations of the applied boundary conditions. A detailed description of the implementation of these boundary conditions can be found in Ref. 23.

The boundary conditions for n and T are held constant on the edge of the domain, at $n=2×1019$ m⁻³ and T = 2 eV. Because the injectors contain inductively driven plasma currents, it is expected that the T_inj > T_wall on the surface of the domain. To account for this, a T profile with the same shape as $|B⊥|$ is applied on the boundary of the domain with $Tmax,inj=13$ eV.

Each simulation begins from the plasma breakdown time of the experimental shot it is attempting to replicate and uses waveforms for I_inj and ψ_inj taken from the discharge and fit to the following function where f(t) and $δϕ(t)$ are piecewise linear functions, with one point per injector period:

A comparison of the driven injector fields with the experimental reference shot for each phasing can be seen in Fig. 3.

Like HIT-SI, a discharge of HIT-SI3 consists of three main stages, which can be seen in Fig. 4. The first stage is the build-up of the non-axisymmetric injector dominated state. While in HIT-SI, this state was primarily n = 1 in structure, and in HIT-SI3 the injector phasings determine the Fourier spectrum of this state. The second stage occurs once a threshold value in the energy of the injector state is met and a relaxation process begins to form the spheromak. The spheromak forms on a fast time scale relative to the L/R time of the plasma, reaching a current amplification of 1 within a few Alfven times. This formation event is a global reconnection event the occurs at the Sweet-Parker timescale, with the Lundquist number in the region being $100<S<1000$⁠. Current filaments form connecting the different injectors together and induce these reconnection events, which play a role in the dynamics similar to those seen in CHI on NSTX^28,29 or LHI on Pegasus.²¹ Due to the 3-D and time-dependent nature of the injector fields, a more detailed study of the relaxation event is a topic for future work and is not discussed in this article. The third stage of the discharge is the slow growth and sustainment of the spheromak, eventually reaching a current amplification of between 2.5 and 3, depending on the parameters. Eventually, the third stage is ended by shutting off the injectors, marking the beginning of the decay stage of the spheromak. Figure 4 shows the evolution of these three stages for each injector phasing, where the differences in the applied non-axisymmetric perturbation can be clearly seen. Figure 5 expands on this, where we see that the poloidal flux profile of the axisymmetric component of the magnetic field forms a similar profile in each case, even with large differences in the perturbation spectra. Late in time, spheromaks with Lundquist number $S∼104$ are observed.

Due to the oscillating nature of the injectors, during all of this activity reconnection events in the vicinity of the injectors occur twice per injector cycles as the fields reverse. Experimentally, these events are thought to occur inside the injector while the boundary conditions of NIMROD mandate they appear in the confinement region.

The temperature evolution of the simulation expands further on the differences between stage 1 and stage 3 of the discharge. As seen in Figs. 6–8, the temperature profile during the injector-dominated stage of the simulation is defined by current channels that travel between the various injector mouths. These current filaments are thought to play a significant role in the relaxation dynamics. During the period of the simulation when the hot regions of the plasma connect to the injectors, we note that the hot-regions connect the various injectors together. Figure 8 shows this most clearly, where the hot filaments connect each injector mouth to another injector. During the spheromak sustainment stage, the profile is more axisymmetric, with a hot channel that is disconnected from the injectors.

Validation of these simulations is done by comparing measurements of synthetic diagnostics with experimental measurements, focusing primarily on comparisons of magnetic signals. A total of nine simulations are analyzed, with three simulations performed at each injector phasing of interest. A typical experimental discharge is selected for each injector phasing, then two zero-β simulations with isotropic and anisotropic viscosity models and one finite-β simulation with the anisotropic viscosity model were performed. Additionally, validation is done with velocity measurements from HIT-SI3’s Ion Doppler Spectroscopy (IDS) system and density measurements a Far-InfraRed (FIR) interferometer. The diagnostics used for the experimental measurements are discussed in more detail in Refs. 30–32 and 33.

There are two sets of probes used for magnetic measurements on HIT-SI3, as seen in Fig. 9. The surface probes³⁰ are used to measure the poloidal and toroidal dependence of B, as well as measure the toroidal plasma current, while the internal probe³¹ is used to measure the radial profile.

The first point of comparison for the magnetics is the global evolution of the spheromak. Experimentally this is examined by looking at the toroidal plasma current, measured by creating Amperian loops using the four poloidal arrays of surface probes. Shown in Fig. 10 is I_tor for both the experimental and simulation shots. As can be seen, each injector phasing exhibits different trends with respect to the viscosity model. While the 0–120-240 phasing sees $∼10%$ increase in I_tor in the anisotropic viscosity case, the 0–120-60 phasing sees little difference and the 0–0-0 phasing sees a $∼10%$ decrease in I_tor. Additionally, for all three injector phasings, the finite-β case sees I_tor reduced early in time and enhanced later compared with both the experiment and zero-β simulations. This is consistent with $η∼T−3/2$ as the plasma heats up the gain increases.

Another point of interest is the magnetic profile of the plasma for these three injector phasings. The toroidal mode structure of the plasma is obtained by a Fourier decomposition of the 16 midplane toroidal probes, allowing a decomposition of n = 0–7. For reducing the effect of the geometric shaping in the area on the analysis, the decomposition of the toroidal magnetic field is used. Direct comparison between simulation and experiment of the toroidal mode structure of the plasma measured by these probes is seen in Fig. 11. There are fine experimental features of the device near these probe locations which are expected to account for the differences in absolute magnitude seen between the simulation and the experiment. In particular, there are a number of holes in the flux conserver for midplane diagnostic access, which can lead to flux loss that the simulation does not capture.

The surface probes can also be used to calculate the current centroid of the spheromak. At each poloidal set of probes, the centroid can be calculated of the poloidal field. The formula used to compute this centroid is given by

This calculation can be performed at each of the four toroidal locations that have a poloidal array of probes. By comparing these four sets of probes, we can also obtain a description of the toroidal symmetry of the spheromak object. These comparisons are seen in Fig. 12. We can see that the finite-β calculations capture a ∼1 cm shift in the location of this center that the zero-β calculations do not.

The perturbations in the core of the device and the internal magnetic profile of the spheromak can be directly compared using the internal magnetic probe. The magnetic field signals from these probes can be seen in Fig. 13, with agreement seen on the phasing of oscillations at the injector frequency and disagreement seen with the magnitude of the field near the magnetic axis. The spheromak magnetic profile seen in Fig. 14, with an equilibrium fit to Grad-Shafranov following the method used in Ref. 13, is obtained by averaging the signal from the probe at f_inj.

In addition to magnetic comparisons, the fluid properties of the plasma can be compared. The IDS system is used to measure the flow velocity (⁠$V$⁠) and the ion temperature (T) in a chord-averaged sense. A Far-InfraRed (FIR) interferometry system is used to measure the line averaged electron density (⁠$〈n〉$⁠). These diagnostics are located on the midplane and Fig. 15 shows the chord locations relative to the injector mouths.

Previous investigations with the IDS system on the HIT-SI3^32,33 device have focused on measurements of the axisymmetric toroidal flow, phase coherence of oscillations at f_inj, and the temperature profile. As seen in Fig. 15, two sets of IDS chords are positioned looking in opposite directions of the plasma. An absolute toroidal flow can be calculated from the difference in mean velocity at a given impact parameter. Due to limitations on data gathered during experimental operations, comparisons can only be done with the 0–120-240 phasing. The comparison of this toroidal flow from simulations with experimental results is seen in Fig. 16. We see little difference in the axisymmetric flow profile observed by the three models, instead the different models change the amplitudes of the non-axisymmetric flows. As discussed later, the flows are dominated by $J×B$ forces related to the injector perturbation.

The temperature profile between the experiment and the simulation can be compared with results from the IDS system. Shown in Fig. 17 is the comparison of the measured and simulated temperatures late in time of the discharge. We see that NIMROD resolves a much flatter and smooth profile than the experimental measurement, potentially related to the spatial variation of the spectral intensities. We can similarly compare the effect of line-averaging by comparing the synthetic IDS signal with the radial temperature profile in NIMROD. Figure 18 shows that while NIMROD observes radial dependance on T, the line-averaged profile is flat. Of note is that while the simulations begin cold with $T=Twall$ everywhere and heat up throughout the discharge, the experiment tends to observe a relatively flat T_ion evolution in time. This can be seen in Fig. 19. The simulation heats up to the temperature observed by the CIII spectral line throughout the length of the discharge, indicating reasonable modeling of the total thermal energy of the plasma.

The FIR measurements³¹ of $〈n〉$ can be extracted from the three finite-β simulations and the results can be seen in Fig. 20. Unfortunately, due to the large particle diffusivity used for numerical convenience, the oscillations are significantly damped in the simulation. A possible improvement would be the implementation of a boundary condition on $V⊥$ to refuel the particles lost to the wall when diffusivity is decreased.

As can be seen in the current centroid and internal magnetic probe comparisons, the location of the magnetic axis undergoes a shift when pressure gradients are included in the model. This shift improves agreement with the experimental measurements indicating that these pressure gradients have a significant effect on the plasma. We can calculate $〈β〉vol=∫pdV∫B22μ0dV$ and find that at each injector phase the simulations exhibit $〈β〉vol∼5%−8%$ towards the end of the discharge, with the value oscillating at f_inj. This can be compared with the experimental value obtained from equilibrium fits to the Grad-Shafranov equation in Fig. 14, using the same method in Ref. 13, which tend to give values between $10% and 15%$ late in the discharge. Additionally, the value of $〈β〉vol$ obtained from taking synthetic experimental measurements from the simulation and fitting to a Grad-Shafranov equilibrium is close to the value obtained from volume integrating the fields directly, giving confidence in the usage of equilibrium fits on these time-averaged fields.

This disagreement of a factor of two in $〈β〉vol$ could be caused by a combination of increased particle diffusivity, smoothing the pressure profile and the simplified model of assuming T_i = T_e. Experimental evidence suggests that $Ti∼2Te$ though numerical studies thus far have proven too computationally difficult to carry out.

An interesting feature observed when comparing the toroidal mode activity seen from both the experimental and simulation results is the large difference in spectrum between the total volume energy and the surface probe measurement. In particular, while the 0–120-240 phasing sees primarily n = 1 perturbation activity from the surface probes in both experiment and simulation, there is an equivalent magnitude n = 2 shaped perturbation concentrated closer to the core of the device. In addition, the 0–120-60 phasing shows a similar discrepancy concerning the magnitude of n = 2 activity relative to other modes. As would be expected, when the phasing is dominant in only one toroidal fourier mode, as in both HIT-SI (n = 1) and 0–0-0 HIT-SI3 (n = 3), this effect is not observed and the spectrum of the surface probes matches the spectrum of the volume integrated energy.

Additionally, the comparisons of magnetic profile between the experiment and the simulations show the importance of finite-β effects. An outward shift of the magnetic axis is seen between the zero and finite-β simulations, which are consistent with $β∼10%$⁠. An analogous outward shift from rotation would require $V∼100$ km/s which is not observed in either the experiment or the simulation. While these effects lead to a change in the magnetic profile, the evolution of global parameters (such as I_tor) changes little, resulting primarily from changes in energy dissipation from the plasma resistivity. This indicates that the mechanism for creation and sustainment of the axisymmetric component of the equilibrium is a zero-β phenomenon.

The additional diagnostic capability provided by the simulation allows us to examine questions that have been raised in previous studies of HIT-SI and make predictions for future measurements. The magnetic profile of the equilibrium is constructed by taking a running average of the internal magnetic probe at f_inj. Because of the large non-axisymmetric nature of the perturbation and the probe only located at one toroidal location, it is unclear whether this time-averaged field is axisymmetric. We are able to produce synthetic internal magnetic probes throughout the device and see the deviations this time-averaging has from axisymmetry. The results are that the radial position of the magnetic axis varies by ±1 cm as can be seen in Fig. 21. This is consistent with the toroidal deviations seen in the current centroid measurement from the surface probes and can be verified with an additional internal magnetic probe that has been installed for the next run campaign.

While IDS comparisons were only made for 0‐120-240 phasing, we can make predictions about the axisymmetric flow profile expected at the other phasings. Figure 22 shows the n = 0 component of $Vtor$ during the sustainment phase. In all three phasings, the angular velocity $Vtor≪ωinjR≈25$ km/s at R = 0.3 m, suggesting that there is no plasma rotation at f_inj. The 0‐120-60 phasing appears to have the most significant variation based on viscosity model, with the anisotropic model having a flat angular velocity profile at $frot≈3.5$ kHz.

As seen the IDS flow velocity comparison with Fig. 16, the magnitude of axisymmetric flow does not depend strongly on the viscous model used. The axisymmetric flow in the simulation ends up being a minor component of the ion flows, with the majority having a toroidal spectrum similar to the injector magnetic perturbation, as seen in Fig. 23. We see that while the anisotropic viscosity model damps away $n≠1$ kinetic energy from the system, the finite-β model more closely resembles the isotropic kinetic energy spectrum, despite using the anisotropic viscosity model.

While no current experimental diagnostics exist to measure the temperature at specific locations in the device, we can look at the profile produced for each injector phasing in the simulation. Figure 24 shows the temperature evolution and profile late in time, where temperatures of T > 30 eV are seen at the magnetic axis. From these, we see that the temperature of the plasma evolves similarly to the toroidal current in the simulations, which is consistent with the heating power being proportional to $Itor2$ from resistive heating.

The primary result of this study is the understanding that (1) HIT-SI3 is able to form an axisymmetric spheromak with any non-axisymmetric helicity injection mode structure and (2) NIMROD is able to capture the general evolution of the formation and sustainment of HIT-SI3 plasmas. These sustained spheromaks show hints of energy confinement, with T_max existing at the magnetic axis during the high gain portion of the simulation. Additionally, while there are changes to the flow profile between isotropic and anisotropic viscosity models, they tend to have little influence on the bulk evolution of the plasma. The primary points of agreement lie in the toroidal modal structure of the perturbation, and measurements of fluid velocity and ion temperature. The main point of disagreement with the experiment lies in the equilibrium profile of the poloidal magnetic field, notably the outward shift that is related to increased β. It is likely that the approximations used for the density and electron temperature contributed to the reduced β seen in the simulation. Future measurements that will serve as a test of the predictions of these NIMROD simulations include a second internal magnetic probe to provide measurements of the non-axisymmetry of the field, a Thomson scattering diagnostic to provide measurements of T_e, and further usage of the IDS diagnostic at more injector phasings.

Comparing the different injector phasings to see which is better, we see that the two constant helicity injection cases appear to outperform the 0‐0-0 phasing in both current gain (⁠$Itor/Iinj$⁠) and temperature, which is consistent with experimental results.¹² The two constant helicity injection cases (0‐120-240 and 0‐120-60) reach similar performance levels, with current gain $∼2.5$ and peak temperatures T > 30 eV.

Future numerical studies to expand upon this work are likely to focus on the modeling of the sustainment phase of the discharge, with emphasis on the evolution of n and T. In particular, the additional particle diffusivity applied in the evolution of n likely has unphysical consequences on energy transport, so this will be a focus of future studies. In addition to this, experimental evidence indicates that $Ti≈2Te$⁠, suggesting that a two-temperature model may be more appropriate. Additionally, verifications of the results of these simulations can be compared with the PSI-Tet code, a code which is able to solve the same equations with an expanded geometry to include dynamics inside of the helicity injectors. These dynamics may play a role in the global structure of the non-axisymmetric perturbations central to HIT-SI3’s operation.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award No. DE-FG02-96ER54361. Simulations presented here used resources of both the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system and funded by the STF at the University of Washington.

Additionally, many thanks to members of both the HIT team and the PSI-Center who have provided multiple discussions and insight into the interpretations of HIT-SI3 results.