Comparison between the analytic theory of n = 0 vertical displacement modes in magnetically confined plasmas of fusion interest and numerical simulations using the extended-MHD code NIMROD is presented. Agreement between analytic and numerical results is highly satisfactory. Differences are interpreted to be caused mostly by the different wall shape and by the presence of a halo plasma surrounding the hot plasma adopted in NIMROD. A numerical study of vertical displacement oscillatory modes [Barberis et al., J. Plasma Phys. 88, 905880511 (2022)] is presented. Axisymmetric X-point currents supported by the halo plasma are discussed. The article provides a successful benchmark and a useful starting point for future numerical investigations of n = 0 modes using more realistic tokamak geometry and plasma equilibria.

Vertical, axisymmetric modes with toroidal mode number n = 0 are of concern for the safe operation of tokamak plasma discharges, as these modes, if growing uncontrolled, may lead to vertical displacement events (VDE) and plasma current disruptions. Therefore, a considerable amount of work has been dedicated to the study of this problem1–8 (references quoted above are a non-exhaustive sample). In general, active feedback control systems have been developed, which under normal circumstances can prevent the occurrence of VDE.9–14 

Recently, there has been renewed interest in the theory of n = 0 modes, prompted by the observation of saturated, high frequency n = 0 oscillations in the JET tokamak.15,16 Analytic theory work, which was possible on the basis of simplified plasma equilibria and geometry, indicated that magnetic X-points, which are intrinsic in modern tokamak plasmas with divertor configurations, can have a significant impact on the dynamics of these modes.17,18 In addition, a new type of n = 0 mode, dubbed19 vertical displacement oscillatory mode (VDOM), which could be driven unstable by resonant interaction with fast ion orbits,20 was proposed. Clearly, the treatment of more realistic equilibria requires numerical work. It must be pointed out that treating X-point effects on n = 0 perturbations is challenging from the numerical point of view, since the toroidal field line going through a magnetic X-point is resonant, in the sense that B eq · k = 0 for n = 0 at magnetic X-points, with k being the perturbation wavevector and B e q being the equilibrium magnetic field.

The purpose of this article is to provide a numerical benchmark between analytic theory and linear, initial value simulations performed by means of the extended-MHD code NIMROD.21 Some preliminary tests of the numerical model for n = 0 vertical modes with NIMROD have already been presented in Ref. 22, with this manuscript providing a more complete and more detailed benchmark. The benchmark will include scenarios where n = 0 oscillations are either unstable according to ideal-MHD, or executing stable oscillations with a frequency close to the poloidal Alfvén frequency. In these numerical simulations, we will mostly stay away from the difficult numerical treatment of X-point effects, in the sense that the hot plasma density is assumed to vanish before the magnetic separatrix is actually reached. Nevertheless, evidence of X-point effects can be seen in the simulations presented in this article, as the X-points actually lie in a low-density, high-resistivity halo plasma, rather than in vacuum. Thus, X-point effects manifest themselves through the appearance of X-point currents carried by the halo plasma.

In order to better clarify the scope of the present article, we point out that four different scenarios come to mind for the occurrence of n = 0 modes in tokamak plasmas. These modes can appear as follows: (1) Unstable vertical displacement perturbations, potentially leading to VDE;1–8 (2) Stable, vertical displacement oscillatory modes (VDOMs),14,19 with a frequency of the order of the poloidal Alfvén frequency (a few hundred kHz in a tokamak such as JET), which could be driven unstable by fast ions;20 (3) Global Alfvén eigenmodes23 (GAE), also with a frequency of the order of the poloidal Alfvén frequency, typically somewhat higher than the VDOM frequency, depending on details of the q profile; and (4) Geodesic acoustic modes (GAM),24,25 with a relatively low frequency of the order of the sound wave frequency (e.g., a few tens of kHz in a tokamak plasma such as JET). GAE and GAM may also be driven unstable by fast ions.15,26 In this article, we will be dealing only with the first two types of n = 0 perturbations, i.e., unstable vertical displacements and VDOMs. Also, fast ions physics, which could lead to VDOM destabilization, is not included in this article.

This article is organized as follows. In Sec. II, we briefly review analytic theory results on unstable vertical displacements and VDOMs that are up for comparison with the numerical simulations presented in Sec. IV. In Sec. III, we recall the main features of the NIMROD code. Differences between the assumptions and approximations used in the analytic model and the physics underlining the NIMROD code are highlighted. The straight tokamak equilibrium adopted in analytic theory and in NIMROD simulations is presented. In all linear NIMROD simulations discussed in this article, the confining wall is taken as ideally conducting. Results are presented in Sec. IV, subdivided into five subsections. Subsection IV A considers the no-wall limit, where vertical displacements are ideal-MHD unstable, with a growth rate that depends on plasma ellipticity. In Subsection IV B, we present wall position scans at fixed values of plasma ellipticity. As the wall is moved closer to the plasma, the ideal-MHD vertical instability is suppressed (passive wall stabilization) and the n = 0 perturbation transitions into stable VDOMs. Subsection IV C presents the results obtained by scanning plasma elongation (plasma triangularity is not considered in this article) at fixed wall position, showing again the transition from ideally unstable vertical displacements to stable VDOMs. Subsection IV D focuses on oscillatory modes. Subsection IV E examines the question of X-point currents carried by the halo plasma. Conclusions are presented in Sec. V.

In this section, we briefly summarize the analytic results on vertical modes and VDOMs obtained in Refs. 17–19, 22, and 27, limiting our considerations to those results that are relevant for the numerical benchmarking presented in this article. The magnetic field is represented as B = e ϕ × ψ + B ϕ e ϕ, where e ϕ is the unit vector along the ignorable toroidal direction, and B ϕ is constant. The plasma flow is represented as v = e ϕ × φ + v ϕ e ϕ. Neglecting toroidal curvature effects, an idealized, low-beta, straight tokamak equilibrium is assumed, with an equilibrium, toroidal current density of the type J e q = J 0 H ( ψ b ψ e q ), where H(x) is the Heaviside (unit step) function, ψ e q ( x , y ) is the equilibrium magnetic flux function, and ψ e q = ψ b = const is an elliptical flux surface with minor semiaxis a and major semiaxis b representing the plasma boundary. Beyond the plasma boundary, the density also drops to zero, and, therefore, in analytic work, the magnetic X-points associated with this elongated equilibrium are assumed to lie in vacuum. The plasma and vacuum regions are bounded by a resistive wall, which is also modeled as an ellipse, with major and minor semiaxes bw and aw, respectively, confocal with the elliptical plasma boundary, i.e., b w 2 a w 2 = b 2 a 2. It is also assumed that equilibrium flows are absent.

The normal mode analysis for n = 0 modes is carried out on the basis of the linearized, reduced ideal-MHD model.28 In the low-beta limit that is standard for tokamak plasmas, velocity flows and magnetic field perturbations along the toroidal direction decouple from the perpendicular dynamics. Thus, the model equations further reduce to a two-field model for the variables ψ and φ. It is convenient to introduce elliptical coordinates ( μ , θ ), which are related to Cartesian coordinates through the transformation relations x = A sinh ( μ ) cos ( θ ) and y = A cosh ( μ ) sin ( θ ), with A = b 2 a 2. The elliptical plasma boundary corresponds to μ = μ b, such that a = A sinh μ b and b = A cosh μ b.

In the plasma region inside the elliptical boundary, the relevant solution for the equilibrium flux function, satisfying 2 ψ e q = J e q ( ψ e q ), after proper normalization is best written in terms of Cartesian components as
ψ e q ( x , y ) = 1 2 ( x 2 b 2 + y 2 a 2 ) .
(1)
The elliptical boundary μ = μ b corresponds to the magnetic surface with constant ψ e q = ψ b = 1 / 2. In the vacuum region outside the elliptical boundary, where μ > μ b, the equilibrium flux ψ e q = ψ e q + satisfies 2 ψ e q + = 0. The superscripts and + indicate the plasma and vacuum regions, respectively. As we assume no equilibrium current sheets, ψeq and its derivative along the normal to the boundary must be continuous across the boundary. The relevant analytic solution is29,30
ψ e q + ( μ , θ ) = 1 2 + α 2 { μ μ b + e 0 2 sinh [ 2 ( μ μ b ) ] cos ( 2 θ ) } ,
(2)
with α 2 = a b / r 0 2,
e 0 = κ 2 1 κ 2 + 1
(3)
the ellipticity of the plasma boundary, and κ = b / a its elongation. Magnetic flux surfaces ψ e q ( μ , θ ) = const exhibit a magnetic separatrix at ψ e q ( μ , θ ) = ψ X = μ b α 2, with X-points located at μ = μ X = 2 μ b and θ = θ X = π / 2 ± n π.
Small perturbations, denoted by an over-tilde, are assumed to vary as
φ ̃ ( μ , θ , t ) = φ ̂ ( μ ) cos θ exp ( γ t ) , ψ ̃ ( μ , θ , t ) = ψ ̂ ( μ ) sin θ exp ( γ t ) ,
with γ = i ω a complex eigenvalue. After straightforward derivation,18 we find the following solution for the stream function within the region bounded by the elliptical surface μ = μ b:
φ ̂ ( μ ) = γ ξ a sinh μ sinh μ b ,
(4)
which corresponds to a rigid vertical shift of the plasma column with constant amplitude ξ. The corresponding perturbed magnetic flux is
ψ ̂ ( μ ) = ξ b cosh μ cosh μ b .
(5)
The eigenvalue γ in the so-called thin wall limit satisfies the cubic dispersion relation19,27
γ 3 + γ 2 ( 1 e ̂ 0 D ) τ η w + γ γ 2 D 1 1 e 0 ̂ D γ 2 ( 1 e 0 ̂ D ) τ η w = 0 ,
(6)
where
γ = 2 κ κ 2 + 1 ( κ 1 κ ) 1 / 2 τ A 1 = ( 1 e 0 ) 1 / 2 [ 1 + e 0 ( 1 e 0 2 ) 1 / 2 ] 1 / 2 τ A 1 ,
(7)
where τ A = ( 4 π ρ m ) 1 / 2 / B p is the relevant Alfvén time, with B P being the radial derivative of the poloidal magnetic field on the magnetic axis, τ η w is the resistive wall time, e ̂ 0 = e 0 κ / ( κ + 1 ), and
D ( κ , b b w ) = κ 2 + 1 ( κ 1 ) 2 { 1 [ 1 κ 2 1 κ 2 ( b b w ) 2 ] 1 / 2 }
(8)
is a geometrical wall parameter, which depends on the elongation κ and on the distance between the plasma and the wall represented by the parameter b / b w. Three relevant limits for D ( κ , b / b w ) are (i) the no-wall limit, where b / b w 0 and D 0; (ii) the circular plasma limit,22  κ 1, where D in Eq. (8) diverges as D ( κ , b / b w ) 2 ( b / b w ) 2 / ( κ 1 ); finally, the limit where the wall approaches the plasma boundary, b / b w 1, and D reaches its maximum value (for a given elongation κ), D ( κ , b / b w ) D max = ( κ 2 + 1 ) / [ κ ( κ 1 ) ], which is always larger than unity.
In the ideal wall limit, τ η w , two terms of the dispersion relation (6) vanish, one solution of the cubic dispersion relation is simply γ = 0, and the other two solutions have real γ 2, which is positive for D < 1, corresponding to the regime where vertical displacements are ideally unstable, and negative when D > 1, corresponding to the VDOM regime. Thus, in this limit,
γ 0 2 = 1 D 1 e ̂ 0 D γ 2 .
(9)
Clearly, γ corresponds to the ideal-MHD growth rate of the vertical displacement in the no-wall limit, D 0. When D > 1, it is convenient to introduce ω = i γ. Then, vertical displacements growing on the ideal-MHD time scale are suppressed by passive wall stabilization, and the VDOM solution is recovered, whose oscillation frequency is
ω = ± ω 0 = ± [ D 1 1 e ̂ 0 D ] 1 / 2 γ .
(10)
Ideal-MHD marginal stability corresponds to D = 1. For the particular wall and plasma geometry assumed by our analytic model, D equals unity when the elliptical wall intercepts the X-points of the equilibrium flux function (corresponding to the up–down symmetric double-null divertor configuration). Indeed, it can be easily checked27 that D ( κ , b / b w ) = 1 when bw = bX, where
b X = κ 2 + 1 κ ( κ 2 1 ) 1 / 2 b
(11)
is the vertical distance of the X-points from the magnetic axis.

Plots of ω 2 τ A 2 as function of κ and b / b w are shown in Sec. IV here below.

Considering finite wall resistivity in the relevant limits D > 1 and ω 0 τ η w 1, one finds that the two oscillatory roots, ω ± ω 0, are weakly damped by wall resistivity, while the zero-frequency root is purely growing on the resistive wall time scale. In tokamak experiments, active feedback stabilization by means of external currents is used to suppress the n = 0 resistive wall mode. However, in this article, we focus our attention on the two other roots of the n = 0 dispersion relation in the limit of an ideal wall.

So far, the analytic results described above assume that the X-points of the equilibrium flux function lie in vacuum. If the hot plasma density extends to the magnetic separatrix, axisymmetric X-point currents can be driven. It was shown in Ref. 17 that these currents can suppress the growth of vertical displacements on the ideal-MHD Alfvén time scale, even if the ideal wall is moved to infinity (D = 0). From the numerical point of view, this requires a very accurate treatment of the plasma response in the vicinity of the magnetic X-points, which will be the subject of a future publication. Nevertheless, in Sec. IV of this article, where the X-points are modeled to lie in a low-density, high-resistivity halo plasma, we will show evidence of the existence of X-point currents carried by the halo plasma.

NIMROD is a parallel three-dimensional initial value code capable of simulating various ideal and non-ideal MHD phenomena in magnetically confined plasmas for both toroidal and linear configurations.21 It uses high order quadrilateral finite element for the poloidal plane and pseudospectral technique in the periodic axisymmetric direction. NIMROD can advance both linear and nonlinear extended MHD equations using implicit/semi-implicit time-advance methods to handle multi-time scales temporal stiffness. The simulations presented in this article advance the linearized version of the single fluid resistive MHD equations
n t + · ( n v ) = · D n ,
(12)
m n ( t + v · ) v = 1 μ 0 ( × B ) × B p + · ν ρ v ,
(13)
3 2 ( t + v · ) T = n T · v ,
(14)
B t = × [ η e ( × B ) v × B ] + κ divb · B ,
(15)
with particle density n and ion mass m, v is the center-of-mass velocity of plasma, p is the total pressure of electrons and ions, T is the single fluid temperature, and the magnetic field is denoted by B. Several explicit diffusive terms are also included: the density diffusivity D, kinetic viscosity ν, electrical diffusivity coefficient ηe (=resistivity over the vacuum permeability μ0), and in the induction equation, κdivb, a diffusivity used to control the · B error.

For simplicity of NIMROD implementation, the analytic straight tokamak is represented by a Cartesian rectangular domain in the poloidal plane and a single toroidal n = 0 mode in the axisymmetric z-direction. This differs from the confocal elliptical wall in the analytic model. Instead, the four sides of the rectangular simulation domain are constrained by imposing the confocal condition, b w 2 a w 2 = b 2 a 2, so that the half-height, bw, and the half width, aw, of the rectangular boundary match the major and minor axes of a confocal elliptical boundary. The rectangular simulation domain circumscribes the analytic confocal ellipse and as such, for equal values of the parameter b / b w, the simulation domain is larger than the analytic one. As we will discuss in Sec. IV, this difference in domain geometry is a likely source of discrepancy between the numeric and analytic results, particularly as the wall is brought closer and the difference in area more pronounced.

The analytic equilibrium for a straight tokamak reviewed in Sec. II is implemented as the NIMROD equilibrium fields. The force balance relation in a linear plasma column takes the following form in the flux coordinate:
d d ψ [ B z 2 2 μ 0 + p ( ψ ) ] = J z ( ψ ) = J 0 .
(16)
The axial magnetic field Bz is uniform and constant throughout the domain. The Heaviside plasma current density J z ( ψ ) = J 0 is a uniform finite value within the hot plasma and zero outside. The pressure profile is a linear function of flux: p ( ψ ) = p 0 ( 1 ψ / ψ b ) = J 0 ψ b ( 1 ψ / ψ b ). Note that ψb (the hot plasma boundary) is not the separatrix. The magnetic geometry is elongated by two external currents Iext located equidistant above and below the plasma: I ext [ δ ( x , y l ) + δ ( x , y + l ) ], located outside the simulation domain. Plasma density is n o = 1.0 × 10 19 m 3 inside the elliptical plasma boundary and ramped down to the halo density nh with a hyperbolic tangent function: n ( ψ ) = n h + ( n o n h ) [ 1 tanh [ σ ( ψ ψ b ) / ψ b ] ] / 2, where σ controls the width of the transition. The equilibrium satisfies the force-balance relation (16) up to the input numerical tolerance and is preserved throughout the linear simulations, i.e., only the perturbed n = 0 component is advanced in time.

NIMROD does not support a vacuum model and instead utilizes a “halo” plasma model to represent the analytic vacuum. The “halo” plasma is a cold region with low density and high resistivity, ηhalo; a few orders of magnitude difference from hot core plasma. This “halo” plasma fills the entire domain outside of the hot plasma (see Fig. 1), both between the hot plasma and separatrix and throughout the open field line region. A measure of the halo resistivity is represented by the dimensionless Lundquist number, S halo = τ R / τ A η halo 1, where the resistive diffusion time τ R = a b / η halo with a being the semi-minor axis and b the semi-major axis, and the Alfvén time τ A = μ 0 ρ m 0 a b / B 0 with B0 and ρ m 0 being the magnetic field and mass density at the magnetic axis.

FIG. 1.

Flux contours in the rectangular simulation domain with the Heaviside current density in red bound by an ellipse with semiaxes (a,b), showing the extent of the hot plasma. Note that the hot plasma ends before the separatrix, with X-points located above and below the hot plasma. A cold, highly resistive halo plasma fills the region between the hot plasma boundary and the perfectly conducting rectangular wall. The three different types of wall used in our work, for b / b w = 0.4, are also shown in the figure a confocal elliptical wall, denoted by the blue dashed ellipse; the green dotted rectangle, with b w / b = a w / a; and the outer black solid rectangle with b w 2 a w 2 = b 2 a 2 as for the confocal elliptical wall. For this particular figure, κ = 2.0.

FIG. 1.

Flux contours in the rectangular simulation domain with the Heaviside current density in red bound by an ellipse with semiaxes (a,b), showing the extent of the hot plasma. Note that the hot plasma ends before the separatrix, with X-points located above and below the hot plasma. A cold, highly resistive halo plasma fills the region between the hot plasma boundary and the perfectly conducting rectangular wall. The three different types of wall used in our work, for b / b w = 0.4, are also shown in the figure a confocal elliptical wall, denoted by the blue dashed ellipse; the green dotted rectangle, with b w / b = a w / a; and the outer black solid rectangle with b w 2 a w 2 = b 2 a 2 as for the confocal elliptical wall. For this particular figure, κ = 2.0.

Close modal

In this section, we present the benchmark comparison between the analytic theory results summarized in Sec. II and NIMROD numerical simulations. In numerical work, the wall is assumed to be ideal. In NIMROD, ideal conducting walls are implemented by enforcing vanishing normal components of the perturbed magnetic field at the wall.

First, we will compare unstable vertical displacements and scans of wall position and ellipticity. Then, we will consider vertical displacement oscillatory modes (VDOMs). In analytic work, where vacuum fills the space between the plasma boundary and the wall, and in the absence of dissipation and fast ions, VDOMs are undamped, purely oscillatory modes. In NIMROD, where finite dissipation is usually required for numerical stability, and a low-density, high-resistivity halo plasma is used instead of vacuum, some degree of damping is always found. This damping is shown to saturate to small values when the halo plasma approach vacuum parameters.

In most simulations presented, the grid resolution is (mx,my) = (360,360) with polynomial degree = 3 and the time step Δ t = 0.5 τ A. The density diffusion τ D / τ A O ( 10 2 ) ( τ D = a b / D); viscosity τ ν / τ A O ( 10 4 ) ( τ ν = a b / ν); resistivity S O ( 10 6 ) inside the elliptical plasma; and · B diffusion τ κ divb / τ A O ( 10 1 ) ( τ κ divb = a b / κ divb) are found sufficient for numerical convergence. However, in some cases, as indicated below, a more refined mesh was necessary to obtain convergence.

This first scan examines the growth rate of the plasma unbound by any conducting wall, i.e., the no-wall limit, as the plasma ellipticity is varied. NIMROD with conducting wall boundaries does not simulate the analytic no-wall limit directly. Instead, for each value of ellipticity, NIMROD estimates the no-wall growth rate by performing successive simulations moving the domain boundary further away until a saturated value for the growth rate is obtained. As shown in Fig. 2(b), the no-wall growth rate for e 0 = 0.3 is reached in simulation for values of b w 6 b. For e 0 = [ 0.2 0.6 ], the wall parameter b / b w necessary for convergence are in the range (0.13, 0.19).

FIG. 2.

(a) Comparison of the normalized growth rate vs plasma ellipticity shows agreement between NIMROD simulations and the analytic no-wall limit. (b) The growth rate gradually saturates to the analytical no-wall value with increasing wall-plasma distance ( e 0 = 0.3).

FIG. 2.

(a) Comparison of the normalized growth rate vs plasma ellipticity shows agreement between NIMROD simulations and the analytic no-wall limit. (b) The growth rate gradually saturates to the analytical no-wall value with increasing wall-plasma distance ( e 0 = 0.3).

Close modal

Figure 2(a) shows good agreement between NIMROD linear simulations and analytic results. However, agreement is better at larger values of ellipticity. For e 0 = 0.1, the growth rate is not quite converged. The highest value b / b w = 0.13 was used. Convergence is constrained by the need to include and resolve the X-points in the ever increasing domain. As ellipticity is decreased and approaches 0, the X-points retreat away from the plasma and approach infinity, requiring ever larger simulation domains, the majority of which is then composed of halo plasma and open field lines. Values of e 0 < 0.1 become computationally prohibitive; the e0 = 0 circular limit is numerically out of reach.

The growth rate in every run is calculated only after the perturbed fields have reached a pure exponentially growing phase after an initial transient. A typical example for the perturbed x-component of the magnetic field, B ̃ x, from an unstable run is plotted in Fig. 3.

FIG. 3.

Time history of the perturbed field B ̃ x for the unstable no-wall case with e 0 = 0.3. After an initial transients, pure exponential growth is observed.

FIG. 3.

Time history of the perturbed field B ̃ x for the unstable no-wall case with e 0 = 0.3. After an initial transients, pure exponential growth is observed.

Close modal

Figure 4(a) shows perturbed pressure contours and Fig. 4(b) 2D arrow plots indicating the direction of the perturbed velocity projected onto a poloidal cross section. The pressure contours show a clear m = 1 signature caused by the upward movement of the plasma column. The arrow plot indicates the direction of flow in the poloidal plane. Within the hot plasma, all arrows point upward, consistent with a nearly rigid vertical shift of the entire plasma and in very good agreement with the mode structure obtained analytically. Outside the hot plasma, in the halo region, the arrows show the direction of the return flow. In analytic theory, where vacuum is assumed instead of halo plasma, the return flow is concentrated on the plasma boundary, giving rise to localized vorticity sheets.

FIG. 4.

Perturbed pressure and plasma flow arrows plotted on equilibrium flux contours show the vertical rigid shift of the plasma. (a) Perturbed pressure. (b) Direction of perturbed flow.

FIG. 4.

Perturbed pressure and plasma flow arrows plotted on equilibrium flux contours show the vertical rigid shift of the plasma. (a) Perturbed pressure. (b) Direction of perturbed flow.

Close modal

We next examine the impact of moving the wall closer to the plasma. The results of a confocal wall scan is presented for two plasma elongations, κ = 2.0 and 1.4, corresponding to ellipticity values e 0 = 0.6 and e 0 = 0.32, respectively. For these two scans, we preserve the confocal constraints ( b w 2 a w 2 = b 2 a 2) applied to the rectangular domain boundary. Figure 5, the confocal wall scan, plots the square of the normalized frequency vs wall parameter, b / b w, and shows agreement between NIMROD and analytic theory. Positive values (shown in circles) indicate oscillating modes. Negative values (shown in triangles) indicate growing modes. The no-wall limit corresponds to values of b / b w 0, while the wall touches the plasma when b / b w = 1. The zero crossing of the analytic theory curve (green line) indicates transition between purely oscillating and purely growing mode, which occurs at bw = bX, with bX the vertical distance of the X-points from the magnetic axis. In other words, marginal stability ( ω 2 = 0) occurs when the confocal elliptical wall intercepts the X-points. Starting from small values of b / b w, the growth rate, γ = i ω, decreases as the wall is brought closer to the plasma (increasing b / b w), until the wall approaches the X-point. With the proximity of the wall to the X-point, the mode becomes marginal, and then transitions to a purely oscillating mode for values of b / b w > b / b X, as the wall is brought within the X-points, i.e., as the X-points lie outside the simulation domain.

FIG. 5.

Confocal wall scan ( b w 2 a w 2 = b 2 a 2) plots the square of the normalized frequency vs wall parameter b / b w for κ = 2.0 (left panel) and κ = 1.4 (right panel), showing agreement between NIMROD and analytic theory. Positive values (circles) indicate oscillating modes. Negative values (triangles) indicate growing modes. The zero crossing of the analytic theory curve (green line) occurs for b / b w = b / b X, where the domain boundary intersects the X-points.

FIG. 5.

Confocal wall scan ( b w 2 a w 2 = b 2 a 2) plots the square of the normalized frequency vs wall parameter b / b w for κ = 2.0 (left panel) and κ = 1.4 (right panel), showing agreement between NIMROD and analytic theory. Positive values (circles) indicate oscillating modes. Negative values (triangles) indicate growing modes. The zero crossing of the analytic theory curve (green line) occurs for b / b w = b / b X, where the domain boundary intersects the X-points.

Close modal

In the simulations, the transition from the unstable growing mode to stable oscillatory mode does not occur exactly when the rectangular wall intercepts the X-points, as in analytic theory. We also note that numerical points lie slightly below the analytic curve, and that this discrepancy is more pronounced at larger values of b / b w, where the wall is closer to the plasma boundary. We believe that this is due to the use of the rectangular simulation domain, which has a larger area than the analytic confocal elliptical boundary for equal values of the parameter b / b w. Therefore, the rectangular wall is on average a bit further away from the plasma boundary than the equivalent confocal elliptical wall. Hence, the stabilizing effect exerted by the rectangular wall is a little weaker and the growth rate found numerically is indeed a bit larger (i.e., smaller negative ω 2) than the analytic result.

To test this idea, in Fig. 6 we present another wall scan for κ = 2.0, where we adopted a self-similar wall, as shown in Fig. 1 as green dashed-dotted rectangle, i.e., b w / b = a w / a, instead of the confocal constraint. Again, we see agreement between simulation and analytic theory. As expected, the numerical results for the self-similar wall lie slightly above the analytic curve, whereas in the prior two wall scans with confocal walls, the simulation results lie below the analytic curve. Again, the discrepancy is larger at larger values of b / b w. We believe that this is due to the fact that, in this case, the self-similar wall is on average closer to the plasma than the confocal wall for equal values of b / b w, and so the stabilizing influence of the wall is larger in this case.

FIG. 6.

A self-similar wall scan ( b / b w = a / a w) plots the square of the normalized frequency vs wall parameter b / b w, showing agreement between analytic theory and NIMROD simulations.

FIG. 6.

A self-similar wall scan ( b / b w = a / a w) plots the square of the normalized frequency vs wall parameter b / b w, showing agreement between analytic theory and NIMROD simulations.

Close modal

A scan of the plasma elongation is performed for values of κ between κ = 1 and κ = 2 (e0 between 0 and 0.6), with wall parameter fixed at b / b w = 0.25. Figure 7 shows the squared normalized frequency plotted against κ; oscillatory cases are denoted by circles and unstable cases by triangles. The zero crossing of the analytic curve occurs for the value of κ = κ X corresponding to the elliptical wall intersecting the X-points. This value can be obtained from Eq. (11), where specifically we set b / b X = b / b w = 0.25, yielding κ X 1.12. Agreement between analytic theory and numerical results is excellent, as the wall at b w = 4 b is relatively far from the plasma, and so the effect due to the difference between elliptical and rectangular wall is negligible in this case.

FIG. 7.

Plasma elongation scan, with wall parameter fixed at b / b w = 0.25, plots the square of the normalized frequency vs elongation, showing agreement between analytic theory and NIMROD.

FIG. 7.

Plasma elongation scan, with wall parameter fixed at b / b w = 0.25, plots the square of the normalized frequency vs elongation, showing agreement between analytic theory and NIMROD.

Close modal

As the conducting wall is brought closer to the plasma, so that the X-points lie beyond the wall and outside the integration domain, the mode transitions from unstable to a stable oscillatory mode. These stable oscillations were dubbed VDOMs, which stands for vertical displacement oscillatory modes.19 

Here, we examine one such stable oscillatory case. Figure 8(a) shows the time history of perturbed fields: velocity ( v ̃ y), magnetic field ( B ̃ x), density ( n ̃), temperature ( T ̃), and Fig. 8(b) midplane profiles of vertical momentum per unit mass, n 0 v ̃ y, at several times during a single period, for κ = 1.4 and b / b w = 0.55. The single oscillation period is denoted by the yellow bar in the time history. The normalized oscillation frequency is ω τ A = 0.31.

FIG. 8.

Time history of perturbed fields: velocity ( v ̃ y), magnetic field ( B ̃ x), density ( n ̃), temperature ( T ̃); midplane profiles of vertical momentum per unit mass at several times during a single oscillation period (denoted by the yellow bar in the time history); κ = 1.4 and b / b w = 0.55. The plasma motion is well approximated by a vertical rigid shift.

FIG. 8.

Time history of perturbed fields: velocity ( v ̃ y), magnetic field ( B ̃ x), density ( n ̃), temperature ( T ̃); midplane profiles of vertical momentum per unit mass at several times during a single oscillation period (denoted by the yellow bar in the time history); κ = 1.4 and b / b w = 0.55. The plasma motion is well approximated by a vertical rigid shift.

Close modal

The time slices of the midplane profile show that the n = 0 oscillation is mostly confined to within the hot plasma, with a perturbed mass flow that well approximates the rigid-shift vertical oscillation found in analytic work. Note that the plasma extends to x / b = a / b 0.7.

Recall that NIMROD does not support a true vacuum, but instead employs a low-density, high-resistivity halo plasma in the region surrounding the hot plasma. As a result, consistent with the flow patterns shown in Fig. 4(b), a return flow is carried by the halo plasma, which guarantees near-incompressibility, · v = 0. However, owing to the low halo density, the return mass flow carried by the halo plasma is very small compared with the mass flow in the hot plasma.

Simulations show that the oscillation is damped due to the finite halo density and resistivity. Figure 9(a) is a scan of halo plasma density nh and Fig. 9(b) Lundquist number Shalo vs damping. For the halo density scan, the Lundquist number is fixed at Shalo = 0.4. A converged damping rate is approached for n h < O ( 10 17 ) m 3. For the Lundquist number scan, we fix the density at n h = 1 × 10 17 m 3. A converged damping rate is approached for S halo < O ( 10 1 ).

FIG. 9.

A scan of halo plasma density nh and Lundquist number Shalo vs damping rate, λ = I m ( ω ).

FIG. 9.

A scan of halo plasma density nh and Lundquist number Shalo vs damping rate, λ = I m ( ω ).

Close modal

These simulations required a more refined mesh with (mx,my) = (480,480) and polynomial degree = 3.

Here, we discuss an ideal-MHD unstable n = 0 vertical displacement mode with κ = 2.0 and b / b w = 0.4 and focus on the currents that form in response to the growing vertical motion of the plasma. Figure 10(a) shows perturbed surface currents at the plasma boundary, as predicted by analytic work. The surface current density is very narrow, a consequence of having chosen a step-function equilibrium current density with a sharp drop at the plasma boundary.

FIG. 10.

(a) Contour plots of the perturbed current density, J ̃ z, show current formation at the X-points and surface currents along the hot plasma boundary. (b) Profiles of normalized J ̃ z / v ̃ y 0 along a vertical line segment through the lower X-point [shown by black dashed line in part (a)]; the X-point location is denoted by the dotted line, for different values of Shalo; v ̃ y 0 is the value of perturbed v ̃ y at the origin. The peak of the X-point current density increases with increasing Shalo. In both panels, κ = 2.0 and b / b w 0.4.

FIG. 10.

(a) Contour plots of the perturbed current density, J ̃ z, show current formation at the X-points and surface currents along the hot plasma boundary. (b) Profiles of normalized J ̃ z / v ̃ y 0 along a vertical line segment through the lower X-point [shown by black dashed line in part (a)]; the X-point location is denoted by the dotted line, for different values of Shalo; v ̃ y 0 is the value of perturbed v ̃ y at the origin. The peak of the X-point current density increases with increasing Shalo. In both panels, κ = 2.0 and b / b w 0.4.

Close modal

Of greater interest to us are the X-point currents, as we believe that this is an interesting result with potential consequences on X-point topology and plasma edge MHD activity. Figure 10(a) shows a 2D color plot and Fig. 10(b) 1D profiles along the lower vertical midplane at x = 0 of the perturbed axial current density, J ̃ z normalized by the value of perturbed v ̃ y at the origin. The color plot shows both surface and X-point currents, having opposite polarity, as expected. The 1D profiles of J ̃ z are drawn for different values of the halo Lundquist number. These profiles become sharper and narrower as halo resistivity is reduced, which is in qualitative agreement with analytic predictions.17 More specifically, as shown in Ref. 17, the perturbed X-point current density becomes a delta function peaked at the resonant X-point in the limit of vanishing resistivity.

There is, however, an important difference between our NIMROD simulations and the analytic work. In these NIMROD simulations, the X-points lie in a halo plasma and the X-point current is sensitive to halo resistivity and halo density. By contrast, in analytic work,17 the X-points lie in a vacuum.

Figure 11 shows the behavior of the growth rate as a function of halo resistivity. In a relatively narrow range, we find that the growth rate decreases as the Lundquist number increases. This suggests a stabilizing effect from the X-point currents on vertical displacements, but this effect is found to be modest in these simulations. We may expect a more pronounced effect in situations where the hot plasma extends to the X-points and the magnetic separatrix.17 The growth rate saturates both at the lower end and at the upper end of the resistivity scan. The saturated growth rate at low values of the Lundquist number (high resistivity), where the halo resistivity approximates vacuum, is in good agreement with analytic prediction.18,19 The saturated growth rate at the higher end has no analytic comparison, as halo effects, which were not considered in analytic work, are important in determining that saturated values.

FIG. 11.

Halo resistivity vs normalized growth rate. κ = 2.0 and b / b w = 0.4.

FIG. 11.

Halo resistivity vs normalized growth rate. κ = 2.0 and b / b w = 0.4.

Close modal

The main result of this article is that we have successfully reproduced recent analytic results on n = 0 modes using the extended-MHD code NIMROD. The agreement between analytic theory and numerical results is very satisfactory, providing a successful benchmark and a useful starting point for future numerical investigations of n = 0 modes using more realistic tokamak geometry and plasma equilibria. Differences between analytic and numerical results have been interpreted as due mostly to the different shape of the wall used in analytic work vs NIMROD simulation, and the fact that NIMROD adopts a low-density, high-resistivity halo plasma in lieu of vacuum assumed in analytic work.

When ideal-MHD vertical displacements are passively stabilized due to the proximity of the ideal wall, vertical displacement oscillatory modes, dubbed VDOMs in Ref. 19, are found. In the present simulations, VDOMs are weakly damped, due to low values of plasma and halo resistivity. Viscosity and other dissipative terms also affect the damping rate, but at the low values of these parameters used in the simulations (see Sec. IV, before Sec. IV A), their effects on damping is negligibly small. The VDOM frequency is Alfvénic, and for realistic values of the wall parameter D defined in Eq. (8), it falls slightly below the poloidal Alfvén frequency, which should make these modes immune from continuum damping when an actual toroidal geometry with realistic equilibrium current density profiles is considered. VDOMs are of great interest as they can be driven unstable by their resonant interaction with fast ions. Indeed, we have proposed20 that recent observations15,16 of n = 0 modes on JET tokamak experiments with a significant energetic ion population can be interpreted as fast-ion-driven vertical modes. Work is in progress studying recent JET discharges to test this hypothesis. Ongoing linear toroidal simulations, using the NIMROD code, aim to extend our analysis to more realistic plasma geometries and profiles, starting from reconstructed equilibria of the relevant JET discharges. We point out that these JET discharges are single-null configurations, and yet, especially as far as VDOMs are concerned, preliminary NIMROD simulations are in good qualitative agreement with what is reported in the present article. These results will be presented in a future publication.

This article also discusses for the first time numerical results related to the occurrence of axisymmetric X-point currents. In this article, X-point currents are supported by the halo plasma. In future work, where we plan to extend the hot plasma to the magnetic separatrix, X-point currents may become more pronounced, leading to a more significant impact on the topology of the X-point region, as the analytic work of Refs. 17 and 18 suggests. We point out that treating X-point effects on n = 0 perturbation is numerically challenging, since the toroidal field line going through a magnetic X-point is resonant (in the sense that B eq · k = 0 at magnetic X-points), which is the reason why axisymmetric current sheets localized near X-points are to be expected.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom Research and Training Programme (Grant Agreement No. 101052200–EUROfusion). A part of this work was carried out using the EUROfusion High Performance Computer (Marconi-Fusion) under the Project No. FUA 37 _ NSAM. The views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

This study was carried out within the Ministerial Decree No. 1062/2021 and received funding from the FSE REACT-EU - PON Ricerca e Innovazione 2014–2020. This manuscript reflects only the authors' views and opinions, neither the European Union nor the European Commission can be considered responsible for them.

We would also like to gratefully acknowledge Eni S.p.A. for the use of their High Performance Computing facility (HPC4), where NIMROD simulations were carried out, and Claudio Carati of Eni S.p.A. for his interest in our work. The author D.B. expresses his gratitude to the partial support from Eni S.p.A for his current researcher position.

The authors have no conflicts to disclose.

Debabrata Banerjee: Conceptualization (equal); Data curation (equal); Formal analysis (supporting); Funding acquisition (supporting); Investigation (lead); Methodology (equal); Software (equal); Validation (equal); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Charlson C. Kim: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Tommaso Barberis: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (equal); Methodology (equal); Software (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Francesco Porcelli: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
G.
Laval
,
R.
Pellat
, and
J. L.
Soule
, “
Hydromagnetic stability of a current-carrying pinch with noncircular cross section
,”
Phys. Fluids
17
,
835
(
1974
).
2.
M.
Okabayashi
and
G.
Sheffield
, “
Vertical stability of elongated tokamaks
,”
Nucl. Fusion
14
,
263
265
(
1974
).
3.
E.
Rebhan
, “
Stability boundaries of tokamaks with respect to rigid displacements
,”
Nucl. Fusion
15
,
277
(
1975
).
4.
F. A.
Haas
, “
Stability of a high-β tokamak to uniform vertical displacements
,”
Nucl. Fusion
15
,
407
(
1975
).
5.
M.
Perrone
and
J.
Wesson
, “
Stability of axisymmetric modes in JET
,”
Nucl. Fusion
21
,
871
(
1981
).
6.
R.
Fitzpatrick
, “
A simple ideal magnetohydrodynamical model of vertical disruption events in tokamaks
,”
Phys. Plasmas
16
,
012506
(
2009
).
7.
L. E.
Zakharov
,
S. A.
Galkin
,
S. N.
Gerasimov
, and
JET-EFDA Contributors,
Understanding disruptions in tokamaks
,”
Phys. Plasmas
19
,
055703
(
2012
).
8.
C. F.
Clauser
,
S. C.
Jardin
, and
N. M.
Ferraro
, “
Vertical forces during vertical displacement events in an ITER plasma and the role of halo currents
,”
Nucl. Fusion
59
,
126037
(
2019
).
9.
E. A.
Lazarus
,
J. B.
Lister
, and
G. H.
Neilson
, “
Control of the vertical instability in tokamaks
,”
Nucl. Fusion
30
,
111
(
1990
).
10.
J. B.
Lister
,
E. A.
Lazarus
,
A. G.
Kellman
,
J.-M.
Moret
,
J. R.
Ferron
,
F. J.
Helton
,
L. L.
Lao
,
J. A.
Leuer
,
E. J.
Strait
,
T. S.
Taylor
, and
A. D.
Turnbull
, “
Experimental study of the vertical stability of high decay index plasmas in the DIII-D tokamak
,”
Nucl. Fusion
30
,
2349
(
1990
).
11.
I.
Hutchinson
, “
Simplified models of axisymmetric magnetohydrodynamic instabilities
,”
Nucl. Fusion
29
,
2107
(
1989
).
12.
M.
Albanese
,
R.
Mattei
, and
F.
Villone
, “
Prediction of the growth rates of VDEs in JET
,”
Nucl. Fusion
44
,
999
(
2004
).
13.
F.
Villone
,
V.
Riccardo
,
F.
Sartori
, and
JET-EFDA Contributors,
Position dependence of the unstable vertical movement of JET plasmas triggered by ELMs
,”
Nucl. Fusion
45
,
1328
(
2005
).
14.
A.
Portone
, “
Active and passive stabilization of n = 0 RWMs in future tokamak devices
,”
Nucl. Fusion
57
,
126060
(
2017
).
15.
H. J. C.
Oliver
,
S. E.
Sharapov
,
B. N.
Breizman
, and
L.-J.
Zheng
, “
Axisymmetric global Alfvén eigenmodes within the ellipticity-induced frequency gap in the Joint European Torus
,”
Phys. Plasmas
24
,
122505
(
2017
).
16.
V.
Kiptily
,
M.
Fitzgerald
,
Y.
Kazakov
,
J.
Ongena
,
M.
Nocente
,
S.
Sharapov
,
M.
Dreval
,
Š.
Štancar
,
T.
Craciunescu
,
J.
Garcia
,
L.
Giacomelli
,
V.
Goloborodko
,
H.
Oliver
,
H.
Weisen
, and
J. Contributors,
Evidence for Alfvén eigenmodes driven by alpha particles in D-He3 fusion experiments on JET
,”
Nucl. Fusion
61
,
114006
(
2021
).
17.
A.
Yolbarsop
,
F.
Porcelli
, and
R.
Fitzpatrick
, “
Impact of magnetic X-points on the vertical stability of tokamak plasmas
,”
Nucl. Fusion
61
,
114003
(
2021
).
18.
A.
Yolbarsop
,
F.
Porcelli
,
W.
Liu
, and
R.
Fitzpatrick
, “
Analytic theory of ideal-MHD vertical displacements in tokamak plasmas
,”
Plasma Phys. Controlled Fusion
64
,
105002
(
2022
).
19.
T.
Barberis
,
A.
Yolbarsop
, and
F.
Porcelli
, “
Vertical displacement oscillatory modes in tokamak plasmas
,”
J. Plasma Phys.
88
,
905880511
(
2022
).
20.
T.
Barberis
,
F.
Porcelli
, and
A.
Yolbarsop
, “
Fast-ion-driven vertical mods in magnetically confined toroidal plasmas
,”
Nucl. Fusion
62
,
064002
(
2022
).
21.
C.
Sovinec
,
A.
Glasser
,
T.
Gianakon
,
D.
Barnes
,
R.
Nebel
,
S.
Kruger
,
S.
Plimpton
,
A.
Tarditi
,
M.
Chu
, and
NIMROD Team,
J. Comput. Phys.
195
,
355
(
2004
).
22.
A.
Yolbarsop
,
F.
Porcelli
,
D.
Banerjee
,
C. C.
Kim
, and
L.
Hong
, “
Axisymmetric oscillatory modes in cylindrical magnetized plasma bounded by a conducting wall
,”
Phys. Lett. A
479
,
128940
(
2023
).
23.
L.
Villard
and
J.
Vaclavik
, “
Alfven frequency modes and global Alfven eigenmodes
,”
Nucl. Fusion
37
,
351
(
1997
).
24.
N.
Winsor
,
J. L.
Johnson
, and
J. M.
Dawson
, “
Geodesic acoustic waves in hydromagnetic systems
,”
Phys. Fluids
11
,
2448
2450
(
1968
).
25.
P. H.
Diamond
,
S.-I.
Itoh
,
K.
Itoh
, and
T. S.
Hahm
, “
Zonal flows in plasma-A review
,”
Plasma Phys. Controlled Fusion
47
,
R35
(
2005
).
26.
G. Y.
Fu
, “
Energetic-particle-induced geodesic acoustic mode
,”
Phys. Rev. Lett.
101
,
185002
(
2008
).
27.
F.
Porcelli
,
T.
Barberis
, and
A.
Yolbarsop
, “
Vertical displacements close to ideal-MHD marginal stability in tokamak plasmas
,”
Fundam. Plasma Phys.
5
,
100017
(
2023
).
28.
H. R.
Strauss
,
Phys. Fluids
19
,
134
(
1976
).
29.
R.
Gajewski
,
Phys. Fluids
15
,
70
74
(
1972
).
30.
F.
Porcelli
and
A.
Yolbarsop
, “
Analytic equilibrium of “straight tokamak” plasma bounded by a magnetic separatrix
,”
Phys. Plasmas
26
,
054501
(
2019
).