Integrated modeling of high β_N steady state scenario on DIII-D

The normalized plasma pressure β_N is a key figure of merit for high fusion gain operation in steady-state tokamaks,¹ where β_N = β_TaB_T/I_p, a (m) is the minor radius, B_T (T) is the toroidal magnetic field, I_p (MA) is the plasma current, $βT=2μ0〈p〉/BT2$ is the toroidal β, and $〈p〉$ (Pa) is the volume average of pressure. Steady-state reactors must sustain plasma current non-inductively without a transformer and minimize the recirculating power needed for external current drive. The fusion gain Q, which is defined as the ratio of the fusion power to the auxiliary heating and current drive power, is proportional to β_N and the self-generated bootstrap current fraction (f_BS) scales² with β_N. High β_N operation at high f_BS reduces the necessary fusion power for the same net electricity generation if energy confinement is good enough, which alleviates material challenges such as neutron wall loading and divertor heat flux control. Fusion power plant studies based on steady-state tokamak operation suggest that β_N in the range of 4–6 is needed for economic viability,³ while the optimum choice of β_N may depend on B_T.⁴ 

High β_N steady-state operation needs simultaneous optimization of the plasma current and pressure profiles. The upper limit to β_N is determined by the ideal MHD stability limit of low-n kink modes, which is a strong function of the current and pressure profiles and plasma shape.^5–9 Previous extensive modeling and experimental studies also indicate that optimizing the current profile is a primary tool to reduce turbulent energy transport resulting from microinstabilities.^10–13 Not only the current profile but also the local gradient of the pressure profile modifies the transport characteristics.^14,15 The relative effects of the local pressure gradient on confinement depend on the current profile as well.¹¹ There exist significant challenges in the simultaneous optimization of the current and pressure profiles since they are linked tightly to each other through the nonlinear dependency of the local transport and bootstrap current to the total plasma current and pressure profiles recursively.

DIII-D is investigating a range of potential high β_N steady-state solutions including high q_min,^5,16 elevated q_min,^17,18 steady-state hybrid,^19,20 and high $ℓi$⁠,^21,22 which are distinguished primarily by the plasma current profile. The high q_min scenario employs a very broad current profile with the minimum safety factor q_min > 2 and utilizes weak positive or negative magnetic shear in the plasma core in order to improve energy confinement. β_N is sustained well above the no-wall ideal-stability β_N-limit by placing more plasma current close to the conducting wall to maximize coupling and wall stabilization of low-n ideal modes. The other extreme is the high $ℓi$ scenario, which relies on a very peaked current profile with $qmin≃1$ for the improved confinement and stability but in different ways from the high q_min scenario. The large positive magnetic shear from the peaked current profile leads to enhanced confinement and results in a high no-wall β_N limit. The high β_N condition can be achieved without the need for the wall stabilization, which makes the high $ℓi$ scenario attractive for the power plant reactors.²³ In both scenarios, simultaneous optimization of the current and pressure profiles is essential to realizing the high β_N operation, which determines requirements for the external heating and current drive to sustain the desired profiles.

Scenario modeling plays a key role in the development of self-consistent high β_N steady-state scenarios on DIII-D in a repeated cycle of scenario design, experimental implementation, and modeling validation. Recently, remarkable progress has been made in improving the individual elements of the integrated scenario modeling, such as a comprehensive theory-based core transport model, TGLF²⁴ for all transport channels (particle, energy, and momentum), and EPED^25–27 for edge pedestal to provide the boundary condition of the core transport. Integration of such high fidelity core transport and edge pedestal models along with the well-established modeling of ideal kink stability and external heating and current drive is an important step towards developing a fully predictive scenario modeling tool with a limited number of free parameters.

This paper describes a promising path to a self-consistent $βN∼5$ steady-state equilibrium found by the theory-based integrated scenario modeling. In Sec. II, a new iterative numerical procedure is introduced, which integrates theory-based models of core transport, edge pedestal, heating and current drive, equilibrium, and stability to self-consistently find the target high β_N steady-state (d/dt = 0) solutions. Experimental validation against DIII-D experiment is presented for high β_N, high q_min discharges in stationary conditions. In Sec. III, we discuss optimization of the current and pressure profiles with an emphasis on inherent correlation between the current and pressure profiles. Section IV shows simulations to determine the detailed heating and current drive requirements for production of an MHD stable β_N ∼ 5 scenario based on the optimization direction identified in Sec. III. Section V summarizes the main conclusions.

A new scenario modeling tool, IPS-FASTRAN, is an efficient and robust iterative solution procedure to find a steady-state solution (d/dt = 0) of core transport, self-consistent with external heating and current drive (H/CD), MHD equilibrium, ideal MHD stability, and edge pedestal. The IPS-FASTRAN modeling has been benchmarked successfully against time-dependent simulations.²⁸ Time-dependent simulations starting with a range of different initial conditions converge to the same final state found by the IPS-FASTRAN modeling. Although the system of equations is highly nonlinear, no bifurcation in the final state depending on initial conditions has been observed so far. The steady state solution procedure requires only a few iterations, typically less than 5, which is equivalent to roughly 10 time steps of time-dependent simulation in terms of computational cost. Current profile relaxation time estimated by $τR=0.17R0/R$²⁹ is about 2 s for typical DIII-D high beta discharges, where R₀ is the major radius in meters and $R$ is the plasma resistance in $μΩ$⁠. The time-dependent simulation to reach close to the fully relaxed state, at least $2τR$ (∼4 s) in a time step of 10 ms, requires an order of 400 time stepping. So, the steady state solution procedure is extremely efficient for the target scenario development. Note that the current profile relaxation time for burning plasmas like ITER is even many orders longer than DIII-D high β_N plasmas.

IPS-FASTRAN has been developed based on a modern integrated modeling framework, Integrated Plasma Simulator (IPS),^30,31 which couples a modular 1-D radial core transport solver, FASTRAN³² using theory-based transport models with a set of external H/CD codes, MHD equilibrium, and stability codes. The IPS modeling employs a framework/component architecture using existing codes and file-based communication with Plasma State for data transfer mechanism between components. The IPS framework provides a flexible capability for the new workflow development and accelerates the repeated cycle of modeling, experimental validation, and scenario design/development. The IPS-FASTRAN scenario modeling employs the DAKOTA-enabled IPS³⁰ for the target design optimization and sensitivity analysis to execute multiple simulations simultaneously with efficient high performance computing (HPC) resource utilization on massively parallel computer. The modeling presented in this paper employs TGLF²⁴ for core transport model with the SAT0 saturation rule,³³ EPED1^25–27 for edge pedestal, NUBEAM³⁴ for Monte-Carlo beam ion slowing down calculation of neutral beam injection (NBI), TORAY-GA³⁵ for ray-tracing modeling of electron cyclotron heating (ECH) and current drive (ECCD), EFIT³⁶ as MHD equilibrium solver, and DCON³⁷ for ideal MHD stability calculation. The Chang–Hinton³⁸ and Sauter³⁹ models are employed for the neoclassical transport and bootstrap current calculation, respectively.

The steady-state solution procedure at each iteration includes the following sequence:

1. The modular transport solver, FASTRAN,³² finds a steady-state solution of electron density (n_e), electron temperature (T_e), ion temperature (T_i), and toroidal rotation (Ω) with turbulent radial fluxes predicted by TGLF for a given MHD equilibrium and source profiles, where a simple nonlinear relaxation method derived by root-finding algorithms is employed to solve the highly nonlinear set of the 1-D diffusion-type transport equations. Boundary conditions are applied at the pedestal top, $ρ=1−3/2wped≡ρpedtop$ with the values predicted by IPS-EPED1, where ρ is the normalized minor radius proportional to the square root of the toroidal flux and $wped$ is the full width of the edge pedestal.

2. IPS-EPED1⁴⁰ is an HPC-enabled nested IPS workflow called by the parent workflow of the IPS-FASTRAN steady-state solution procedure, replicating the original EPED1^25–27 to predict the width and height of the pedestal. The electron density profile for $ρpedtop≤ρ≤1$ is taken from the EPED1 model profile with a hyperbolic tangent shape in the pedestal. The value at the separatrix of ρ = 1 is assumed $nesep=1/2 neped$⁠, where both $neped$ and $nesep$ are inputs of EPED1. The temperatures for $ρpedtop≤ρ≤1$ are updated by assuming

   $neTe=niTi=12pEPED1,$

   where $pEPED1$ is the total pressure predicted by EPED1. Note that the ion density n_i is calculated by the charge balance with the input value of $Zeff$⁠.

3. The plasma current profile is fully relaxed in the entire region inside of the separatrix using FASTRAN by solving the steady-state solution of poloidal magnetic flux evolution equation with given external driven and bootstrap current profiles, where the nonlinear dependency of the bootstrap current on the total current profile is taken into account.

4. The source profiles from the external heating and current drives are calculated by the NUBEAM³⁴ and TORAY-GA³⁵ components. The NUBEAM Monte-Carlo calculation of the beam ion slowing down is time-dependent and reaches to the stationary state within a few slowing down time. In order to minimize Monte-Carlo noise in the calculated beam-driven current and source profiles (particle, torque and energy), an average of the time-dependent profile is taken for about one slowing down time in the stationary state.

5. MHD equilibrium is updated using EFIT³⁶ with the current and pressure profiles as input. The pressure is the sum of the thermal pressure calculated in step 1 and the beam ion pressure calculated in step 4.

For robust convergence in the steady-state solution procedure, a simple Picard-type under-relaxation is employed between iteration: $fk=αf*+(1−α)fk−1$⁠, where f denotes unknown such as n_e, T_e, T_i, Ω, the superscript k is the iteration step, and f* is the solution obtained in step 1.

The IPS-FASTRAN steady-state solution procedure reproduces most features of DIII-D high β discharges in the stationary condition. Figure 1 shows an example of IPS-FASTRAN application to stationary discharge with q_min > 1.5 at $βN∼3$⁠, which was sustained by the off-axis NBI and ECCD longer than $2τR$⁠.^5,16 The calculated radial profiles of n_e, T_e, T_i, and Ω reproduce the measurement reasonably well as shown in Fig. 1. The predicted plasma current profile [Fig. 1(e)], to which the discharge would eventually evolve, is close to the estimation by kinetic EFIT equilibrium reconstruction obtained with MSE and the calculated edge bootstrap current as constraint. The discrepancy in the current profile (and q) near the magnetic axis between the measurement and the modeling indicates that the inductive current would penetrate further inside when the loop voltage reaches a fully relaxed condition. These simulations are largely based on the theory-based models with the adjustable input free parameters limited to the choice of

where $Ωpedtop$ is the plasma rotation at the pedestal top and D_b is the anomalous beam ion diffusion coefficient. For the simulation shown in Fig. 1, a uniform anomalous fast ion diffusion of $Db=0.3 m2/s$ was used in the NUBEAM calculation to make the calculated total stored energy match the stored energy determined by equilibrium reconstruction.

The toroidal field ripple may affect the fast ion confinement and toroidal rotation.⁴¹ In this work, the ripple effects are not taken into account. The high n non-axisymmetric perturbed magnetic field tends to be localized at the plasma periphery,⁴¹ where the boundary condition for toroidal rotation is applied and the beam ion density is very low. Including the ripple effects in the IPS-FASTRAN modeling, however, is left for future study.

DIII-D has demonstrated fully non-inductive operation (⁠$fNI≡INI/IP=1$⁠) with the elevated q_min scenario. An example of such scenario reported in Ref. 18 employs relatively broad current profile at $qmin∼1.5$ with insufficient off-axis current drive to sustain higher q_min > 2. This discharge was operated at $βN∼3.8$⁠, well above the n = 1 no-wall β_N stability limit but close to the ideal-wall β_N stability limit. This section describes a promising path for optimization of the pressure and current profiles identified by the theory-based integrated modeling in order to increase β_N toward a DEMO reactor relevant $βN∼5$ by improving both confinement and stability.

The modeling predicts that higher thermal energy confinement is obtainable by broadening the current and pressure profiles. Figure 2 shows the calculated profiles of the electron and ion temperatures and the thermal diffusivities for a range of different current profiles by varying minimum q (q_min) and its location (⁠$ρqmin$⁠), otherwise in the same conditions such as B_T, I_P, and injection power. The integrated scenario modeling for the steady-state solution procedure described in Sec. II is used with the TGLF model for core transport along with the EPED model for edge pedestal. However, in order to isolate the effects of the MHD equilibrium on thermal confinement, the plasma current profile is prescribed rather than calculated self-consistently by solving the poloidal magnetic flux evolution to steady-state. This also mimics the transient target current profile formation in limited pulse length high beta experiments. In the edge pedestal (⁠$ρ>ρpedtop$⁠), the model current profile is given by the bootstrap current. In the core region, the current density profile is characterized by the peak location of the current density ρ₁ and the ratio $J(0)/J(ρ1)$ to make a variation in q_min and $ρqmin$⁠, as shown in Fig. 2. The simulation conditions are B_T = 1.6 T, I_p = 1.2 MA, and q₉₅ = 5.2 in a double null shape with R₀ = 1.7 m, a₀ = 0.6 m, κ  =  2.0, and δ = 0.6. A total NBI power of 15.6 MW is applied in this simulation, with 9.6 MW of this directed off-axis with the same geometry as described in Ref. 28. In this specific scan, the electron density is fixed with the experimental profile from one of the DIII-D high q_min discharges.

Systematic differences in the electron and ion temperature profiles are predicted as the current profile changes. The weak magnetic shear generated by the broader current profile leads to improved confinement with reduced χ_e and χ_i as shown in Fig. 2. Accordingly, the calculated β_N increases as $ℓi$ decreases [Fig. 3(a)] at the same injection power. Here we use $ℓi$ as a measure of the current profile broadness. The broader current profile also results in broader T_e and T_i. The pressure peaking, which is defined as $fP=p0/〈p〉$⁠, increases with $ℓi$ as shown in Fig. 3(b).

The correlation between the current and pressure profiles is evident from the experiment reported in Ref. 5. Figure 4 plots the pressure peaking as a function of $ℓi$⁠, showing that the discharges with broader current have broader pressure. This correlation appears to result from the internal transport process as shown in Fig. 2 as well as the way in which the broader current is sustained in the experiment. For example, off-axis NBI for the generation of the broader current produces broader beam ion pressure, leading to the reduced total pressure peaking. The modeling shows that broader plasma current inherently broadens the thermal pressure profile, which is independent of the heating location and the beam ion pressure profile, because the heating profile is the same for all cases shown in Figs. 2 and 3.

Figure 3(a) shows the low-n kink mode ideal MHD β_N stability limit with conducting wall calculated by DCON for a series of equilibria with variations of current profile, otherwise in the same condition. The minimum value of the β_N limit for $1≤n≤3$ was plotted, and the n = 1 limit is the lowest limit in most cases. The predicted β_N limit increases as $ℓi$ decreases since broader current places more plasma current close to the wall and maximizes the conducting-wall stabilization of low-n ideal modes. Broader pressure associated with broader current also plays an important role in the high β_N >5 limit at low $ℓi$⁠. Previous modelings with independent variation of the current and pressure profiles indicate that both broader current and pressure lead to higher β_N stability limit.^5,6 This means that broader current does not necessarily have a higher β_N stability limit if the pressure peaking is higher. Simultaneous optimization of the pressure and current profiles for high β_N is promising in high q_min scenarios since broader current tends to broaden pressure inherently. Techniques to sustain broad current self-consistently by external H/CD actuators are discussed in Sec. IV.

The goal of the modeling presented in this section is to determine the detailed heating and current drive (H/CD) requirements for the production of an MHD stable β_N = 5 scenario. Obtaining fully noninductive condition is challenging because it requires an exact balance among the confinement, heating, and current drive.⁵ We will describe the effect of possible modifications to the DIII-D neutral beam heating configuration along with increased ECCD power on discharge parameters to guide optimum choice of the H/CD system upgrade.

DIII-D is equipped with four positive-ion-based neutral beamlines with three co-current injection and one counter-current injection. Each beamline consists of two ion sources having different tangency radii. One of the beamlines has been modified for downward vertical steering to provide substantial off-axis current drive. The measured off-axis NBCD is consistent with NUBEAM modeling without anomalous fast ion transport for a range of beam injection and discharge conditions⁴² unless there are significant MHD activities present in the plasma such as Alfven eigenmodes (AE). The DIII-D experiments on off-axis neutral beam injection into high β_N discharges have demonstrated broadening of the current and pressure profiles, leading to the increased β_n limit determined by ideal low-n MHD instabilities.^5,16

Figure 5 shows the calculated β_N and the β_N stability limit in d/dt = 0 condition as a function of the ratio of the off-axis NBI power to the total NBI power, f_OANB = P_OANB/P_NB. Total 15.0 MW of the co-current NBI is injected, which is approximately the maximum available NBI power from three co-current beamlines. The highest f_OANB for the present NBI configuration is ∼0.3 for P_NB = 15 MW. An increased ECCD power P_EC = 9 MW is used aiming the peak ECCD location off-axis at ρ ∼ 0.6. The β_N value stays about constant with f_NI ∼ 1 for the variation of f_OANB, while the β_N stability limit increases rapidly with f_OANB mostly due to the broadened current and pressure profiles as shown in Figs. 5(b) and 5(c). Broader current achieved by more off-axis NBCD broadens the thermal pressure profile as well, as shown in Sec. III. In addition, the beam ion pressure profile is broader for higher f_OANB, resulting in further reduced peaking of the total pressure profile. Stable operation can be achieved at the β_N ≈ 4.5 if f_OANB > 0.4 for this specific condition. This modeling suggests that more co-current and off-axis NBI power are needed for stable operation at $βN≥5$⁠, which has motivated the planned modification of the present on-axis counter beamline to off-axis, toroidally co- or counter-steerable beamline.

A self-consistent fully noninductive operation space at β_N > 4 was explored by the IPS-FASTRAN scenario modeling with the planned upgrade of DIII-D heating and current drive tools—second off-axis beamline (toroidally steerable) and increase of ECCD power. Figure 6 shows the predicted discharge parameters as a function of q₉₅ at the same B_T = 1.6 T. Each data point in the scan satisfies an exact current balance at f_NI = 1. A maximum off-axis beam power of P_OANB = 9.6 MW is injected along with P_EC = 6 MW for off-axis current drive, and additional on-axis NBI power is supplied to maintain the necessary β_N value. It is predicted that lower q₉₅ requires higher β_N and more NBI power for f_NI = 1, which is consistent with the 0-D experimental scaling obtained from a range of the relevant high β_N DIII-D discharges.⁵ The β_N stability limit, however, decreases as q₉₅ decreases. This decrease is due to a combination of lower q_min, higher $ℓi$⁠, and increased $p0/〈p〉$⁠, as shown in Figs. 6(b)–6(d). The profiles are more peaked at lower q₉₅ because more on-axis NBI power (so lower f_OANB) is needed to maintain f_NI = 1. The highest stable β_N is ∼4.64 at q₉₅ = 6.07, which is just below the ideal β_N stability limit at 4.83.

The predicted profiles for the self-consistent solution at the highest β_N are plotted in Fig. 7. The current profile has a weak magnetic shear inside 0.2 < ρ  <  0.6 with $|q−qmin|<0.5$ and $qmin∼2.14$⁠. The narrow shear reversed at ρ < 0.2 could be eliminated if a small amount of seed current is provided near the magnetic axis. Compared with the high β_N, high q_min discharges produced with the present H/CD systems^5,16 (dominant ion-heating, T_e/T_i < 1), the volume average temperature ratio $〈Te〉/〈Ti〉∼1.02$ is close to 1 since the total external heating delivered to electrons and ions is comparable at 8.4 MW and 8.8 MW, respectively. Note that the ratio of electron heating to ion heating is greater than 4 where the ECCD is localized as shown in Fig. 7(f).

One of the main motivations for q_min > 2 is to eliminate the q = 2 rational surface to avoid the most harmful m/n = 2/1 tearing mode. Although there is no q = 2 surface in the solution shown in Fig. 7, operation at β_N close to the ideal stability limit may develop a low-n tearing mode such as m/n = 3/1, 5/2 modes, since the classical tearing stability index $Δ′$ can increase exponentially as β_N approaches to the ideal limit.⁴³ In contrast, the slightly higher q₉₅ ∼ 6.2 solution has a sufficient margin for the β_N limit, but with ∼4% lower operation β_N.

The high power ECCD provides far off-axis current drive to generate a weak magnetic shear q profile with q_min > 2, while still maintaining a reasonable current drive efficiency at the larger minor radius. The optimum location of ECCD (ρ ∼  0.6) is determined by maximizing the highest stable β_N. Figure 8 shows that increased ECCD power up to P_EC = 9 MW moves the highest stable β_N solution to lower q₉₅, extending the f_NI = 1 operation space at $βN>4.5$⁠. At the same q₉₅, higher ECCD power increases the β_N stability margin.

In high q_min DIII-D steady-state scenario discharges, Alfven eigenmode-type instabilities have been observed, resulting in additional fast ion transport which reduces the driven current and the power transferred to the thermal plasma.^5,44 Simulation and experiment show that moving the location of q_min to larger outer minor radius, which is an optimization direction for thermal confinement and ideal MHD stability as shown in Fig. 2, is effective in minimizing fast ion transport due to AE.⁴⁵ There is not yet a model integrated into IPS-FASTRAN that can provide a quantitative prediction of fast ion loss as a result of AE modes. All of the simulations presented so far in this paper assume a modest value of the anomalous fast ion diffusion D_b = 0.3 m²/s. Figure 9 shows a sensitivity of the predictions to the assumed D_b. The simulation conditions including q₉₅ and ECCD power are the same as the stable highest β_N solution (solution shown in Fig. 7). It is predicted that the highest stable β_N is not sensitive to D_b, while the required beam power increases with D_b.

The IPS-FASTRAN integrated scenario modeling validated against DIII-D experiment guides development of high β_N steady-state scenario. A new efficient steady-state solution procedure found a range of fully noninductive β_N > 4.5 solutions, self-consistent with core transport, edge pedestal, external heating and current drive, and low-n ideal MHD stability. This high β_N scenario with high q_min > 2 relies on very broad current density profile to improve both the ideal MHD stability and thermal plasma confinement. A strong correlation between the current and pressure profiles is found. Broader plasma current broadens thermal pressure, independent of the heating location and the beam ion pressure profile. This inherent correlation facilitates simultaneous optimization of the current and pressure profiles in the high q_min scenario.

The predictive modeling presented in this paper determines the detailed heating and current drive requirements for the production of an MHD stable $βN=5$ scenario, guiding the DIII-D upgrade path including the 2nd off-axis beamline and increased ECCD power. These planned upgrades will enable optimization of steady-state scenarios to better inform the development of steady state fusion toward DEMO.

This work is based upon work supported by the U.S. DOE, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility, under Award Nos. DE-AC05-00OR22725, DE-FC02-04ER54698, DE-AC52-07NA27344, DE-FG02-95ER54309, and DE-SC0012656. This research used resources of the National Energy Research Scientific Computing Center (NERSC), 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.