Computational mathematics and physics of fusion reactors

Equilibrium

Most configurations for magnetic fusion use toroidal geometry to confine the orbits of ions and electrons on closed flux surfaces inside the plasma. The Kolmogorov–Arnold–Moser theory gives conditions for the existence of nested families of these magnetic surfaces (4). It has been shown by similar methods that continuously differentiable solutions of the equilibrium problem cannot be found when the geometry of the magnetic field is genuinely three-dimensional (5). Nevertheless, stellarator experiments have become increasingly successful in confining fusion plasmas. Mathematical analysis shows that an appropriate response to this anomaly is to look for weak solutions of the equilibrium equations in a conservation form that captures current sheets and islands but does not require the magnetic field to be either differentiable or continuous (6).

The variational principle of magnetohydrodynamics (MHD) leads to a method of accelerated steepest descent to compute weak solutions of the equilibrium problem for both stellarators and tokamaks (7, 8). This has been implemented in a family of high-performance codes including BETA, VMEC, and NSTAB. The VMEC code is widely accepted by the fusion community and has become an integral part of the diagnostics in major stellarator experiments (9). Work remains to be done on free-boundary versions to reconstruct the plasma from observations and calculate solutions accurately enough to detect the presence of significant islands among the flux surfaces.

The codes have been run to study the LHD experiment, which has performed surprisingly well because of high shear near the edge of the plasma. They have played an important role in the design of the W7-X experiment at Greifswald in Europe and the National Compact Stellarator Experiment at the Princeton Plasma Physics Laboratory (10, 11). The codes are mathematical tools that have enabled physicists to exploit a concept of two-dimensional quasisymmetry of the magnetic field providing good drift surfaces of trapped particles. Moreover, bifurcated solutions with three-dimensional asymmetries have been calculated for tokamaks that simulate observed physical phenomena (12). In many cases, three-dimensional computations seem to model tokamak physics better than results from two-dimensional theory (12, 13).

Stability

The failure of differentiability in solutions of the three-dimensional equilibrium equations leads to difficulties in applying linear stability theory to stellarators because calculation of coefficients that are singular is required. The mathematics are harder than the two-dimensional analysis of tokamaks, and it is not always easy to distinguish between equilibrium and stability limits on β. Local estimates of stability such as the Mercier criterion or the ballooning theory turn out to be unrealistic for stellarators, and experiments have been found to exceed predictions by those methods (1–3).

Because it is based on the MHD variational principle, the NSTAB code furnishes tests for linear and nonlinear ideal MHD stability that do not depend on unjustified integration by parts and differentiation (8). An intuitively plausible mountain-pass theorem can be applied that asserts that, between any two minima of the energy landscape, a saddle point can be found representing an unstable equilibrium. Runs of the NSTAB-iteration scheme can be used to search for bifurcated solutions that indicate whether a configuration is in linearly or nonlinearly unstable equilibrium, because the unexpected appearance of several solutions of the problem establishes instability.

This method has been applied to the LHD experiment to explain why surprisingly high values of β were observed after the plasma passed through states that were linearly unstable but nonlinearly stable. Examples of such calculations over 5 or 10 field periods, with the pressure p given as a function of the toroidal flux s, are displayed in Figs. 1 and 2. The theory is also consistent with recent measurements (3) of β in the W7-AS experiment (compare with Fig. 3). It seems to be more realistic than computations of local criteria for the LHD displaying the dependence of a ballooning eigenvalue on the poloidal transit of a magnetic line that was examined (compare with Fig. 4).

Accuracy of the NSTAB code adequate to address finer issues of stability has been achieved by coupling a delicately arranged difference scheme in the radial flux coordinate s to a spectral treatment of dependence on poloidal and toroidal angles within the flux surfaces (7, 8). Convergence studies and comparisons with physical observations have been performed to validate the theory, and the method has been applied to design configurations for a new stellarator experiment (14).

Transport

The three-dimensional problem of transport in stellarators is not easy to handle by asymptotic expansions. However, it can be formulated quite effectively as a linearized Fokker–Planck equation in many independent variables. One can apply the Monte Carlo method in that analysis, which solves a drift kinetic equation by alternating between a random walk to model the collision operator and an ordinary differential equation solver to track the orbits of charged particles (15). Separation of variables shows that there are solutions that decay exponentially in time, and they lead to a fast algorithm determining the confinement times of the ions and electrons independently (16).

Quasineutrality requires that there be no difference between the confinement times and distribution functions of the two species. That can possibly be achieved at low collision frequencies by iterating in a systematic fashion on the choice of the electrostatic potential, which may be allowed to vary within the magnetic surfaces (13). A flow in the background could be adjusted to make the collision operator conserve momentum, but numerical experiments show that, if the flow is not thermal, it does not have enough effect on the collision operator to account for quasineutrality (1). Small oscillations of the electrostatic potential along the magnetic lines serve that purpose and also model the ambipolar electric field and turbulence tolerably well. A parallel version of the TRAN computer code implementing this concept has been programmed for the IBM Special Purpose computer at the National Energy Research Computing Center by using openmp software.

Recent high-temperature shots in the LHD experiment have produced observations of the energy confinement time at low collision frequencies that can be simulated by the Monte Carlo calculations we have described (2). Careful comparisons establish that the numerical method gives realistic results. Oscillations of the electric potential along the magnetic lines model anomalous transport, so there is good agreement with measured values (1). Similar agreement with experimental data for tokamaks is obtained when bifurcated equilibria without perfect two-dimensional symmetry are used in the computations (13).

For axially symmetric tokamaks, steady solutions obtained from asymptotic expansions are used in the two-dimensional theory of transport to predict intrinsic ambipolarity. This is not easy to reconcile with steady solutions we have calculated by the Monte Carlo method from exponential decay of the distribution function, which have unequal sources of ions and electrons (1). The TRAN code estimates how long an ion or an electron stays in the plasma. Tracking orbits of the electrons is what reveals the importance of quasineutrality in these calculations of transport. Because the mean free path of the particles is orders of magnitude longer than the circumference of the torus, it is hard to see how conservation of momentum could make the numbers of ions and electrons agree to several significant digits in a real volume of plasma.

Coils

When running equilibrium, stability, and transport codes to design stellarator experiments, one usually starts by looking for Fourier coefficients specifying the boundary of a plasma with good physical properties (7). After that it becomes necessary to use a line-tracing code such as NESCOIL to find a set of coils that generate a corresponding external magnetic field according to the Biot–Savart law (17, 18). This is not a well posed mathematical problem, because the coils may not satisfy physically reasonable constraints unless care is exercised in choosing the coefficients that define the original equilibrium. It is advisable to limit the degree of the harmonics that are used and to avoid shapes for the separatrix that are unnecessarily complicated so that the final coils will furnish a robust equilibrium that continues to have good Kolmogorov–Arnold–Moser surfaces after realistic changes are made in the geometrical parameters of the configuration. The coils become excessively twisted within the control surface where they are wound if that is located too far from the plasma. This problem is as important as physics issues that sometimes receive more attention.

Unwanted islands are suppressed by manipulating resonant harmonics in the Biot–Savart law. This seems to succeed if the term in the magnetic spectrum with the same resonance as an undesirable island turns out to be small in the original equilibrium. Such islands can be controlled by auxiliary shaping coils or by modification of the larger modular coils that generate the external magnetic field, which should be only moderately twisted. Numerical calculations show that the quality of the Kolmogorov–Arnold–Moser surfaces may improve if the coordinates defining the coils are smoothed. One can also address these issues for finite β by running free-boundary codes (19).

Design

The modular Helias-like Heliac 2 (MHH2) is a compact stellarator designed by computational methods to perform well at reactor conditions (1, 14). Sixteen moderately twisted modular coils define an external magnetic field with rotational transform near 2/5 that has robust magnetic surfaces. There are two field periods, and the aspect ratio of the plasma is 3.5. The nominal width of a square cross section of the coil structure is a quarter of the plasma radius. Quasiaxial symmetry of the magnetic spectrum provides excellent confinement at low collision frequencies. The shape of the separatrix produces a good magnetic well and assures nonlinear stability up to an average β of 4% for ideal MHD modes (14). According to our transport theory, an MHH2 experiment could be constructed with a major radius of 2 m that might achieve a peak ion temperature of 10 keV in a magnetic field of 2 T at an electron density of 1.5 × 10¹³ cm^–3. The configuration can be operated effectively as a hybrid and can be conceived of as an advanced tokamak.

It is a challenge to predict the performance of future experiments by making numerical calculations. With that in mind, we examined the National Compact Stellarator Experiment proposal for a proof of principle stellarator experiment at the Princeton Plasma Physics Laboratory, which will have a plasma of aspect ratio 4.5 with three field periods. To get adequate space between the separatrix and the first wall, the filaments defining the modular coils become kinky, which may lead to deterioration of the magnetic surfaces in the edge region. Consequently, we ran the NSTAB and TRAN computer codes for an example with less complicated geometry. Nonlinear stability of the alternate configuration is favorable, and the β limit exceeds 4%, but islands may appear when the rotational transform ι crosses 3/7, 3/6, and 3/5 near the separatrix.

Free-Boundary Code

We describe next a free-boundary algorithm for the design and analysis of stellarator and tokamak equilibria. An important application is the reconstruction of plasma geometry from available diagnostics and observations. Correct mathematical formulations of the problem are required not just to meet acceptable standards of self-consistent physics but also to arrive at computer simulations of nonlinear phenomena that can be correlated with experimental observations to make predictions going beyond those provided by standard scaling laws (20).

Some of the mathematical difficulties in understanding plasma confinement in three dimensions have been explained already. A reliable free-boundary code is needed for the engineering analysis of coil sets, and a satisfactory treatment of islands and current sheets at resonant surfaces in discontinuous solutions of the equilibrium problem is required. We present a mathematical method to attack these questions by solving the free-boundary value problem numerically.

The NSTAB code applies the MHD variational principle to calculate weak solutions of the equilibrium equations that satisfy the relations for the magnetic field B in a plasma with fixed boundary defined by the formula where r, v, and z are cylindrical coordinates and u is a poloidal angle (8). In this theory the functions s, θ, and ϕ comprise an invariant system of flux coordinates and ζ is an additional Clebsch potential that specifies the distribution of current on the family of nested toroidal surfaces s = constant. The strength of the field in the plasma has a Fourier series representation the coefficients B_mn = B_mn(s) of which comprise the magnetic spectrum. Elementary transformation of the equilibrium equations shows that ζ has a related representation with small denominators that vanish at resonant surfaces corresponding to rational values of the rotational transform ι. In this formulation, ∇s ×∇ζ is the current density (16). The singularities occurring at rational surfaces explain why it is important to work with weak solutions of the problem.

On an outer control surface, let ζ₀ stand for a distribution of current with level curves that are filaments for a system of modular coils that generate the external magnetic field controlling a stellarator equilibrium, and let ζ₁,..., ζ_k be corresponding values of the Clebsch potential ζ that describe distributions of surface current on a discrete set of magnetic surfaces s = s_j inside the plasma. The NSTAB code provides an accurate calculation of the magnetic spectrum B_mn when the shape factors Δ_mn defining the separatrix are known, and therefore it also provides a way to compute ζ_j. Moreover, the Biot–Savart law gives an approximation of the magnetic field that is valid both inside and outside the plasma, where the double integrals are evaluated over the control surface and the discrete flux surfaces s = s_j, and δs_j is an increment in s_j.

To solve the free-boundary problem, what is required is to adjust the Fourier coefficients Δ_mn such that the shape of the plasma becomes compatible with magnetic lines in the solution. In a preparatory code implementing this concept, which has features in common with NESCOIL, the surface integrals from the Biot–Savart law are approximated by line integrals over filaments modeling the Pfirsch–Schlueter current and the bootstrap current (cf. ref. 18).

In practice the method exploits line tracing to recalculate the shape of the separatrix and eliminate islands from the solution. This problem is not well posed because of the Kolmogorov–Arnold–Moser theory, which means that only an approximate answer should be sought (5). However, the approach is effective for the design and analysis of quasisymmetric stellarators if just a few terms in the Fourier series for the Clebsch potential ζ and a sparse set of flux surfaces s = s_j suffice to give reliable results. Fig. 5 displays a Poincaré section of the magnetic lines in a computation for the MHH2 stellarator by a version of the method that uses just one internal flux surface s = s ₁ in the Biot–Savart law and one coefficient B ₁₀ in the spectral representation of ζ. In the reconstruction of the free boundary, the location of the magnetic axis determined by the NSTAB code plays a significant role. Finally, let us observe that the form of the Biot–Savart law we have presented can also be applied more directly to find modular coils for a prescribed equilibrium configuration when β is positive.

Reactor Studies

A reactor study of compact stellarators is currently being undertaken by the US Department of Energy. Principal configurations under consideration as candidates for the study are the National Compact Stellarator Experiment with three field periods and the MHH2 with two field periods. Here we review results from runs of our computer codes that were performed in connection with this project.

As we have seen, NSTAB calculations show that the LHD stellarator is linearly unstable but remains nonlinearly stable at the β of 3.2% that has been achieved experimentally. At a peak ion temperature of 3 keV observed in the LHD, the TRAN code predicts an energy confinement time of 160 ms that agrees with the measured value (1). This validates our model of anomalous transport based on considerations of quasineutrality. On the other hand, computations of the energy confinement time τ_E in an LHD reactor study call for major radius 25 m and plasma radius 4.5 m at average T = 12 keV, n = 1.2 × 10¹⁴ cm^–3, and B = 5 T. That is consistent with the low ion temperatures achieved thus far for a hydrogen plasma in the experiment (2).

Good correlation of our computations with observations in the LHD has been applied to assess equilibrium, stability, and transport for a compact MHH2 stellarator. According to our transport theory, an MHH2 reactor could be constructed with a major radius of 7.5 m and a plasma radius of 2.2 m that would achieve average temperature T = 16 keV in a magnetic field of 5 T at average electron density n = 10¹⁴ cm^–3. The energy confinement time is >2 s if the drop in the radial electric potential is twice the temperature. Dependence of ion, electron, and α-particle transport on the radial electric field and the mirror coefficients Δ₀₁ and B ₀₁ has been examined, and transport can be improved by using the mirror term to raise the radial electric field.

The β limit of the MHH2 power plant is 4%, and 16 only moderately twisted coils provide robust magnetic surfaces. Above the β limit a ballooning mode appears in the equilibrium that is wall-stabilized and can be detected by long runs of the NSTAB code (compare with Fig. 6). The gap between the separatrix and central filaments defining the coils exceeds 140 cm everywhere, and there is room for a maintenance port between any pair of adjacent coils. This fusion reactor is designed to deliver 1 GW of electric power.

A principal issue for stellarators is the quality of the magnetic surfaces, especially near the edge of the plasma. For compact stellarators the rotational transform per field period is usually large, so low-order resonances may become dangerous. In the MHH2 reactor the rational surface ι = 2/5 appears in the scrape-off layer at low values of β, and the rational surface ι = 2/4 may emerge near the magnetic axis when β becomes larger so that there is significant bootstrap current. The resonant terms B ₄₂ and B ₅₂ in the magnetic spectrum are both relatively small, which means that the widths of the corresponding islands are not expected to be large (1).

We have described results that illustrate the role of mathematical and computational methods in magnetic fusion. To complete the research it would be desirable to construct an experiment for a quasiaxially symmetric stellarator such as the MHH2 with two field periods and only moderately twisted coils.