1 Introduction
Magnetically confined plasmas aimed at thermonuclear fusion best perform close to pressure- and current-driven magnetohydrodynamic (MHD) stability boundaries and are thus subject to a wide spectrum of potentially dangerous large-scale relaxation phenomena, e.g. the tearing modes, the vertical instabilities, the kink modes, the edge localised modes (see, e.g. White Reference White2006; Wesson & Campbell Reference Wesson and Campbell2011). The effect of such instabilities is generally deleterious (White Reference White2006), unless stable nearby equilibria may be attained by nonlinear saturation or through externally imposed magnetic perturbations. This is the case in hybrid (Piovesan et al. Reference Piovesan, Bonfiglio, Cianciosa, Luce, Taylor, Terranova, Turco, Wilcox, Wingen, Cappello, Chrystal, Escande, Holcomb, Marrelli, Paz-Soldan, Piron, Predebon and Zaniol2017; Burckhart et al. Reference Burckhart, Bock, Fischer, Pütterich, Stober, Günter, Gude, Hobirk, Hölzl, Igochine, Krebs, Maraschek, Reisner, Schramm and Zohm2023) and advanced (reversed shear) tokamak scenarios (Gormezano et al. Reference Gormezano, Challis, Joffrin, Litaudon and Sips2008; Zhang et al. Reference Zhang, Wu, Li, Ding, Liu, Qian, Gong, Gao, Gao, Wu and Li2021), or in quasi-single-helicity states in the reversed field pinch (RFP) dynamo/flux-pumping effect (Bonfiglio, Cappello & Escande Reference Bonfiglio, Cappello and Escande2005; Marrelli et al. Reference Marrelli, Martin, Puiatti, Sarff, Champan, Drake, Escande and Masamune2021).
Among the most dangerous MHD instabilities, an important class is represented by the free-boundary modes. Characterised by a finite fluid velocity at the plasma edge, they may lead to disruptions in the tokamak configuration (Granetz et al. Reference Granetz, Hutchinson, Sorci, Irby, LaBombard and Gwinn1996; Wesson & Campbell Reference Wesson and Campbell2011); this is the case of vertical instabilities and external kink modes.
Advanced numerical modelling is crucial for understanding such complex and typically nonlinear dynamics. Starting from the 1980s, an increasing attention has been devoted to the formulation of realistic magnetic boundary conditions (BCs), mostly relying on resistive-wall modules (see, e.g. Gimblett Reference Gimblett1986; Hender, Gimblett & Robinson Reference Hender, Gimblett and Robinson1989; Haines, Gimblett & Hastie Reference Haines, Gimblett and Hastie2013) typically implemented at the numerical boundary of MHD codes (as in Schnack et al. Reference Schnack, Barnes, Mikic, Harned and Caramana1987; Strauss et al. Reference Strauss, Pletzer, Park, Jardin, Breslau and Sugiyama2004; Marx & Lütjens Reference Marx and Lütjens2017; Bonotto et al. Reference Bonotto, Liu, Villone, Pigatto and Bettini2020; Bunkers & Sovinec Reference Bunkers and Sovinec2020; Isernia et al. Reference Isernia, Schwarz, Artola, Ventre, Hoelzl, Rubinacci and Villone2023), but not necessarily (Ferraro et al. Reference Ferraro, Jardin, Lao, Shephard and Zhang2016).
The problem of a free (fully consistent) normal velocity boundary has deserved renewed interest, also prompted by some theoretical studies (see Zakharov Reference Zakharov2008; Zakharov, Galkin & Gerasimov Reference Zakharov, Galkin and Gerasimov2012; Zakharov & Li Reference Zakharov and Li2015), suggesting that its absence in most present-day codes may prevent them from reliably modelling the vertical displacement events in tokamak plasmas. Two numerical studies have specifically investigated the impact of the velocity boundary formulation, performed, respectively, with the M3D code (Strauss Reference Strauss2014) and with the NIMROD code (Bunkers & Sovinec Reference Bunkers and Sovinec2020), both demonstrating its critical implications for numerical predictions. Still more recently, realistic velocity BCs have been indicated among the unresolved demanding challenges for the modelling of tokamak vertical displacement events (see Artola et al. Reference Artola, Schwarz, Gerasimov, Loarte and Hoelzl2024).
The nonlinear code-to-code verification study of a resistive wall module with a fully consistent velocity boundary has been recently performed (Spinicci et al. Reference Spinicci, Bonfiglio, Chacón, Cappello and Veranda2023) between the SPECYL (Cappello & Biskamp Reference Cappello and Biskamp1996) and the PIXIE3D codes (Chacón Reference Chacón2008), extending a previous nonlinear verification study between the same two codes, where ideal wall BCs were adopted (see Bonfiglio, Chacón & Cappello Reference Bonfiglio, Chacón and Cappello2010). The plasma interface was modelled as an axisymmetric thin shell (as in Hender et al. Reference Hender, Gimblett and Robinson1989), separated from an outer ideal wall by a vacuum region. The purely perpendicular velocity boundary featured a free normal-flow (NF) and allowed a fully self-consistent coupling within the MHD model. The quantitative reliability of both SPECYL's and PIXIE3D's NF boundary implementations was verified by a remarkably general mutual agreement on several test cases, both in the tokamak and the RFP configurations.
This paper employs the SPECYL code in a comparative study of the fully consistent BCs with a free NF, just outlined in the previous paragraph, against an approximated fixed NF formulation that does not allow finite velocity perturbations at the numerical boundary, often employed in nonlinear codes (analogous BCs are adopted, e.g. in Bonfiglio et al. (Reference Bonfiglio, Chacón and Cappello2010), McAdams (Reference McAdams2014), Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016), Marx & Lütjens (Reference Marx and Lütjens2017) and Ramasamy et al. (Reference Ramasamy, Ramirez, Hoelzl, Graves, López, Lackner and Günter2022)). The reliability of both is measured against well-known figures of merit from the ideal MHD linear theory of external kink modes (see, e.g. Freidberg Reference Freidberg2014). The thin shell of SPECYL's BCs can be set to a high enough resistivity for the magnetic field diffusion through it to be almost instantaneous on the characteristic time scales of the plasma dynamics. In this condition, here dubbed the transparent wall limit, free-boundary modes driven by plasma current are expected to be unstable in the presence of a fully consistent velocity field. We show in this study that only SPECYL's free NF formulation allows proper modelling of the plasma interface through a transparent shell, in optimal agreement with the linear theory.
Remarkably, external kink modes simulations in the transparent-wall limit substantially differ from the typical approach adopted by other codes in the past, that leverages a high-resistivity and low-density pseudovacuum region (PV) to displace the numerical boundary (where BCs are enforced) away from the plasma interface. In the PV, high resistivity and vanishing density strongly damp the plasma flow, so that simplified (even not fully consistent) velocity BCs may be enforced with very low impact on simulations outcomes.
The paper is organised as follows. In § 2 we start with an overview of SPECYL's boundary module. Section 3 presents the numerical set-up, leveraged to pursue both the free NF velocity boundary approach and the simplified velocity boundary with PV. Section 4 contains the main body of our results: we study the external kink mode stability of two axisymmetric equilibria, characterised by a uniform current density and a shaped current density, respectively. Section 5 gives conclusions and outlooks. Four appendices are also present, dedicated to, respectively, an overview of linear kink stability analysis (Appendix A), the analytical profiles and growth rates for the simplified initial equilibrium characterised by a flat current (Appendix B), an overview of the semianalytical tool employed to produce the theoretical figures of merit for the Wesson-like current density profile (Appendix C) and a discussion in relation to the remarkable benchmark illustrated in Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016) (Appendix D).
2 The SPECYL code and its BCs
This section aims at a general overview of the SPECYL code and its boundary module, which is described in detail in Spinicci et al. (Reference Spinicci, Bonfiglio, Chacón, Cappello and Veranda2023).
2.1 The SPECYL code
The SPECYL code computes the numerical solution of nonlinear three-dimensional (3-D) MHD equations by evolving in time $t$
the magnetic field $\boldsymbol {B}$
and the flow velocity field $\boldsymbol {v}$
according to the following viscoresistive scheme:
In the equations above, $\rho$
is the axisymmetric mass density, $\boldsymbol {J}$
is the current density, $\nu$
and $\eta$
are a uniform viscosity and an axisymmetric resistivity, respectively, and $\boldsymbol {E}$
is the electric field.
All quantities appear in dimensionless units: $\rho =\hat {\rho }/\rho _0$
(where $\rho _0=\hat {\rho }|_{r=0}$
); $\boldsymbol {B}=\boldsymbol {\hat {B}}/B_0$
(where $B_0=|\boldsymbol {\hat {B}}|_{r=0{,t=0}}$
); $\boldsymbol {v}=\boldsymbol {\hat {v}}/v_A$
($v_A=B_0/\sqrt {\mu _0\rho _0}$
being the Alfvénic velocity, with $\mu _0$
vacuum permeability); $t=\hat {t}/\tau _A$
($\tau _A=r_{\text {BC}}/v_A$
being the Alfvénic time and $r_{\text {BC}}$
the radial span of the computational domain, coincident with the plasma radius $a$
unless a PV region is included); $\boldsymbol {J}=\hat {\boldsymbol {J}}\mu _0 r_{\text {BC}}/B_0$
; $\boldsymbol {E}=\hat {\boldsymbol {E}}/v_A B_0$
. Hatted quantities are in physical units. In dimensionless units, the scalar kinematic viscosity corresponds to the inverse viscous Lundquist (or Reynold's) number, $\nu =\tau _A/\tau _\nu \equiv M^{-1}$
(where $\tau _\nu =r_{\text {BC}}^2/\nu$
is the viscous time scale); and resistivity corresponds to the inverse Lundquist number, $\eta =\tau _A/\tau _R\equiv S^{-1}$
(where $\tau _R=\mu _0 r_{\text {BC}}^2/\eta$
is the resistive-diffusion time scale).
The geometry consists of a periodic cylinder: a spectral approach is adopted for the periodic coordinates, $0\leq \theta \leq 2{\rm \pi}$
and $0\leq z\leq 2{\rm \pi} R$
($R$
being the plasma major radius), while adopting a finite-differences staggered scheme in the radial direction, $0\leq r\leq r_{\text {BC}}$
. The pressure gradients are neglected, and the axisymmetric density $\rho (r)$
and resistivity $\eta (r),$
along with the uniform scalar viscosity $\nu$
are kept constant in time. No pressure or temperature dependence is assumed for $\rho$
and $\eta$
, thus presently limiting SPECYL's ability to accurately model highly nonlinear processes.
2.2 The boundary module
A thin cylindrical shell of arbitrary uniform resistivity is modelled at the plasma edge, separated from a coaxial ideal conductor by a tuneable-width vacuum region. In our notation, we call $r_{\text {BC}}\geq a$
the numerical boundary radius (always located at the thin shell position), and $b>r_{\text {BC}}$
the outer ideal wall radial position.
The magnetic boundary leverages the consolidated thin-shell relations (see e.g. Gimblett Reference Gimblett1986), descending from the assumption that the radial-$\boldsymbol {B}$
is continuous across the wall, and prescribing
where $\boldsymbol {B}_{r,\text {BC}}=\boldsymbol {B}_r\vert _{r=r_{\text {BC}}}$
, $\tau _w=\mu _0 r_{\text {BC}} \delta _w/\eta _w$
is the resistive time scale of the magnetic field diffusion through the shell (with $\delta _w$
shell thickness and $\eta _w$
shell resistivity), $\boldsymbol {E}_w$
is the 3-D electric field inside the wall, $E_z^{0,0}$
is the uniform electric field representative of the inductive loop voltage and $\boldsymbol {B}_t=B_\theta \hat {\boldsymbol {e }}_\theta+B_z\hat {\boldsymbol {e}}_z$
is the tangential magnetic field, $\hat{\boldsymbol{e}}_\theta$
and $\hat{\boldsymbol{e}}_z$
being the azimuthal and axial unit vectors. Here, square brackets indicate the difference of the enclosed quantity between just outside the wall (in vacuum) and just inside (plasma edge); the vacuum solution is obtained in SPECYL from the analytical solution of Poisson's problem, and from the self-consistent value of $B_{r,\text {BC}}$
, as thoroughly described in Spinicci et al. (Reference Spinicci, Bonfiglio, Chacón, Cappello and Veranda2023).
Equation (2.7) provides a boundary condition for $\boldsymbol {B}_{t}$
by imposing the continuity of the tangential electric field between wall surface and plasma edge, via Ohm's law,
where the subscript ‘bc’ indicates that quantities are evaluated at numerical boundary.
This is a general set-up, capable of reproducing diverse experimental conditions, from an ideal wall attached to the plasma (when $\tau _w\gg \tau _{\text {sim}}$
, $\tau _{\text {sim}}$
being the simulation duration), to a generic plasma-facing shell of arbitrary resistivity. In addition, the proximity of the outer ideal conductor can freely vary, from attached to the plasma to a virtually infinite distance.
If the magnetic diffusion through the resistive wall is instantaneous with respect to the time scales of plasma dynamics ($\tau _w\ll \tau _{\text {dyn}}$
), the wall acts as if transparent to the magnetic field. In this regime, a fully self-consistent boundary module must be capable of reproducing free-boundary current-driven instabilities, including the external kink modes. Good quantitative accuracy is expected, as long as SPECYL's cylindrical geometry can realistically describe the evolution of a large-aspect ratio plasma, including the linear growth of the mode, possibly up to low-amplitude saturating kink modes (see Eriksson & Wahlberg Reference Eriksson and Wahlberg1997).
A 3-D velocity at the plasma edge, including a free NF, is mandatory for numerical self-consistency of resistive-wall boundary conditions. Moreover, allowing finite flow into the wall is required for a hot fusion plasma, which must be set free to impinge on the wall and annihilate (see Zakharov et al. Reference Zakharov, Galkin and Gerasimov2012).
Two sets of velocity BCs hereinafter described are sketched in figure 1; a schematic view of the poloidal cross-section illustrates the plasma (pink), the thin resistive shell at $r=r_{\text {BC}}$
(claret) and the outer ideal wall (grey). For the case of SPECYL's free NF formulation (already described in Spinicci et al. (2023)), $r_{\text {BC}}=a$
, and a dotted black contour suggests fully consistent modulation of the velocity boundary in reproducing an $m=2$
external kink mode. The fixed NF implementation enforces Dirichlet boundary conditions on velocity fluctuations and is mostly stable to the free-boundary modes, unless enforcing a suitable PV region (shaded in lighter pink) to displace the plasma interface away from the numerical boundary.
Figure 1. Schematic representation of the two boundary formulations. The plasma is represented in pink, surrounded by a thin resistive shell (claret) at $r=r_{\text {BC}}$
and by an outer ideal wall at $r=b>r_{\text {BC}}$
(in grey). In the free flow formulation (a), the plasma edge corresponds to the numerical boundary ($a=r_{\text {BC}}$
), and a 3-D velocity is permitted: a black dotted contour suggests self-consistent edge flow modulation compliant with an $m=2$
external mode. Fixed flow boundary conditions are illustrated in (b): in this case, the numerical boundary is separated from the plasma interface by a PV region (in a lighter pink shade), and no flow perturbation is present at $r=r_\text {BC}$
.
The free NF velocity boundary of SPECYL aims at an $\boldsymbol {E}_w\times \boldsymbol {B}$
drift at the plasma edge, thus purely perpendicular to the magnetic boundary,
A finite parallel velocity could be enforced by means of, e.g. a Bohm speed at the plasma boundary, even if its value would still be irrelevant to self-consistency with (2.6)–(2.8), as they only depend on $\boldsymbol {v}_{\perp }$
.
For the purpose of this comparative study, we will also consider an approximated fixed normal flow formulation, enforcing Dirichlet boundary conditions for all non-axisymmetric fluctuations, on the top of a purely one-dimensional (radial and axisymmetric) pinch flow velocity (analogous to the one in Bonfiglio et al. (Reference Bonfiglio, Chacón and Cappello2010)), defined as
where $\hat {\boldsymbol {e}}_r$
is the radial unit vector. Equation (2.10) describes the leading term of (2.9) in the standard RFP operation regime, while it is typically small in tokamak configuration; for the sake of this study, (2.10) represents almost a Dirichlet BC also for the axisymmetric harmonic of NF.
The theoretical figures of merit used throughout this work reproduce well-known results from the linear stability theory of external kink modes, and are derived leveraging the energy principle, briefly discussed in appendix A (see also Freidberg Reference Freidberg2014). The resulting linear differential equation for the plasma displacement $\boldsymbol {\xi }$
is numerically solved, unless an analytical solution can be attained (as for the case reported in Appendix B), for any given initial axisymmetric equilibrium utilising the semianalytical LENS code (i.e. linear Euler and Newcomb solver, from the names of the two main differential equations). The LENS is implemented in the IDL (interactive data language) coding language and has been widely tested, both against analytical figures of merit and manufactured solutions (see Spinicci Reference Spinicci2023). An overview of LENS is reported in Appendix C.
To comply with SPECYL's implementation, we will enforce the so-called straight tokamak limit theory (corresponding to a very large plasma aspect ratio), assuming a cylindrical plasma with a circular cross-section.
3 Numerical set-up
The external kink mode stability will be studied for the two velocity boundary formulations already reported in figure 1, and for two alternative initial equilibria, each defined by a specific initial current density distribution,
We will refer to the case of $\lambda =0$
as the flat current case, and to the $\lambda =1$
case as the Wesson current case, since this is the typical equilibrium employed in Wesson's theoretical works (see, e.g. Wesson Reference Wesson1978).
Throughout this study we will leverage an initial axisymmetric Ohmic equilibrium (as in Delzanno, Chacón & Finn Reference Delzanno, Chacón and Finn2008), defined as the solution of the force-free condition $\boldsymbol {J}\times \boldsymbol {B}=0$
and the parallel Ohm's law $\boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}=\eta \boldsymbol {J}\boldsymbol {\cdot }\boldsymbol {B}$
. The use of these equations is meant to provide initial equilibrium profiles that are approximate solutions of the full system of equations solved by SPECYL, apart from the small convection and viscous terms in the momentum balance equation. The Ohmic equilibrium equations, as well as (2.1)–(2.5), assume mass-density and a resistivity only depending on $r$
,
The $\eta _{\max }$
in (3.2) is a truncation value preventing numerically demanding edge values of $\eta$
. The two profiles in (3.2) and (3.3), along with the uniform scalar viscosity, are set as a specific input to each simulation. A uniform loop voltage in the plasma volume yields the initial equilibrium current density, roughly proportional to plasma conductivity $1/\eta$
; the parameters $A, B$
and $C$
are chosen accordingly, to reproduce the desired initial current profiles as in (3.1). As for the plasma density, the model parameters $\alpha$
and $\rho _\text {BC}$
are chosen so to approximate a sharp-step-like profile across the plasma interface; for the case of the free flow formulation, where the interface is at numerical boundary, we adopt a uniform density profile ($\alpha =0,\ \rho _\text {BC}=1$
).
Figure 2 illustrates the hereinafter considered Ohmic equilibria. Figure 2(a,b) correspond to the two alternative choices of BCs: the free and fixed NF formulations, respectively. Figures 2(a i) and 2(b i) correspond to the equilibria with a flat current density ($\lambda =0$
), while figures 2(a ii) and 2(b ii) represent the Wesson-like initial current distribution ($\lambda =1$
). For each of the four combinations, two diagrams are displayed: in the upper one the SPECYL's profiles of current density and safety factor (in black) are compared with their respective reference profiles (in pink and teal, respectively), whereas in the lower diagram the two corresponding input profiles of resistivity and mass density are also represented (solid lines), along with the theoretical sharp-step-like density profile (orange dashed).
Figure 2. Initial equilibria employed in this study. The panels correspond to different velocity boundary formulations: (a) free NF and (b) fixed NF with PV (shaded). Subpanels (i) and (ii) report equilibria corresponding to the two alternative initial current profiles enforced. For each of the four combinations two diagrams are reported: the superior one illustrates SPECYL's initial current density and safety factor (in black), in comparison with the analytical model (coloured dashed lines), whereas the lower one displays the two equilibrium-relevant input profiles of SPECYL's density and resistivity (solid lines). The theoretical (infinitely sharp) mass density profile is also reported as an orange colour dashed line. Unlike elsewhere, the radial coordinate on horizontal axes is normalised to the plasma radius $a$
for better graphical clarity.
As anticipated, the fully consistent free normal flow implementation allows full identification of the plasma boundary with the numerical boundary ($a=r_{\text {BC}}$
), whereas the simplified fixed NF boundary requires that we displace the plasma interface away from the numerical boundary ($a< r_{\text {BC}}$
) in order to get it linearly unstable: the resulting high-resistivity, low-density PV region is shaded.
All along, the displayed level of confidence in reproducing the theoretical reference profiles appears to be widely acceptable, as the only visible discrepancy – the smooth transition of the current density profile across $r=a$
for the flat current model with a fixed NF – is insufficient to produce sensible discrepancies, e.g. in the safety factor profile.
Table 1 presents the main simulation parameters, corresponding to the initial equilibria represented in figure 2, always enforced in the following unless otherwise specified.
Table 1. Numerical set-up for the initial equilibria represented in figure 2. For each one we report the coefficients defining the input profiles of plasma resistivity ($A,\ B,\ C, \eta _{\max }/\eta _0$
) and of plasma density ($\alpha,\ \rho _\text {BC}/\rho _0$
), the number of radial points ($N_r$
), yielding comparable resolution in the plasma region, depending on its width ($a/r_{\text {BC}}$
).
The edge current in Wesson's model, as much as both $J_z^{0,0}$
and $\rho$
in the PV region, should be null; this is numerically unfeasible in SPECYL. Instead, table 1 reports the most extreme values for $\eta _{\max }/\eta _0,\ \alpha \ \text {and}\ \rho _\text {BC}/\rho _0$
still compliant with numerical stability of the code time stepping scheme. To recreate more reliably the numerically challenging Wesson-like equilibrium, we evolve the axisymmetric initial state for a short time (approximately spanning the initial $1 \tau _A$
) with magnified viscosity and resistivity ($\nu =\eta _0=10^{-2}$
, $\eta _{\max }=5$
) to substantially damp unwanted residual flow and current in the PV region, before starting the non-axisymmetric dynamics.
The plasma region width ($a/r_{\text {BC}}$
) for the fixed flow case balances two opposing and quantitatively relevant requirements. On the one hand, as will be manifest in § 4, the Dirichlet condition for the velocity fluctuations, enforced by (2.10), produces a non-physical boundary layer in which the internal plasma flow unnaturally adapts to meet the imposed value in $r_{\text {BC}}$
; this takes approximately the outermost $10\,\%$
of the numerical domain, and should be included in the PV region. On the other hand, since the plasma density and current in the PV cannot be exactly null, the PV width should remain significantly thinner than the plasma region width, to provide a negligible contribution to the integral mass and current (further amplified by the cylindrical Jacobian). Also, in equally spaced meshes like in SPECYL, a finite PV width proportionally enlarges the computational cost, by subtracting radial resolution from the plasma core region. For both initial equilibria, the optimal value of $a/r_{\text {BC}}$
had to be determined a posteriori, maximising the agreement with the theoretical expectations; in both cases, this coincides with $a=0.88 r_{\text {BC}}$
, as reported in table 1. The number $N_r$
of radial mesh points is fixed so to have approximately the same resolution (${\sim }900$
points) for the plasma region with both fluid boundary formulations.
The dynamical evolution of such equilibria is here described with a two-dimensional (2-D) spectrum of modes, composed by the axisymmetric one, with the addition of several non-axisymmetric harmonics with a fixed helical pitch, such that $m=2n$
. This configuration is adequate to study the linear stability to the external kink $m=2, n=1$
, already considering just four non-axisymmetric modes with $m>0$
, plus their conjugates. The initial destabilising contribution, driving the system towards relaxation, is applied in the form of a velocity perturbation,
where $\mathcal {A}=10^{-6}$
. In all simulations we operate with a very fast resistive wall time scale, $\tau _w/\tau _A=10^{-6}$
, in compliance with the transparent wall assumption, and enforce a confidently large aspect ratio $R/a=20$
, as prescribed for the straight tokamak limit.
Uniform plasma viscosity, $\nu =M^{-1}$
, and on-axis resistivity, $\eta _0=S_0^{-1}$
, are kept mutually identical throughout this study (with the exception of Appendix D). Their specific value will be reported for each case study, individually. In all cases, the enforced values should tend to ideal plasma conditions (very low dissipation); this, along with the challenging equilibrium configurations already discussed and the inherently fast dynamics of external kink modes, requires very small time steps
Table 2 summarises the hitherto described input parameters underlying the dynamical evolution in our simulations. For some of them (namely, $S_0,~M\ \text {and}~b/a$
) a range of enforced values is provided; for the ideal wall proximity, $b/a$
, we have underlined the most commonly adopted value.
Table 2. Numerical set-up for the dynamical evolution of SPECYL's simulations in this study. We report the adopted values for the 2-D helical pitch ($m/n$
), the number of modes with positive $m$
($N_{m>0}$
), the aspect ratio ($R/a$
), the ideal wall proximity to the plasma edge ($b/a$
), whose typical value is underlined, the transparent shell time scale ($\tau _w/\tau _A$
), the on-axis Lundquist number and uniform viscous Lundquist number ($S_0$
and $M$
) and the simulations time step ($\Delta t/\tau _A$
).
4 Numerical results
We present in this section SPECYL's predictions of linear instability to the external kink mode 2/1 of the two equilibria introduced in § 3, and leveraging both sets of velocity BCs. It will be apparent that not only the free NF velocity boundary provides a better match for the analytical theory, but it also achieves more reliable asymptotic convergence to the expected outcomes.
4.1 The flat current model
The flat current model is widely present in numerical verification studies (e.g. in Hoelzl et al. Reference Hoelzl, Huijsmans, Merkel, Atanasiu, Lackner, Nardon, Aleynikova, Liu, Strumberger, McAdams, Chapman and Fil2014; Ferraro et al. Reference Ferraro, Jardin, Lao, Shephard and Zhang2016; Marx & Lütjens Reference Marx and Lütjens2017), owing to its reduced complexity and to the existence of an exact analytical solution of the Energy principle (see Shafranov Reference Shafranov1970; Wesson Reference Wesson1978). In fact, for a 2/1 external kink mode in the straight tokamak limit, the linear eigenfunction is $\boldsymbol {\xi }^{2,1}(r)=\xi _{r,a}^{2,1} (\hat {\boldsymbol {e}}_r+\text {i}\,\hat {\boldsymbol {e}}_\theta )r/a$
and the dispersion relation reads
where the characteristic time scale is the toroidal Alfvénic time $\tau _A/\varepsilon =R/v_A$
, $\varepsilon =a/R$
being the inverse aspect ratio. A complete analytical derivation is available in Appendix B. For the confidently large aspect ratio enforced and for $b/a\gg 1$
, the characteristic dynamical time scale $\tau _{\text {dyn}}=1/\gamma$
predicted by (4.1) ranges in the tens of $\tau _A$
; hence, our choice of $\tau _w/\tau _A=10^{-6}$
reliably fulfils the transparent-wall condition.
Figure 3 reports a convergence study to the theoretical growth rate from ideal MHD theory for a specific value of the uniform safety factor ($q(r)=q_a=1.5$
) and considering various levels of plasma dissipation. As SPECYL plasma becomes increasingly ideal, we see that the fully consistent free NF boundary produces an asymptotic convergence to the expectation. A significantly different trend can instead be observed for the fixed flow boundary with PV, which non-monotonically oscillates around the target value for rather dissipative SPECYL plasmas, crossing it before $S_0=M=10^6$
and asymptotically plummeting towards stability as the plasma gets increasingly ideal (no visible mode growth can be detected in our simulations above $S_0=M=10^{12}$
). Indeed, a very ideal plasma corresponds to a very ineffective PV, owing to the numerical limitations on $\eta _{\max }/\eta _0$
that keeps vacuum resistivity unphysically finite.
Figure 3. Comparison of linear growth rates for the flat current case study for $q_0=q_a=1.5$
and for various values of $S_0=M$
. Free NF BCs asymptotically reach the ideal MHD theoretical value (dotted horizontal line) for vanishing plasma dissipation. The fixed NF BCs cross the target value around $S_0=M=10^{6}$
and decrease asymptotically towards stability (which is almost reached above $S_0=M=10^{12}$
).
Figure 4 extends the growth rate analysis to a wider range of initial flat current equilibria, for various values of $q(a)$
. The additional normalisation factor $q(a)/\varepsilon$
for $\gamma$
is widely present in the literature (e.g. in Wesson Reference Wesson1978) and it is therefore employed here as well. For the free NF formulation we employ $S_0=M=10^7$
(but similar results could be displayed for $10^8$
or higher), whereas for fixed flow BCs the dissipation level that results in the best agreement with linear theory has been chosen, corresponding to the target-crossing simulation (i.e. $S_0=M=10^6$
; see figure 3). For the latter set of BCs, slightly more ideal simulations with $S_0=M=10^7$
are also shown in green, for comparison. Again, the theoretical dispersion relation (see (4.1)) is reported as a black line. Linear growth rates from SPECYL are numerically estimated with a logarithmic fit, except for $q(a)\leq 1$
and $q(a)\geq 2$
where a fast-oscillating low-amplitude signal is observed, lacking any detectable global trend.
Figure 4. Comparison of linear growth rates for the flat current case study, for various values of $q(a)$
. The theoretical dispersion relation is represented as a black contour. (a) The free NF BCs (blue) stick robustly to the theoretical expectation; (b) the fixed NF BCs with PV are mostly consistent with the ideal MHD theory in the target-crossing set-up (red), besides slightly lower precision for $q(a)\leq 1.4$
and $\gamma <0$
in the ideally stable domain (due to sensible dissipation levels). Such an overall positive result is, however, frail, as choosing a slightly more ideal set-up (green) visibly spoils the previous agreement.
The simulations performed with a 3-D fluid boundary are in robust agreement with theoretical values (within $1\,\%$
). The SPECYL simulations performed with a fixed velocity boundary and PV at the target-crossing dissipation level achieve a mostly comparable agreement with theory in the mode instability window with $1< q(a)<2$
, despite somewhat lower precision on the left-hand edge of the hill. The effect of sensible plasma dissipation is, however, manifest in some simulations with $q\gtrsim 2$
, significantly departing from the ideal mode stability condition $\gamma =0$
with a negative damping rate. Reducing the plasma dissipation level, the ideal stability condition is recovered, at the cost of a much weakened agreement with the theoretical expectation in the unstable region.
The radial velocity and magnetic field eigenfunctions for the three set-ups employed in figure 4 are benchmarked with the ideal MHD predictions, in figure 5, along with the radial profile of the axial current oscillations. The radial velocity is confronted with the radial component of the linear plasma displacement: $\xi _{r}^{2,1}/\xi _{r,*}^{2,1}=v_{r}^{2,1}/v_{r,*}^{2,1}$
(being $v_{r,*}=v_r\vert _{r=r_*}$
, and $r_*=0.9a$
an arbitrarily chosen position for normalisation). Once again, the reliability of the free NF is remarkable. The same consideration applies to the fixed flow BCs for what concerns the plasma domain. Yet its behaviour in the PV is manifestly ill-behaved, in two main aspects: there is finite (and considerable) flow in the vacuum region, and the radial magnetic field pitch variation across the plasma edge reveals exaggerated surface currents. The strong deviation of SPECYL magnetic field in the PV from the vacuum linear eigenfunction (see (B9)) deteriorates with decreasing plasma dissipation and it is both sustained by the finite velocity, through a dynamo contribution, and by finite conductivity effects, allowing stray current fluctuations $\boldsymbol {J}^{2,1}=\boldsymbol {\nabla }\times \boldsymbol {B}^{2,1}$
.
Figure 5. Benchmark of eigenfunctions profiles between SPECYL simulations with different velocity BCs and the analytical theory (black, dashed curves), for the flat current equilibrium. (i,ii) The SPECYL 2/1 radial velocity and magnetic field are compared with the radial displacement $\xi _r^{2,1}$
and the radial magnetic eigenfunction of ideal MHD: (a) the free NF formulation reproduces remarkably the theoretical curves; (b) the target-crossing simulation with the fixed NF BCs and (c) the slightly more ideal simulation with the same BCs are mostly reliable in the plasma region, with significant discrepancies in the PV. All profiles are normalised at $r_*=0.9a$
, indicated by vertical dotted lines. (iii) We also report the (axial) current density oscillations for a unitary magnetic perturbation, compared with the analytical $J_z^{2,1}(r)\approx J_{\text {surf }}^{2,1}\delta (r-a)$
: exceeding currents across plasma interface sustain exaggerated pitch variation of $B_r^{2,1}$
for the fixed flow formulation.
Figure 5(a/,iii,b/,iii,c/,iii) reports the dominant (axial) contribution to the current density fluctuations: for each subpanel, the linear current density is normalised to the amplitude of the magnetic perturbation. The analytical expectation for this case study is nearly proportional to a Dirac's delta distribution, centred across the plasma interface: $J_z^{2,1}(r)\approx J_{\text {surf }}^{2,1}\delta (r-a)$
(as discussed in Appendix B), which is quite well reproduced by the free NF implementation. The fixed NF formulation instead presents a much broader current distribution across plasma edge, much of which is on the PV side. Despite the current peak being lower than for the free NF case, the total surface current is larger, motivating the exaggerated pitch variation of $B_r^{2,1}$
across interface, as we already commented.
Up to this point the outer ideal wall in SPECYL's BCs formulations has always been kept largely distant from the plasma edge ($b/a=100$
; see table 2). In figure 6 we explore the stabilising effect of bringing the ideal shell closer to the plasma edge, for two fixed flat current equilibria with $q(a)\approx 1.1$
and $q(a)\approx 1.5$
(qualitatively equivalent results would be obtained for any other value of $1< q(a)<2$
). The set-up with dissipation values that mostly agree with the analytical theory ($S_0=10^{7}$
and $S_0=10^{6}$
for the free and fixed NF boundaries, respectively), are benchmarked with (4.1) (black line) showing remarkable qualitative agreement; until the ideal wall is as far as $b=3a$
, its effect is mostly negligible and all numerical datasets correctly saturate to the far-wall limit growth rates already discussed for figures 3 and 4. As the ideal wall closes in, SPECYL's modes 2/1 progressively damp, up to the stability threshold proximity ($b/a\approx 1.78$
for $q_0=1.1$
, $b/a\approx 1.19$
for $q_0=1.5$
), which is accurately reproduced by both sets of BCs. As the ideal wall approaches further, the time traces of SPECYL's modes become fast oscillating noisy signals without any measurable global trend.
Figure 6. Scan on the outer ideal wall proximity $b/a$
, within the flat current model and for two different equilibria ($q(0)=1.1$
in (a), $q(0)=1.5$
in (b)): ideal wall stabilisation of MHD activity becomes increasingly effective as it closes in to the plasma edge (shaded). The best performing set-ups for both the free and the fixed NF boundary are compared with the theoretical prescription (black). Two vertical dotted lines highlight the tightest possible proximity with either approach: owing to finite PV width, the fixed NF formulation cannot reproduce a closer proximity than $r_{\text {BC}}/a\approx 1.136$
.
As for the quantitative agreement, the reliability of the free NF BCs lies within $6\,\%$
(ever decreasing towards $1\,\%$
as the ideal wall moves apart), whereas the level of compatibility of the fixed NF boundary lies within $20\,\%$
for $q_0=1.1$
and within $6\,\%$
for $q_0=1.5$
. We highlight with two vertical dotted lines the tightest achievable proximity with either set of BCs, as $b=r_{\text {BC}}$
; the presence of a PV region in the case of the fixed NF prevents simulating an ideal wall closer that $b\approx 1.136a$
. For the present case study this inaccessible parametric region presents reduced interest, but would prevent simulating a portion of the unstable region for any flat current equilibrium having $q(a)\geq 1.6$
, unless picking a larger, less favourable, $a/r_{\text {BC}}$
ratio, or using free NF conditions. Furthermore, this limitation could in principle restrict the ability of the fixed flow boundary to reproduce, e.g. real experimental conditions when a narrow vacuum interval separates the plasma edge from the first wall.
The displayed level of agreement for the fixed NF formulation in figure 6(a) is in line with what previously obtained with the JOREK-STARWALL code enforcing an analogous set-up for the same case study of a flat current, at $q(0)=1.1$
(see McAdams Reference McAdams2014). Substantially improved agreement with the analytical theory was achieved in Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016) with the M3D-C$^{1}$
code, by adopting a resistivity profile strongly discontinuous across the plasma interface; such results can be reproduced also in SPECYL, however, finding reduced reliability as the resonance radius approaches the plasma interface, as documented in Appendix D.
4.2 Wesson's model
Wesson's current model presents no known analytical solution for its eigenfunction $\boldsymbol {\xi }$
and eigenvalue $\gamma$
. Alongside, its self-consistent modelling requires a numerically cumbersome fast-ramping resistivity profile at the plasma edge, corresponding to a vanishing current density at plasma interface, as already visible in figure 2. These are the reasons why it finds almost no application in verification studies of MHD codes.
Conversely, while certainly more realistic than a uniform current channel surrounded by pure vacuum, its numerical complexity makes it an important and challenging test for SPECYL's new BCs.
Also for this case study, the chosen numerical set-up enforces $\tau _w/\tau _A=10^{-6}\ll 1/\gamma$
, and $R/a=20$
. Theoretical expectations are obtained through the LENS code.
Figure 7 presents an analogous convergence study to the one presented in figure 3; the free and fixed NF BCs with PV (in blue and red, respectively) are compared with the expected theoretical growth rate (black dotted line), for a specific value of $q(a)$
and at various levels of plasma dissipation. The qualitative trends closely resemble what is already observed for the flat current equilibrium; only a 3-D and fully consistent free NF boundary can robustly reproduce the expected ideal MHD dynamical behaviour in the limit of low plasma dissipation, whereas the fixed NF presents just a (this time narrow) crossing of the target growth rate and becomes asymptotically stable as the PV becomes too inviscid and conductive to suppress MHD activity.
Figure 7. Convergence study to the theoretical growth rate predicted by the LENS code (black dotted line) of SPECYL's simulations with free and fixed NF boundaries, respectively, employing the Wesson-like equilibrium with $q(a)\approx 1.8$
, for various levels of plasma dissipation.
The quantitative aspect is yet rather different from figure 3, since even the asymptotic convergence of the free NF is manifestly less effective, reaching a $17\,\%$
discrepancy only around $S_0=M=10^{10}$
. The rationale for this is to be sought in the effect of strong and fast-ramping resistivity profiles at plasma edge. In fact, combining (2.3) and (2.4) from SPECYL's MHD model we get
The first term on the right-hand side of (4.2) is the only one prescribed by the ideal MHD and corresponds to the dynamo action in the magnetic field variation. Of the other two terms, the one proportional to the resistivity gradient was vastly negligible in the plasma region in the previous case study with uniform $\eta$
, as it still is now in the plasma core, where $\eta$
is mostly uniform. Conversely, at the plasma edge, both for the free and for the fixed NF BCs, this same term is now roughly one to two orders of magnitude larger than the other resistive term, for any assigned value of $S_0$
.
The two velocity boundaries produce comparable predictions at high plasma dissipation. In fact, in such a regime the very high resistivity produces the effect of a PV region at the plasma edge even for the free NF case, strongly damping the edge velocity and thus making the two models quantitatively alike.
Looking at the fixed flow boundary simulations in figures 3 and 7 we notice a pretty similar qualitative trend, where numerical growth rates exceed the theoretical one at high dissipation and converge to stability for an ideal plasma. This is, however, not the case for the 3-D flow boundary simulations, that asymptotically converge to the expected growth rate in the ideal MHD limit, but they do it from below in figure 3 and from above in figure 7. This effect can be linked to the above mentioned tendency of large edge resistivity gradients (only present in Wesson's current modelling) to behave as PV. Notably, the PV has the function of displacing the effective plasma radius to a more internal position with respect to $r_{\text {BC}}$
. The most direct consequence of this is that the effective edge safety factor must be evaluated at a more internal position and will thus be lower than the nominal value of $1.8$
; as we will see in figure 8, $q(a)$
values slightly below $1.8$
are associated with larger growth rates according to the specific dispersion relation of Wesson's model.
Figure 8. The linear growth rates of several Wesson-like equilibria, corresponding to diverse values of $q(a)$
, are illustrated. We compare the theoretical expectation obtained by the LENS code with (a) SPECYL's simulations with the free NF boundary (blue), at a very reduced plasma dissipation, and (b) SPECYL's simulations with fixed NF, corresponding to the target-crossing dissipation level (red), as discussed in figure 7, and for a mildly more ideal plasma set-up (green). In this challenging numerical test, the free NF formulation retains quantitative and qualitative reliability, unlike the fixed NF formulation.
In figure 8 we extend the verification benchmark to a range of initial Wesson-like equilibria, each defined by a different value of $q(a)$
. The theoretical dispersion relation is obtained through the LENS code and is represented as a black contour. The least dissipative set-up of figure 7 is leveraged in SPECYL's simulations with free NF BCs (blue circles), whereas for the fixed NF BCs with PV the target-crossing dissipation level is represented in red, while a slightly more ideal set-up is reported in green. Also in this quite challenging case, the SPECYL code with its full 3-D velocity boundary displays a remarkable capability of reproducing the correct dispersion relation. Quite at the opposite, the fixed NF implementation seems drastically incapable of achieving even qualitative agreement with the theoretical expectations; not only are single growth rates strongly incorrect (by up to $190\,\%$
for $q(a)<2$
for the best performing set-up), but the instability window is completely inaccurate, predicting as significantly unstable equilibria that should instead be completely stable to the external kink modes. Internal resistive tearing modes for $q(a)>2$
cannot explain such behaviour, as they should present drastically lower growth rates for the given set-ups ($\gamma \tau _A/\varepsilon \sim 10^{-4}$
, as estimated by LENS). Instead, the high resistivity at the plasma edge is once again blurring the interface between plasma and vacuum, as will appear evident from figure 9; the choice of a more conductive set-up (as represented by the green circles) importantly reduces the overstretching of the instability window, however, at the cost of further spoiling the agreement with theory for $q(a)<2$
. Interestingly, SPECYL's numerical results with a fixed NF boundary are very sensitive to the specific value of $a/r_{\text {BC}}$
: slight increase or decrease of such parameter significantly extend the unstable region beyond $q(a)=2$
and may even produce unphysical kink stability where $q(a)\lesssim 1.6$
. In our perspective, such a strong dependence on a set-up parameter that retains no physical meaning would in itself already drastically undermine the reliability of fixed NF BCs in modelling free-boundary instabilities for the large class of equilibria with a vanishing edge plasma current.
Figure 9. Comparison of radial profiles of $v_r^{2,1}$
(i), $B_r^{2,1}$
(ii) and $J_z^{2,1}$
(iii), for free NF (in blue) and fixed NF BCs with PV (target-crossing case in red). A dashed black line represents the theoretical profiles of ${\xi }_r^{2,1}/{\xi }_{r,*}^{2,1}$
, $B_r^{2,1}/B_{r,*}^{2,1}$
and $J_z^{2,1}/J_{z,*}^{2,1}$
($r_*$
being marked by vertical dotted lines), as computed by the LENS solver. The large resistivity gradient at plasma edge produces deviations of SPECYL's eigenfunctions from the theoretical model, especially visible for the flow and the current density profiles. For the fixed NF case, the PV behaves as a higher-resistivity plasma rather than a vacuum: $v_r^{2,1}$
presents the vertical asymptote at resonance ($r/a\approx 1.054$
, indicated by a cyan dotted line) typical of tearing modes, while a variation in the pitch of $B_r^{2,1}$
reveals finite surface currents $\boldsymbol {J}^{2,1}$
at the same resonance layer, also clearly visible in (b iii).
A more complete picture can be drawn by looking at the fluid and magnetic radial eigenfunctions in figure 9. Here, SPECYL's outcomes with the free and fixed NF boundary formulations, in blue and red, respectively, are compared with the radial displacement $\xi _r^{2,1}$
and the radial magnetic fluctuation $B_r^{2,1}$
, as obtained from the LENS code. Despite the good agreement with theoretical radial profiles of both formulations within the plasma core region, high resistivity effects can be seen in all plots at the plasma edge, in the form of a mild reduction of edge radial velocity and magnetic field, correlated to an excess of edge current oscillations. It should be noted that the good performance of fixed NF is in this case, as well as in figure 5, strongly cherry-picked; we already commented on how frail the displayed agreement is upon small variations in physical conditions ($q(a),\ S,\ M$
) or in $a/r_{\text {BC}}$
.
Looking at the PV, shaded in figure 9(b), we may observe some apparent marks of resistive plasma dynamics in a region that we would rather regard as vacuum. In particular, $v_r^{2,1}$
presents a sharp change of sign across $r_0\approx 1.054a$
, paired with a pitch variation of $B_r^{2,1}$
across the same radial position, revealing the presence of a (mostly axial, due to straight tokamak ordering) finite slab current $\boldsymbol {J}^{2,1}$
, well visible in figure 9(b iii). This is the signature of a resistive internal (tearing-like) mode across a resonance layer: $q(r_0)\approx 2$
.
5 Final remarks and future perspectives
New resistive wall BCs with a fully consistent 3-D velocity boundary were recently implemented in the nonlinear MHD codes SPECYL and PIXIE3D. Such a module consists of a thin shell of tuneable resistivity at plasma edge, separated from an outer ideal conductor by an arbitrary-width vacuum region. The magnetic boundary enforced at the plasma edge leverages the resistive Ohm's law across the thin shell and is self-consistently coupled with an $\boldsymbol {E}\times \boldsymbol {B}$
velocity field, dictated by the internal MHD and by Ohm's law at boundary.
In this paper we presented a challenging verification study against ideal MHD external kink modes theory. The conceptually simplest set-up to perform this study consists in modelling a free interface at plasma boundary, possibly also including the effect of an ideally conductive wall at finite distance. In practice, this is achieved by setting SPECYL's thin shell resistivity so high that it becomes transparent for the magnetic field dynamics; such a numerical set-up should be appropriate for the study of linear free-boundary MHD instabilities driven by plasma current.
With the fully consistent free NF formulation recently implemented in the SPECYL code, the linear evolution of external kink modes can be properly reproduced already looking at the transparent shell as a free interface in itself. This is, however, not the case with a not fully consistent set of velocity BCs, imposing null edge velocity perturbations. As the latter is mostly stable with a rigid wall attached to the plasma boundary, we leveraged a high-resistivity, low-density PV region to displace the plasma edge from the numerical boundary; such an approach is routinely employed in most analogous studies already present in the literature (see, e.g. Merkel & Strumberger Reference Merkel and Strumberger2015; Marx & Lütjens Reference Marx and Lütjens2017).
The study considered the instability to external kink modes of two classes of initial tokamak equilibria characterised by circular cross-section and by a very large aspect ratio. The first one is defined by a uniform current channel, coincident with the plasma itself, and the second one features a Wesson-like parabolic current density that vanishes at plasma edge.
The first case study presents a well-known analytical solution for the linear instability eigenfunctions and growth rates; along with its reduced numerical complexity, this is why it is used in most numerical verification studies already present in the literature. This case study could be remarkably reproduced on a wide range of equilibrium current intensities and of external ideal wall positions by both sets of boundary conditions, each for a selected numerical set-up. However, the free NF formulation achieved robust and asymptotic convergence to the theoretical growth rates and eigenfunctions in the correct limit of an increasingly ideal plasma set-up, whereas the fixed NF formulation best fitting set-up was found for intermediate plasma viscoresistive values, while showing asymptotic stability in the ideal MHD limit. We identify the unavoidable numerical limitations in the modelling of the PV region as the most likely candidates to explain such behaviour.
The second case study presents no known analytical solution and needs enforcing a specific semianalytical solver; along with its increased numerical complexity, this makes its employment in this work rather exceptional, with respect to past analogous verification studies of nonlinear MHD codes. Also in this case, a robust and asymptotic convergence of SPECYL's free NF simulations was illustrated, even if reached at considerably lower dissipation levels than in the previous case, owing to increased numerical challenges. On the other hand, the fixed NF formulation proved itself extremely unreliable even in predicting bare modes stability. Growth rates analysis for a fixed equilibrium current intensity found the fixed NF formulation in SPECYL to be substantially incompatible with the expectation, both in the resistive plasma limit and in the ideal plasma limit (where it predicts once more full modes stability). Particularly concerning appears the strong and non-monotonic dependence of simulations outcomes on set-up parameters retaining no physical meaning, namely the distance of the numerical boundary from the plasma interface.
The main findings of this work are the following.
(i) A free NF formulation is essential for complete self-consistency of resistive wall BCs; indeed, it allows reliable modelling of the linear growth of current-driven external modes identifying the plasma interface with the numerical boundary. As far as we know, this result is unprecedented for a nonlinear MHD code.
(ii) The enforcement of a PV region at plasma edge is unavoidable in presence of BCs bearing a velocity formulation in partial consistency with the rest of the leveraged model. Nonetheless, its usage has general applicability in predictive MHD modelling, and it allows to describe the nonlinear evolution of external instabilities (especially in more advanced codes where $\rho$
and $\eta$ can evolve through time); its presence can indeed be justified to account for a scrape-off layer in realistic experimental conditions. Our study has, however, highlighted reasons for concern when its enforcement aims to account for an actual vacuum region, especially in the face of physical configurations that are challenging from the numerical point of view. In particular, the more or less narrow parametric region, where the theoretical models are matched with acceptable fidelity, appears to be hardly predictable a priori (unlike the asymptotic convergence to the right theoretical value, displayed by the fixed NF formulation). Even if a suitable set-up for fixed NF with PV formulations was eventually found in optimal agreement with Wesson's model predictions (possibly, extending the approach explored in Appendix D), an a posteriori determination of such a set-up might not be completely satisfying for a numerical tool that aims at performing predictive modelling.(iii) Albeit the specific set-up featuring a transparent wall, coincident with the plasma interface, has reduced applicability in predictive simulations, we hope to provide with it a remarkably challenging and repeatable test for the self-consistency of free NF BCs in MHD nonlinear codes. The SPECYL's documented reliability on this critical set-up is expected to hold in the general case of a finitely conductive wall, or even considering a cold and sparse scrape-off layer around the plasma core.
Along with the already published verification study against PIXIE3D, this paper completes the numerical verification process of SPECYL's resistive wall module. The code will now be employed in both first-principle studies and predictive numerical modelling in the tokamak and the RFP configurations. For the RFP, preliminary numerical studies performed with the SPECYL code enforcing a simplified formulation of its boundary module with a fixed NF velocity like in (2.10) are reported in Marrelli et al. (Reference Marrelli, Martin, Puiatti, Sarff, Champan, Drake, Escande and Masamune2021). Analogous modelling work is now being carried out with the free NF boundary formulation, to extend the validity of previously achieved results and in view of anticipating the next device RFX-mod2 (Marrelli et al. Reference Marrelli, Cavazzana, Bonfiglio, Gobbin, Marchiori, Peruzzo, Puiatti, Spizzo, Voltolina and Zanca2019; Terranova et al. Reference Terranova2024). A central point of interest for the RFP community concerns the self-organised helical plasma states at high current regimes, with periodicity $m=1,\ n=7$
in the RFX-mod device (see Bonfiglio et al. Reference Bonfiglio, Veranda, Cappello, Escande and Chacón2013; Marrelli et al. Reference Marrelli, Martin, Puiatti, Sarff, Champan, Drake, Escande and Masamune2021), which still lack self-consistent numerical predictability. Stimulated by the promising results presented in this paper we are confident in SPECYL's resistive wall module future findings also in this merit.
Free NF BCs are expected to play a determinant role especially in plasma-limiter configurations, namely including RFP routine operation and tokamak disruptions. The impact of a fully consistent fluid-magnetic formulation could be of especial interest in the modelling of externally imposed resonant magnetic perturbations at the plasma edge, widely used in numerical simulations for the modelling of disparate physical conditions, e.g. including edge localised modes triggering (Cathey et al. Reference Cathey, Hoelzl, Lackner, Huijsmans, Dunne, Wolfrum, Pamela, Orain and Günter2020), hybrid tokamak scenarios (Piovesan et al. Reference Piovesan, Bonfiglio, Cianciosa, Luce, Taylor, Terranova, Turco, Wilcox, Wingen, Cappello, Chrystal, Escande, Holcomb, Marrelli, Paz-Soldan, Piron, Predebon and Zaniol2017) and stimulated QSH states in RFP configuration (Veranda et al. Reference Veranda, Bonfiglio, Cappello, Escande, Auriemma, Borgogno, Chacón, Fassina, Franz, Gobbin, Grasso and Puiatti2017).
First-principle studies are also being pursued for the numerical validation of Zakharov's theory on vertical displacement events in the tokamak configuration. Such a model predicts the existence of external kink-like instabilities flowing into the wall and sustaining surface currents in its passive structures; the combination of these currents with the plasma-confining magnetic field would be at the origin of experimentally observed intense horizontal loads on device support structures (see Zakharov et al. Reference Zakharov, Galkin and Gerasimov2012). Conceptually similar simulations to the ones presented in this work, but considering more realistic wall resistivities, will prove if a fully consistent velocity boundary with a normal velocity perpendicular to the plasma edge is indeed sufficient to produce the so-called wall-touching kink modes and sustain quantitatively comparable currents with respect to the model predictions.
Acknowledgements
The authors want to especially acknowledge the two anonymous referees for raising interesting points, which prompted further important investigation, particularly regarding different approaches for the modelling of the PV region, now reported in Appendix D. Further important acknowledgements go to L. Chacón for actively contributing in the formulation of SPECYL's BCs, P. Zanca for valuable feedback on the manuscript and to M. Hölzl and F.J. Artola for inspiring conversations.
Editor P. Helander thanks the referees for their advice in evaluating this article.
Funding
The present work has been entirely conducted at Consorzio RFX (Padova, Italy) and was partially cofunded by the University of Padova. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (grant agreement no. 101052200 – EUROfusion). 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.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Stability to the kink mode in the straight tokamak limit
The external kink mode is an ideal MHD, free-boundary instability, fruitfully described to linear order with a variational approach dubbed the ‘energy principle’ (e.g. Freidberg Reference Freidberg2014), allowing to find the most energetically favourable plasma shape relaxation, $\boldsymbol {\xi }=\sum _{m,n}\boldsymbol {\xi }^{m,n}(r)\,{\rm e}^{\text {i} m\theta -\text {i} nz/R-\text {i} \omega t}$
, and its linear growth rate $\gamma ={\rm Im}\{\omega \}$
. This corresponds to an extremant of the cost-functional $W(\boldsymbol {\xi })$
, representing the conservative work produced in the deformation against the restoring force,
where the integral is intended over the volume of the whole system $V_{\text {tot}}=V_{\text {plasma}}+V_{\text {vacuum}}$
, and the star indicates complex conjugation. The problem eigenvalue is retrieved from the linearised energy balance equation,
with $\boldsymbol {\xi }$
and $W$
as in (A1). Depending on the sign of $W$
, external kink modes growth rates can either be null (in case of stability) or strictly positive (instability).
In the straight tokamak limit theory, obtained for an asymptotically large aspect ratio, pressure gradients become negligible and the geometry reduces to a cylinder. In such a regime the initial axisymmetric equilibrium and its stability properties under non-axisymmetric perturbations solely depend on the equilibrium axial current density profile $J_{z}^{0,0}$
and on the mass density profile $\rho$
. The assumption of a circular cross-section is further enforced in this study, to achieve full compliance with SPECYL's implementation.
The variational problem in (A1) is solved within a Lagrangian approach, yielding the Euler equation for $\boldsymbol {\xi }$
. In the straight tokamak limit and with the reasonable assumption of an incompressible displacement, $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\xi }=0$
, this reads
with suitable scalar functions $f,~g$
and $\mathcal {C}$
, as follows:
where $B_{\textrm {eq},z}$
is the axial component of the equilibrium field. Substituting (A3) in (A4) and integrating by parts, the minimum work is simply expressed as
which can be employed in (A2) to yield the mode stability analysis.
Equivalent to Euler's equation is solving the so called Newcomb's equation for the magnetic flux perturbation $\psi ^{m,n}=-\textrm {i} rB_r^{m,n}$
, which is related to the radial plasma displacement by the frozen-in law of ideal MHD (see, e.g. Freidberg Reference Freidberg2014),
with $k_n=n/R$
, and with $B_{\textrm {eq},\theta }$
and $B_{\textrm {eq},z}$
cylindrical components of the equilibrium magnetic field. Newcomb's equation is usually derived in a force-free plasma, from the assumption that the linearised magnetic torque produced by the perturbation is negligible (see, e.g. the appendix in Fitzpatrick (Reference Fitzpatrick1999)),
where $H=m^2+ k_n^2 r^2$
and
Equations (A3) and (A8) are completely equivalent for an asymptotically large aspect ratio (see, e.g. § 3.7 in Spinicci (Reference Spinicci2023)).
The Newcomb's equation, combined with Faraday's law and with the solenoidal property of $\boldsymbol {B}$
, yields useful relations for the spatial components of magnetic field and current density linear fluctuations,
Appendix B. Analytical solutions to the flat current model
The flat current model is characterised by a uniform distribution of axial current and mass density with a sharp step profile at plasma edge (see Shafranov Reference Shafranov1970),
For this simplified configuration, Euler's equation reduces to
whose solution is, via regularity conditions on axis, $\xi _r^{m,n}=\xi _{r,a}^{m,n}(r/a)^{m-1}$
. Enforcing now the incompressibility assumption for vanishing $\varepsilon$
(i.e. in the straight tokamak limit), the whole displacement eigenvector results as
Correspondingly, the associated minimum conservative work produced by the plasma deformation is
while the kinetic term in the linearised energy balance (A2) reads
The resulting dispersion relation is hence
where the toroidal Alfvénic time $\tau _A/\varepsilon$
naturally emerges from the analytical computation.
To derive the magnetic fluctuation eigenfunction, it is convenient to solve Newcomb's equation, which also greatly simplifies. Retaining only the leading order in the small parameter $\varepsilon$
we get in both the plasma and the vacuum region
whose solution is, via regularity conditions, $\psi ^{m,n}(r)=\psi _a^{m,n}(r/a)^{m}$
inside plasma and $\psi ^{m,n}(r)=\psi _a^{m,n}(a/r)^{m}$
outside, so that the radial magnetic eigenfunction results from (A10) as
where $B_{r,a}^{m,n}=\textrm {i}\psi _a^{m,n}/a$
.
In § 4, alongside $\xi _r^{2,1}$
and $B_r^{2,1}$
, we also discuss of current density oscillations. These can be obtained from (A11)–(A15). It is, however, immediately manifest, without making any calculation, that the leading terms in equations (A13)–(A15) are those proportional to $\sigma '$
(the prime here means radial derivative). In fact, from (A9) it appears that this term contains the derivative of $J^{0,0}_{z}$
, which in turn is proportional to a Dirac's delta centred at the plasma interface.
From this consideration we can argue that the radial current oscillations will be mostly negligible with respect to the other two spatial components. These two will thus approximately read
Finally, the axial component should largely dominate over the azimuthal, since $B_{\textrm {eq},z}\gg B_{\textrm {eq},\theta }$
via the tokamak ordering, as long as $q_a=B_{\textrm {eq},z} \varepsilon /B_{\textrm {eq},\theta }\gtrsim 1$
.
Appendix C. The LENS code and numerical solutions of Wesson's model
In the general case, (A3) and (A8) present no known analytical solution. A semianalytical solver is hence mandatory to predict the linear dynamics of MHD perturbations.
This appendix aims at presenting an overview of our such code, dubbed LENS (i.e. linear Euler and Newcomb solver), which is implemented in the IDL coding language and has been thoroughly verified against analytical figures of merit and manufactured solutions (see Spinicci Reference Spinicci2023). The LENS code solves (A3) and (A8), along with (A10)–(A15) and $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\xi }^{m,n}=0$
, to determine the radial profiles of linear perturbations in $\boldsymbol {\xi }^{m,n}$
and $\boldsymbol {B}^{m,n}$
. The modes stability analysis for the external kink mode is obtained from (A2) and (A6) in the straight tokamak limit, while suitable dispersion relations are also provided for tearing modes stability, as in Wesson & Campbell (Reference Wesson and Campbell2011), Furth, Rutherford & Selberg (Reference Furth, Rutherford and Selberg1973), Bondeson & Sobel (Reference Bondeson and Sobel1984) and Porcelli (Reference Porcelli1987). A wide variety of built in graphical routines help visualise the obtained results.
We consider in this study a class of initial equilibria (Wesson & Campbell Reference Wesson and Campbell2011), defined by initial current density profiles as in (3.1) (reported as follows for clarity),
of which, only the simplified case with $\lambda =0$
is amenable of analytical solution, as already discussed in Appendix B.
Figure 10 illustrates the modes stability analysis and the 2-D profiles of the perturbation for the external kink mode instability, produced by LENS for a cylindrical plasma with a circular cross-section, surrounded by pure vacuum. In figure 10(a), the dispersion relation is solved numerically for several values of $q(a)$
, considering diverse current-peaking factors ($\lambda$
) in (C1). In general, the instability domain of an $m/n$
mode corresponds to $m-1\leq nq(a)\leq m$
. It is, however, apparent the stabilising effect of finite magnetic shear, positively scaling with $\lambda$
; this becomes increasingly effective as we move to high-$m$
modes, while the $m=1$
mode ($\forall n$
) is completely insensitive of $\lambda$
. For the flat current case, here solved numerically in optimal agreement with (B6), in the absence of magnetic shear, growth rates can only scale with $J_{z}^{0,0}(0)$
; such a dependency is, however, eliminated by the chosen normalisation $\gamma q(a)\tau _A/\varepsilon$
, since $q(a)\propto 1/J_{z}^{0,0}(0)$
. Analogous figures in the literature perfectly compare with what described here.
Figure 10. Example of LENS's graphical elaboration: (a) stability analysis to the external kink mode of several initial equilibria as in (C1), corresponding to various values of $q(a)$
and $\lambda$
; (b) poloidal projections of fluid and magnetic linear eigenfunctions in the case of a flat equilibrium current. Along with linear eigenfunctions, for each case the external resonance radius is indicated with a green dashed contour, and the axial surface currents codirected and counterdirected with respect to $J_{z}^{0,0}$
are indicated in blue and red, respectively. The illustrated predictions are perfectly in line with analogous figures on the past literature and with the analytical predictions described in Appendix A.
In figure 10(b), below each corresponding instability domain, the projection onto the plasma poloidal cross-section (shaded in pink) of the displacement and magnetic eigenfunctions is represented for the flat current profile. Formally similar figures would be obtained for different values of $\lambda$
. Looking at displacement patterns, we see for each mode the characteristic incompressible deformation of corresponding $m$
-periodicity; such modes are destabilised by an external resonance near plasma interface (green dashed contour). The magnetic field fluctuations are offset by ${\rm \pi} /2m$
with respect to $\boldsymbol {\xi }^{m,n}$
, in compliance with the frozen-in law: $B_r^{m,n}=\textrm {i} F(r)\xi _r^{m,n}/r$
(see (A7) and (A10)). Their vorticity across the plasma interface is sustained by intense and predominantly axial surface currents, whose sign is discordant with the equilibrium current (which enters the page) on the side of the growing instability (blue), and concordant where the interface moves inwards (red), in agreement with the analytical theory (see, e.g. Wesson Reference Wesson1978; Zakharov et al. Reference Zakharov, Galkin and Gerasimov2012).
Appendix D. Rendering of PV with a discontinuous resistivity profile
Remarkable benchmark agreement with the linear theory of external kink modes has been achieved in the past with the M3D-C$^{1}$
code (see Ferraro et al. Reference Ferraro, Jardin, Lao, Shephard and Zhang2016), enforcing Dirichlet boundary conditions for the velocity (fixed NF formulation) and PV. The study is conducted against the flat current model, and employs a strongly discontinuous resistivity profile across the plasma interface.
The same approach has been pursued also in SPECYL: the corresponding initial equilibrium is displayed in figure 11. The resulting current density profile indeed more closely resembles the theoretical current channel, with respect to the alternative approach leveraged for this case study in the rest of our work, as it is significantly more damped in the PV region. The ratio of PV to core resistivity is constrained by numerical stability; in this appendix we focus on equilibria where its value is $10^5$
.
Figure 11. Initial flat current equilibrium, employed in Appendix D. (a) The SPECYL's unperturbed current density and safety factor (in black) are compared with the corresponding theoretical profiles (pink and teal, respectively). (b) The two equilibrium-relevant input profiles of density (black) and resistivity (red) are reported, along with the theoretical infinitely sharp mass-density profile (orange). The strongly discontinuous resistivity profile as in Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016) is able to produce an almost vanishing equilibrium current in the PV region.
The optimal benchmark agreement obtained in the figure 2 of Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016) is successfully reproduced with SPECYL in figure 12, using the same physical parameters ($q(a)\approx 1.1$
, $a= 0.88r_{\textrm {BC}}$
). The best performing dissipation level has been chosen to be $S_0=10^8$
, $M=2\times 10^7$
, and is represented in the plot with full red dots. A new set of simulations has been performed with the free NF formulation, enforcing the same dissipation level ($S_0=M=10^8$
) and is reported in blue. Red empty contours mark the positions of fixed NF simulations already reported in figure 6, while a black line illustrates the theoretical expectation (see (B6)). The new set-up is here found in remarkable agreement with the model, especially when $b\geq 3a$
, where it performs comparably, with respect to the free NF simulations. A slight shift in the critical ideal wall proximity is visible, as the mode stability seems to be achieved around $b/a\approx 1.74$
instead of $1.78$
. This is likely related to a slight upwards shift in the $q(r)$
profile inside the PV region with respect to the corresponding analytical profile, as already visible in figure 11, as if the effective $q_a$
was slightly higher; such an effect was previously prevented when enforcing (3.2), as some residual $J_z$
in the PV region compensated for the missing core current due to the unavoidably smooth current decay at the plasma edge. Even this way, compatibility lies within $2\,\%$
for most the plot.
Figure 12. The SPECYL's simulations with a flat current and $q(0)=1.1$
, for various proximities of the ideal wall to the plasma edge. The analytical growth rate is reported in black. The SPECYL's simulations with the most favourable plasma dissipation set-up and discontinuous resistivity across plasma edge (full red dots) lie in remarkable agreement with the expectation, apart from a slight mismatch in the critical proximity. Fixed NF simulations already presented in figure 6 are also reported as empty red circles, while free NF simulations with the same plasma resistivity $S_0=M=10^8$
are reported in blue. This figure successfully reproduces figure 2 in Ferraro et al. (Reference Ferraro, Jardin, Lao, Shephard and Zhang2016).
In figure 13 we extend the analysis just performed to a full scan on the value of $q(a)$
for the flat current model, keeping the ideal wall at large distance from the plasma interface. As we did before, along with the new set of SPECYL's simulations (full red dots), we also report the free NF simulations performed with the same dissipation level (blue dots), and the corresponding values obtained with a continuously varying $S(r)$
across plasma edge (empty red circles). It can be seen that the discontinuous resistivity profile is successful in reproducing the correct stability boundaries, also obtaining significantly more accurate predictions for $0< q(a)<1.5$
than the other fixed NF set-up. The quality of the agreement is, however, spoiled as the magnetic resonance moves closer to the plasma edge, where the resistivity presents a numerically demanding radial gradient.
Figure 13. Comparison of the linear theory growth rates for several values of $q(a)$
for the flat current model (black line), with the best-performing set-ups for the fixed NF formulation, with a continuous resistivity across plasma interface (red empty circles, as in figure 4b) and with a discontinuous resistivity (full red dots). The latter reliably captures the correct stability boundaries and the expected growth rates as $q(a)\gtrsim 1$
, but gradually detaches from the model predictions as the magnetic resonance approximates the plasma interface.
A great variability of the dynamical behaviour as $q(a)\lesssim 2$
is found in SPECYL's simulations on the choice of $S_0$
and $M$
, but also on small variations of the interface radius $a/r_{\textrm {BC}}$
. This seems to suggest that a wisely tailored set-up should be viable for SPECYL to reproduce the linear expectations within reasonable fidelity, as obtained with the XTOR-2F code in Marx & Lütjens (Reference Marx and Lütjens2017) with a tightly fitting ideal wall, enforcing an analogous set-up. Such a configuration would, however, resent of the same general limitations of the PV approach, already discussed in the rest of this paper; reduced robustness (and predictability) of the most suitable set-up, and cumbersome physical relevance of parameters that should retain no physical meaning, like $a/r_{\textrm {BC}}$
.
Open access
