Axisymmetric simulations of vertical displacement events in tokamaks: A benchmark of M3D-C¹, NIMROD, and JOREK

A vertical displacement event (VDE) denotes the vertical movement of a tokamak plasma toward the vessel walls, generally leading to a complete loss of the plasma confinement. In the case of a cold VDE, the loss of the control of the plasma position is caused by a rapid change in the plasma pressure and current density profile due to a thermal quench. In a hot VDE, the plasma initially maintains most of its thermal energy and the loss of control occurs for different reasons, for example, when a stability threshold in elongation is exceeded. The vessel wall currents produced by the VDE can lead to large transient forces on the vessel, including sideway forces that can challenge the mechanical supports of the vessel. Also, these forces have the potential to create significant mechanical stresses in the vacuum vessel, especially in tokamaks with large magnetic fields and large plasma currents. Therefore, comprehensive analysis, including large-scale simulations seeking to identify the worst case VDEs, is necessary as part of the design process.

Although a hot VDE is initially an axisymmetric instability, 3D instabilities develop during the course of the event. In particular, when the edge safety factor drops below some critical value due to scraping-off by the wall or by impurity cooling of the edge, the plasma becomes unstable to an external kink or a resistive wall mode (RWM). These asymmetries can lead to asymmetric forces on the vessel walls. Since asymmetric forces will lead to asymmetric stresses and might even rotate in mechanical resonance with the vessel structures,¹ they can lead to large local stresses in the vacuum vessel.

Employing 3D nonlinear magnetohydrodynamic (MHD) codes to understand and predict the consequences of different types of VDEs for tokamak operation has become an active field of study. Recently, 3D VDE simulations have been performed by Strauss^2,3 with the M3D code, Pfefferlé et al.⁴ with M3D-C¹, Artola⁵ with JOREK, and Sovinec and Bunkers⁶ using NIMROD.

Combining a global 3D MHD model for the plasma evolution with implicit time stepping and a model for a resistive wall, NIMROD,⁷ JOREK,^8,9 and M3D-C¹ (Refs. 10 and 11) are among a small set of codes possessing the necessary capabilities for 3D VDE simulations. Nevertheless, a benchmark between any of the three codes involving VDE calculations had not been performed, and code-to-code comparison is an essential tool in this context. As discussed in Secs. II and IV–VI the codes use significantly different numerical models, in particular for the resistive wall (discretizing vacuum vs Green's function method). Except for idealized problems, there are no analytic solutions to verify the codes' correctness, and experimental data are not available in enough detail to perform true validation. Code-to-code comparison is thus the only way to verify non-trivial calculations in detail.

In the following, we present the setup and results of a benchmark exercise between M3D-C¹, NIMROD, and JOREK, which is based on a vertically unstable NSTX equilibrium. Note that although an experimental equilibrium is used, this work is solely intended to be an inter-code benchmark exercise with the purpose of code verification focused on VDE relevant physics. An experimental validation of the physics models implemented is a separate but equally important and ongoing endeavor.

The goal of this work is to provide the fusion community with a useful set of validated standard benchmark cases to be used to test MHD codes with a resistive wall model, which are to be applied to study VDEs. The benchmark calculations presented in this paper are strictly 2D (axisymmetric) even though the three codes involved are fully 3D. These 2D benchmark simulations represent a crucial first step toward a realistic 3D nonlinear VDE benchmark, which is ongoing.

The 2D VDE problem discussed in this work could, of course, also be solved by evolutionary equilibrium codes such as DINA.¹² However, this does not make it less important to validate 3D nonlinear MHD codes on this problem, in particular since the formalism and numerical representations that these use are quite different from those in DINA. Moreover, DINA makes some assumptions regarding the halo current and the neglect of the plasma inertia, which are not made in the calculations presented here. Our results show that, for the 2D problem, these assumptions are justified.

In Sec. II, a brief description of the M3D-C¹ model used for VDE calculations is given. In addition, it is discussed how response currents in the open field line region can cause a deviation from the linear dependence of the VDE growth rate on the resistivity of the vessel wall. The setup of the benchmark case is described in Sec. III. The differences between the three codes are discussed, and the results of the benchmark are presented in Secs. IV–VI, where the three sections are concerned with a linear benchmark, a benchmark of the linear phase of nonlinear simulations, and a benchmark of the axisymmetric, nonlinear evolution of a full VDE, respectively. Good agreement is found between the linear VDE growth rates and between the nonlinear evolution of the VDE calculated with the three codes. In particular, the evolution of currents in the plasma and in the vessel wall is compared. In Sec. VI, we also analyze why the reduced MHD model used for the JOREK simulations is able to reproduce the results of the full MHD models used in the other two codes. A summary and outlook on future work are given in Sec. VII.

The M3D-C¹ code is a high-order finite element code that solves the nonlinear time-dependent extended MHD equations. It uses a split-implicit time advance in order to enable simulations over transport time scales. For the spatial discretization, triangular wedge finite elements are used.¹³ 

For simulations of VDEs, a three region model (as illustrated in Fig. 1) is used: Within the central region, the nonlinear extended MHD equations as described in Ref. 10 are solved. It is enclosed by a resistive wall region of arbitrary thickness. Between the resistive wall and the outer domain boundary, there is a vacuum region. The mesh resolution can be locally increased where needed, e.g., in the vicinity of the resistive wall.

At the boundary between the plasma domain and the resistive wall domain, no-slip boundary conditions are employed, and the temperature and density are kept at fixed values. Inside the resistive wall domain, the magnetic field is evolved according to $∂tB=−∇×(ηwallj)$⁠, where $ηwall$ is the wall resistivity and the current density is given by $j=(∇×B)/μ0$⁠. In the vacuum region, $j=0$⁠. There are no boundary conditions on the magnetic field at the resistive wall. Halo currents can flow into and out of the wall.

Ideal wall boundary conditions are used at the outer domain boundary. The sensitivity of the VDE growth rates to the location of the ideal wall has been tested, and the influence of the ideal wall has been estimated to result in a deviation below 10% relative to no-outer-wall conditions.

The resistive wall model has been successfully benchmarked against analytic solutions for RWMs.¹¹ The model has recently been extended to provide an option for a spatially varying resistivity, which can be used to model conducting and non-conducting structures around the plasma. Furthermore, impurity and pellet models have been recently added to the plasma model for disruption mitigation studies.¹⁴ These options however have not been used for the calculations presented here.

M3D-C¹ can be used for 2D and 3D linear and for 2D axisymmetric nonlinear and 3D nonlinear simulations. 2D nonlinear simulations can also be restarted in both linear and 3D nonlinear modes.

In some axisymmetric codes that are used for VDE calculations such as DINA¹² and TSC,¹⁵ the width of the halo region is an input parameter. Note that this is not the case in M3D-C¹, NIMROD, and JOREK, where the width of the halo region is not artificially set, but it is determined by the heat transport model and can be adjusted via heat diffusion anisotropy. Figure 2 shows the resulting temperature profiles on the midplane for 2D nonlinear VDE simulations performed with M3D-C¹ with different values of the heat diffusion anisotropy, i.e., the ratio of the parallel heat diffusion coefficient $κ∥$ to the perpendicular heat diffusion coefficient $κ⊥$⁠. (Here, the value of $κ∥$ has been varied while keeping the value of $κ⊥$ fixed.) Since response currents in the halo region slow down the VDE, smaller values of the heat diffusion anisotropy leading to higher edge temperatures cause reduced VDE growth rates (γ). Note that a detailed study in which the thermal conductivity in 2D nonlinear M3D-C¹ simulations is used as a parameter to scan a wide range of cases and identify worst case scenarios for ITER is described in Ref. 16. A separate study examines the effects of plasma thermal conduction modeling and boundary conditions in axisymmetric VDE evolution.¹⁷ 

The temperatures in the open field line region, affected by both the edge temperature boundary condition and by the ratio of thermal conductivities as shown in Fig. 2, can play an important role in VDE simulations. In the initial phase of hot VDEs, the vertical movement of the plasma is inhibited by response currents in the conducting vessel structures. Due to the finite resistivity of these structures, the response currents decay, and thus, the linear VDE growth rate is determined by the resistive time of the vessel in this initial phase. However, if the electron temperatures in the open field line region are large in a simulation, its resistance can become sufficiently small to compete with the vessel resistance, leading to response currents forming in the open field line region.

Figure 3(a) shows the linear VDE growth rates for different values of the wall resistivity and different values of the edge temperature in linear M3D-C¹ simulations. By edge temperature, we mean the value of the electron temperature at the boundary between the plasma domain and the resistive wall domain. This value is changed by changing the value of the equilibrium edge pressure while keeping the same edge density. Note that details on the setup and parameters of these simulations are given in Sec. III and Table I(a).

In the limit of small wall resistivities and small edge temperatures for which the edge plasma resistivity is very large, the expected linear dependence of the growth rate on $ηwall$ is recovered. For larger edge temperatures (i.e., lower $ηedge$⁠) or larger values of $ηwall$⁠, the resistance of the wall becomes comparable to the resistance of the open field line region. In this regime, the VDE growth is slowed down by response currents in the open field line region as illustrated in Figs. 3(b) and 3(c). In the limit where the wall resistance is much larger than the resistance of the open field line region, the VDE growth rate is only determined by the plasma edge resistivity and becomes independent of the wall resistivity.

The equilibrium used for this benchmark The equilibrium used for this benchmark case has been reconstructed with the EFIT²² code from the NSTX discharge #139536 at t = 309ms. It is illustrated in Fig. 4. Instead of using the complicated shape of the NSTX vacuum vessel, an axisymmetric rectangular resistive wall is used to simplify the geometry. The corners of the inner boundary of the resistive wall domain are at (⁠$R=0.24 m, Z=±1.4 m$⁠) and (⁠$R=1.6 m, Z=±1.4 m$⁠). The thickness of the resistive wall is set to $Δw=0.015 m$⁠.

The equilibrium position of the magnetic axis is (⁠$Raxis=1.07 m,Zaxis=−0.015 m$⁠). The toroidal magnetic field on axis is $Btor=0.37 T$⁠, and the total toroidal plasma current is $Itot=5.7×105 A$⁠. The difference between the poloidal magnetic flux at the boundary and at the magnetic axis is $Ψbnd−Ψaxis=−0.059 V s$⁠, where $Ψ=−∫Bpol dA/2π$⁠. The temperature profile is given by $Te(Ψ)=1 keV·(p(Ψ)/paxis)0.6$⁠. The pressure and current density profiles are defined in the geqdsk equilibrium file. Note that the geqdsk file and files containing the coil positions and currents are available in the supplementary material.

Dynamic viscosity, perpendicular and parallel heat diffusion coefficients, and the particle diffusion coefficient are constant in space and time. Their values are given in Table I(a) for the simulations discussed in Sec. II and Table I(b) for those discussed in Secs. IV–VI. The isotropic plasma resistivity is given by the (perpendicular) Spitzer resistivity {i.e., $η(Te)=1.03×10−4·Z· ln Λ·(Te[eV])−3/2 Ω m$⁠, where Z = 1 and $ln Λ=17$}. The ion mass is set to twice the proton mass. A loop voltage is not applied.

While NIMROD and M3D-C¹ have similar physics models, the numerical methods differ, which makes benchmarks between these two codes particularly valuable. In contrast to M3D-C¹, NIMROD uses high-order C⁰ quadrilateral finite elements in the R-Z plane and a Fourier spectral representation for the toroidal discretization. When testing resolution for these computations, bicubic and biquartic elements have been applied. As described in Ref. 6, the NIMROD computations presented here use a thin-wall model that couples dynamics in the plasma subdomain with numerically computed magnetic responses in an outer vacuum subdomain. This differs from M3D-C¹, which represents a thick wall within a single mesh that spans all domains.

The geometry of the outer subdomain used in the NIMROD benchmark computations matches the shape shown in Fig. 1, except that the top and bottom corners at $R=0.02 m$ are not rounded. The computations also used fixed viscosity and thermal conductivity coefficients (same as M3D-C¹) without the dependence on plasma density that is shown in Refs. 6 and 7.

Figure 5 shows a comparison of the linear VDE growth rates obtained from linear M3D-C¹ and NIMROD simulations [using the parameters listed in Table I(b)]. The growth rates agree well over a wide range of wall resistivities. The largest deviation between the growth rates is 13%, and it occurs at the largest values of the wall resistivity, where the results are most sensitive to the representation of the equilibrium and halo responses.

For simulations with JOREK, which include a resistive wall, the JOREK-STARWALL coupling is used.^9,18 Similar to NIMROD and in contrast to M3D-C¹, JOREK uses a spectral representation for the toroidal discretization. Cubic Bézier finite elements are used for the discretization in the R-Z plane. There are a few differences between the model that JOREK uses for the benchmark simulations and the models that M3D-C¹ and NIMROD use: (i) Although JOREK has a full MHD model, it uses a reduced MHD model¹⁹ for the VDE calculations presented here since the JOREK-STARWALL coupling is not yet available for the full MHD model. (ii) In JOREK-STARWALL, the vacuum contribution is implemented by using Green's function method. Therefore, it is not necessary to discretize the vacuum region and apply ideal wall boundary conditions in an outer boundary. This property comes from the fact that the full vacuum response can be expressed as a function of the magnetic field at the plasma boundary. (iii) At the resistive wall, instead of no-slip boundary conditions, only the normal component of the velocity vanishes.

For the JOREK simulations presented here, a polar grid is used with increased resolution in the region surrounding the point of contact between the plasma and the wall. The number of Bézier elements used is 22 000, and the number of linear triangular elements used for the representation of the wall is 48000.

Since JOREK does not have an option that allows linear simulations with toroidal mode number n = 0, we compare the VDE growth rates in the early, linear phase of the evolution obtained in 2D axisymmetric nonlinear simulations.

In order to be able to run benchmark cases in the regime where the VDE growth rate is not influenced by response currents in the open field line region, the value of the edge temperature has to be sufficiently small. Since in nonlinear simulations too small values of the edge temperature can lead to numerical problems, we use a small temperature offset only within the calculation of the Spitzer resistivity such that $η(Te)=ηSpitzer(Te−Te,off)$ in all three codes. Here, the edge temperature is $Te,edge=14.65 eV$ and the offset is $Te,off=13.65 eV$⁠, which results in an effective edge resistivity corresponding to a temperature of $Te,eff=1 eV$⁠. For simplicity, the Ohmic heating term in the temperature equation is switched off.

Figure 6 shows a comparison of the resulting VDE growth rates. They have been obtained from the 2D nonlinear simulations by fitting $Zaxis=a+b· exp (γt)$ to the time trace of the vertical position of the magnetic axis, where γ is the growth rate. Only the early, linear phase of the evolution (vertical position of the magnetic axis between $Zaxis=−1.64 cm$ and $Zaxis=−3.04 cm$⁠) has been taken into account.

All three codes find the expected linear relationship between the VDE growth rate and wall resistivity, and the results agree well. The deviation between the obtained growth rates is around 3% or less for most wall resistivities and does not exceed 12% in the other cases, except for a deviation of 15% between the M3D-C¹ and the NIMROD result for the smallest wall resistivity. We also show the growth rates obtained from linear M3D-C¹ simulations, and they agree well with the results of the nonlinear calculations. However, it should be noted that such good agreement between linear M3D-C¹ results and the nonlinear results is only achieved if the linear M3D-C¹ simulation is not directly started from the original equilibrium. First, the simulation is run in nonlinear mode for a few time steps (in this case until the plasma has been drifted by $∼2 mm$⁠), and then, it is restarted as a linear simulation. There are several possible reasons for this to be necessary. First, we find that the geqdsk equilibrium has a layer of reversed current density that increases in magnitude when approaching the inboard side of the separatrix. When the equilibrium is re-solved on the meshes of the respective codes, the treatment of the sharp termination of current density at the separatrix involves mesh-scale numerics, leading to discrepancies among the codes. A second possibility is that this might be due to a slight difference between the initial, ideal equilibrium and the stationary state the system goes into in the presence of non-ideal terms. In either case, the small difference relaxes quickly when run in nonlinear mode, i.e., when the transport equations are applied to the equilibrium in the time advance.

In the following, the results obtained by JOREK, NIMROD, and M3D-C¹ on the further axisymmetric nonlinear evolution of a VDE are compared. The setup and parameters of these simulations are the same as for the simulations discussed in Sec. V, except that $Te,off=0$ such that the edge resistivity corresponds to an edge electron temperature of $Te,edge=14.65 eV$⁠. The resistivity of the wall has been set to $ηw=3×10−6 Ω m$⁠. While it should be emphasized that a comparison with experimental measurement results is out of the scope of this work, we note that VDE growth times in NSTX are tens of milliseconds⁴ and the chosen wall resistivity leads to a growth time of the same order of magnitude.

In addition, a thermal quench has been artificially initiated during the course of the evolution. In 3D nonlinear MHD simulations, e.g., Ref. 4, the decrease in the edge safety factor during the course of a VDE causes non-axisymmetric instabilities to develop. These 3D instabilities cause the magnetic flux surfaces to break up, which leads to greatly increased thermal transport. Since this effect cannot occur in axisymmetric simulations, an artificial thermal quench is initiated by increasing the perpendicular heat diffusion coefficient by a factor of 500 when the plasma becomes limited by the wall. Also, the particle diffusion coefficient is multiplied by a factor of 20. The coefficients for the post-thermal quench phase are chosen to help match experimental disruption times in NSTX (⁠$∼10 ms$⁠).

The poloidal magnetic flux at this time when the plasma becomes limited by the wall in the M3D-C¹ simulation and the time traces of the thermal energy in the M3D-C¹, JOREK, and NIMROD simulation are shown in Fig. 7.

In order to enable a meaningful comparison of the results, the signals are slightly shifted in time such that the points in time of the first plasma-wall contact, e.g., the start of the thermal quench, coincide. This compensates for differences caused by the exponential dependency on the initial conditions. (The plasma first touches the wall at $t≈126 ms$ in the JOREK simulation, at $t≈87.4 ms$ in the NIMROD simulation, and at $t≈91.5 ms$ in the M3D-C¹ simulation.) Figure 8 compares the time traces of the vertical and radial positions of the magnetic axis, the toroidal current enclosed by the last closed flux surface (LCFS), the total toroidal current inside the LCFS and in the open field line region, and the net toroidal current in the resistive wall. [The NIMROD time traces in Figs. 8(d) and 8(e) are missing because NIMROD does not currently have the corresponding diagnostics.]

In addition, the halo current at the plasma-wall interface, i.e., the component of the current density perpendicular to the wall at the wall, is shown for a point in time during the late evolution (when $Zaxis≈−1.23 m$⁠). The halo current is plotted against the distance along the wall, measured counterclockwise, starting at the low-field side midplane. For the JOREK simulation, the halo current is calculated from $j×B=∇p$⁠, assuming that the plasma is in equilibrium. Note that the location of the halo current spikes resulting from the M3D-C¹ simulation appears to be slightly shifted with respect to the other two traces. This is an artifact caused by the M3D-C¹ resistive wall having a slightly larger circumference since its corners are less rounded and then the ones of the resistive walls used for the JOREK and NIMROD simulations.

As expected, the halo current flows into and out of the wall in a narrow region surrounding the contact point of the last closed flux surface and the wall. Despite the differences in physics models and numerical implementation between the three codes, the results agree well. This implies that the reduced MHD model that JOREK uses reproduces the results of the full MHD models in M3D-C¹ and NIMROD.

The good agreement between the reduced MHD model and the full MHD models originates from the presence of a large vacuum toroidal field.²⁰ The formulation of the energy principle for n = 0 reveals that important stabilizing terms involving the large toroidal magnetic field (⁠$Bϕ$⁠) can be minimized by choosing the following form of the plasma velocity v (or plasma displacement $ξ$⁠):

where γ is the growth rate, $ϕ$ is the toroidal angle, and u is the velocity stream function. This condition is called the “slip motion condition,”²¹ in which the plasma moves across the large toroidal field without doing work against it.

The velocity representation of the slip motion condition is the main assumption used in reduced MHD, which is equivalently justified whenever $F≡RBϕ≈Fvacuum$ and explains the excellent agreement with the full MHD models. The nature of the problem makes the full velocity to be well represented by the reduced MHD velocity. In the presented simulations, the total toroidal field differed from the vacuum toroidal field by a maximum of 5% despite the small aspect ratio of the configuration.

A VDE benchmark case for nonlinear MHD codes has been presented. It is based on a vertically unstable NSTX equilibrium and uses an axisymmetric rectangular model for the resistive wall. The 3D nonlinear MHD codes M3D-C¹, NIMROD, and JOREK are applied to the benchmark case, but they run in a linear and a 2D (axisymmetric) nonlinear mode. Linear simulations show that the expected linear dependence of the VDE growth rate on the resistivity of the vessel wall (⁠$ηw$⁠) is recovered for small values of $ηw$ and a sufficiently large resistivity in the open field line region. If the temperature in the open field line region becomes too large, response currents build in this region and slow down the VDE growth.

Agreement within approximately 10% is found between the results of the linear benchmark simulations obtained by NIMROD and M3D-C¹, as well as between the VDE growth rates in the linear phase of axisymmetric nonlinear simulations performed with NIMROD, M3D-C¹, and JOREK. Where the plasma response is most important, the agreement of growth rates based on fits from the early nonlinear computations (Fig. 6) is better than the agreement of growth rates from the true linear computations (Fig. 5). We checked a number of possible sources for the discrepancy in the linear computations, including normalizations, outer-wall geometry, thermal conductivity modeling, and electrical resistivity. By the process of elimination and by noting the improved agreement after short nonlinear evolution, we infer that a likely cause is the sharp layer of edge current density in the equilibrium. The equilibrium is not smoothed in the true linear computations, and the NIMROD and M3D-C¹ representations differ. This will affect the linear computations, particularly when meshes are not aligned with flux surfaces, and such alignment is not suited to the nonlinear VDE computations that are the focus of this benchmarking.

The further axisymmetric nonlinear evolution of a selected case has been calculated using the three codes, and the time traces of the position of the magnetic axis, the toroidal plasma current, and the toroidal current in the resistive wall, as well as the resulting halo current, have been compared. Despite the differences in physics models and numerical methods, the results agree well. This implies that for the benchmark cases presented here, the reduced MHD model that JOREK employs reproduces the results of the full MHD models of NIMROD and M3D-C¹. In particular, it is shown that in 2D calculations, the reduced MHD model can calculate the halo currents correctly, which has been an open question in the past. Whether the reduced MHD model is still reproducing the results of the full MHD model in 3D simulations will need to be investigated.

The linear and 2D nonlinear benchmarks described in this work are a necessary step to form a solid basis and serve as a starting point for a 3D benchmark. They demonstrate that the codes maintain force balance, where profiles are determined by rudimentary transport models, over times that are far longer than Alfvénic propagation times. In addition, features relevant for VDEs that are already important in 2D, such as different resistive wall models and the use of reduced MHD for the calculation of halo currents, have been tested. Future work will include the extension of this benchmark toward fully 3D nonlinear simulations with the three codes. Furthermore, we hope that other linear, axisymmetric nonlinear, or 3D nonlinear codes used for VDE calculations might be applied to the benchmark case presented here as well.

See the supplementary material for the geqdsk file defining the equilibrium and files containing the coil positions and currents for the discussed benchmark cases.

This work was supported by U.S. DOE Contract Nos. DE-AC02-09CH11466 and DE-SC0018001 and the SciDAC Center for Tokamak Transient Simulations. Some of the computations presented in this article used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231. The support from the EUROfusion Researcher Fellowship programme under Task Agreement No. WP19-20-ERG-DIFFER/Krebs is gratefully acknowledged. Part of this work was carried out within the framework of the EUROfusion Consortium and received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

ITER is the Nuclear Facility, INB No. 174. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization. This publication is provided for scientific purposes only. Its content should not be considered as commitments from the ITER Organization as a nuclear operator in the frame of the licensing process.