Nonlinear MHD simulations of external kinks in quasi-axisymmetric stellarators using an axisymmetric external rotational transform approximation
Rohan Ramasamy, Matthias Hoelzl, Erika Strumberger, Karl Lackner, Sibylle Günter
NNonlinear MHD simulations of external kinks inquasi-axisymmetric stellarators using anaxisymmetric external rotational transformapproximation
R Ramasamy , , M Hoelzl , E Strumberger and K Lackner and S G¨unter Max-Planck Institut f¨ur Plasmaphysik, Boltzmannstraße 2, 85748 Garching beiM¨unchen Max-Planck Princeton Center for Plasma Physics, Princeton, New Jersey 08544USAE-mail: [email protected]
Abstract.
Reduced magnetohydrodynamic (MHD) equations are used to study thenonlinear dynamics of external kinks in quasi-axisymmetric (QA) stellarators withvarying fractions of external rotational transform. The large bootstrap currentsassociated with high beta plasmas typically make QA configurations susceptible tolow n external modes, limiting their operational space. The violence of the nonlineardynamics, and, in particular, when these modes lead to a disruption, is not yetunderstood. In this paper, the nonlinear phase of external kinks in an unstable QAconfiguration with an edge safety factor below two is simulated.The nonlinear MHD code JOREK is used, validating the linear dynamics againstCASTOR3D. The external rotational transform is approximated axisymmetrically, andvaried from the tokamak limit to understand its influence on the nonlinear dynamics.The violence of the kink instability is quantified, and shown to reduce with theincreasing external field. At the same time, the external kink triggers (9, 5) and(7, 4) internal modes that exacerbate the loss in confinement during the nonlinearphase, such that it remains large over much of the parameter space. It is only with asignificant fraction of external rotational transform that these subsequent modes arestabilised.
Keywords : quasi-axisymmetry, stellarator, nonlinear magnetohydrodynamics, externalkink modes
Submitted to:
Nucl. Fusion a r X i v : . [ phy s i c s . p l a s m - ph ] F e b
1. Introduction
In recent years, there have been ongoing efforts to produce a compact, quasi-axisymmetric (QA) stellarator that can be scaled up to a fusion reactor [1]. Animportant question that remains unanswered in this work is the susceptibility of suchdevices to disruptions.Previous studies of 3D equilibria with Ohmic current profiles have shown thatstellarators can be considered vertically stable, and free from tearing mode drivendisruptions when a sufficiently large external rotational transform is applied [2, 3, 4].However the potential instabilities for QA configurations are different, and at the timeof writing, no experimental QA devices have been constructed to provide insight intotheir nonlinear stability properties. For this reason, the current understanding of theseconfigurations relies on existing theory and numerical codes.Linear ideal magnetohydrodynamic (MHD) codes have been applied in theoptimisation of QA stellarators [5, 6]. When these devices make use of a significantbootstrap current to maintain the equilibrium conditions, they are typically susceptibleto similar MHD instabilities as ITER advanced operational scenarios [7]. Such candidateinstabilities include infernal modes, double tearing modes, and external kinks. Theseinstabilities will be milder than those in tokamaks, which do not have the stabilisinginfluence of an external rotational transform, magnetic well, and 3D shaping, but theviolence of the nonlinear dynamics that they can drive is still unclear.The linear, viscoresistive MHD code, CASTOR3D [8], was recently applied to aQA configuration with an edge safety factor below two [9]. As such, the stellaratorwas known to be unstable to violent kink instabilities. CASTOR3D provides a detailedanalysis of these instabilities, characterising the mode structure, and the influence ofseveral stabilising effects, such as resistive walls. While these studies give a detaileddescription of the linear dynamics, there are remaining questions regarding the nonlineardynamics. In particular, it is unclear at what amplitude the observed kink instabilitieswill saturate, and whether they can trigger internal modes, leading to a disruption, ashas been observed in kink driven disruptions in tokamaks [10].Nonlinear, stellarator capable, MHD codes are necessary to answer these questionswith a high fidelity. There are presently only a few nonlinear codes that can be appliedto stellarator geometries, such as MEGA [11] and FLUXO [12]. The nonlinear MHDcode, JOREK [13], is currently being extended to model stellarators. A hierarchy ofstellarator capable models have been derived [14, 15], and are currently being validated.Assessments of the nonlinear dynamics of global current driven instabilities require anevolution of the magnetic field in time. For advanced stellarators, where the magneticfield description requires many toroidal harmonics, simulations are expected to becomputationally demanding.As such, it follows to consider supplementing the high fidelity models beingdeveloped in JOREK with simplified simulations of the nonlinear dynamics ofstellarators. A common simplification to choose is to modify configurations such thatthey can be modeled axisymmetrically. Axisymmetric approximations of stellaratorshave a long history, beginning with the first stellarator equilibrium expansion [16]. Theordering used in this expansion is similar to the high beta tokamak expansion, with theinclusion of a large helical field to represent the external rotational transform, leading toa modified Grad Shafranov equation for the equilibrium [17]. Soon after the derivation ofthis expansion, it was applied to linear MHD stability [18]. In the 1970s, common toolsfor assessing tokamak stability, such as the tearing mode equation, and reduced MHDmodels were adapted to look at stellarator instabilities, both linearly and nonlinearly[19, 20, 21]. Such models have provided insight into MHD phenomena experimentallyobserved on W7-A and W7-AS, as well as being applied to predict the stability propertiesfor NCSX [2, 3, 22].This paper outlines an implementation of the axisymmetric stellarator approxima-tion in JOREK, called the virtual current model. The approach makes use of currentequilibrium codes to provide an appropriate description of the external rotational trans-form. This is achieved by representing the externally generated field by a virtual current inside the plasma that is parallel to the total equilibrium magnetic field. Past studiesusing a similar approach have attributed this additional current to a lower hybrid cur-rent drive [22]. Herein, we prefer to use the term virtual current to emphasise that thiscurrent is not a true plasma current, but a model to incorporate the effect of 3D coilson the rotational transform axisymmetrically.The virtual current model in JOREK has already been validated for simple tearingmode, and vertical displacement event test cases, showing the expected stabilisation ofboth instabilities [23]. Linear growth rates of double tearing modes and resistive kinksin Wendelstein 7-X (W7-X) have been compared against CASTOR3D, using differentaxisymmetric approximations. The results of this validation are shown in the Appendix.In this study, the virtual current model is applied to external kinks in a compact,high beta QA configuration. These instabilities prove to be more challenging for themodel to capture, due to the global nature of the modes, and the influence of 3Dshaping and the magnetic well, which cannot be captured axisymmetrically. It is shownthat there is reasonable agreement between the linear growth rates and eigenfunctionsobserved in the virtual current model, and CASTOR3D. Simulations are then continuedinto the nonlinear phase, varying the fraction of external rotational transform to see howthe dynamics are influenced by its presence. While the external kink is stabilised bythe external field, internal modes are triggered that lead to a loss of confinement acrossthe middle and outer region of the plasma. The results give a first indication of thepotential dynamics that might occur in QA configurations.The rest of the paper is arranged as follows. In Section 2 the implementation of thevirtual current model in JOREK is described. The equilibria considered for this studyare discussed in Section 3. The virtual current model is validated in the linear regimefor external kink cases in Section 4. The nonlinear dynamics of these external kinks isthen assessed in Section 5, and the paper concluded in Section 6.
2. The virtual current model in JOREK
JOREK is a fully implicit, 3D, nonlinear extended MHD code, capable of being run withreduced and full MHD models of varying complexity [13]. In cylindrical coordinates,the poloidal planes are discretised using bi-cubic C Bezier finite elements, with a realFourier series in the toroidal direction [24]. A Crank-Nicholson or Gears time steppingscheme is used to advance the physics variables, leading to a large sparse matrix thatis iteratively solved using GMRES [25]. The solution to the block diagonal matricesformed by individual toroidal harmonics is typically solved using PASTIX [26], orSTRUMPACK [27], and used as a pre-conditioner for the solution to the full system ofequations.Unlike early axisymmetric approximations which use an analytic representation ofthe vacuum poloidal flux, modern axisymmetric approximations of stellarators normallytake a numerical approach, which represents the external field as a toroidal current[9, 22, 28]. During the CASTOR3D assessment of an unstable QA configuration, anaxisymmetric approximation to the stellarator was found using the Variational MomentsEquilibrium Code (VMEC) [29]. The n = 0 Fourier components of the stellaratorboundary were used to generate an axisymmetric equilibrium with the same pressureand safety factor profiles. The poloidal field generated externally is then replaced byan increase in the toroidal current in the device. The enhanced destabilising currentdensity gradients lead to more unstable kink modes with similar structure. This typeof axisymmetric approximation is refered to as the tokamak approximation because noattempt is made to capture the external rotational transform of the original stellarator.In the following study, to replicate the external rotational transform, the additionalcurrent in the tokamak approximation is removed from the plasma, as virtual currents,leading to models refered to as virtual current approximations .In order to do this in JOREK, the solution to the Grad Shafranov equilibriumequation is modified, such that∆ ∗ ψ t = − µ R p (cid:48) ( ψ t ) − F F (cid:48) ( ψ t )= − µ R p (cid:48) ( ψ t ) − (cid:104) F p F (cid:48) p ( ψ t ) + F v F (cid:48) v ( ψ t ) (cid:105) . (1)Therefore, the solution to the Grad Shafranov equation does not need to change.Once the solution for ψ t is found, the current is split into plasma and virtualcurrent components, both of which remain functions of ψ t . As shown in Figure 1,the modification does not break the force balance of the equilibrium, because thediamagnetic current represented by F F (cid:48) is parallel to the total equilibrium magneticfield. The virtual current represented by F v F (cid:48) v is calculated from VMEC computations,in such a way as to preserve the radial profile of the plasma current density from thestellarator. This approximation is considered to be suitable for capturing similar linearstability properties because the linear stability depends on the local gradients of thecurrent density and magnetic shear. Figure 1.
Sketch of virtual currents in a cylindrical plasma. The virtual current isalong the total equilibrium magnetic field and as such does not modify the solution tothe Grad Shafranov equation.
During the time evolution, the virtual current is fixed in space and time. JOREK’ssingle fluid, reduced MHD equations [30] are used, removing the parallel momentumequation to improve numerical stability. The total poloidal flux, electric potential,density and temperature, are evolved in time, subject to the governing equations ∂ψ t ∂t = R [ ψ t , u ] + η ( j t − j v ) − F ∂u∂φ (2) ∂ρ∂t = −∇ · ( ρ v ) + ∇ · (cid:16) D ⊥ ∇ ⊥ ρ + D (cid:107) ∇ (cid:107) ρ (cid:17) + S ρ (3) ∂ ( ρT ) ∂t = − v · ∇ ( ρT ) − γρT ∇ · v + ∇ · (cid:16) κ ⊥ ∇ ⊥ T + κ (cid:107) ∇ (cid:107) T (cid:17) + 23 R η ( T ) ( j t − j v ) + S p (4) R ∇ · (cid:32) R ρ ∇ ⊥ ∂u∂t (cid:33) = 12 (cid:104) R |∇ ⊥ u | , R ρ (cid:105) + (cid:104) R ρω, u (cid:105) + [ ψ t , j t − j v ] − F R ∂j t ∂φ + (cid:104) ρT, R (cid:105) + Rµ ( T ) ∇ ω (5) j t = ∆ ∗ ψ t = j p + j v (6) ω = 1 R ∂∂R (cid:32)
R ∂u∂R (cid:33) + ∂ u∂Z (7)where the Poisson bracket, [ a, b ] = ∇ φ · ( ∇ a × ∇ b ), has been used, and the reducedMHD ansatz is assumed for the plasma velocity, v , and magnetic field, BB = ∇ ψ t × ∇ φ + F ∇ φ (8) v = − R F ∇ u × ∇ φ (9)where ψ t is the total poloidal flux generated by the total current, j t , from theplasma, j p , and virtual current, j v . F = RB φ represents the toroidal field, which isfixed in time. The necessary modifications for the virtual current model are markedin blue, and amount to the inclusion of a virtual current, which is subtracted fromthe terms responsible for releasing free magnetic energy. The change to equation 2 isequivalent to the addition of a current source, fixing the background virtual current.In equations 4 and 5, the virtual current is removed from the Ohmic heating term andLorentz force, as the plasma should only experience the influence of its own current.As external kinks are modeled in this study in the no wall limit, keeping the n = 0component fixed, the influence of virtual currents on the boundary conditions also needsto be considered. The free boundary condition for n > n · B JOREK = n · B STARWALL (10) n × B JOREK = n × B STARWALL (11) n · J JOREK = n · J STARWALL (12)on the JOREK simulation boundary. Because the virtual current is a function ofthe initial ψ t , and the simulation boundary is initially a flux surface, equation (12) isindependent of the virtual currents in the domain. In JOREK, the boundary conditionsare implemented in terms of the total magnetic field [33], and so there is no need tochange the coupling between the two codes. As shown in Section 3, it is importantto note that the correct rotational transform profile of a stellarator in the vacuumregion can only be achieved axisymmetrically by having virtual currents in the vacuumregion, to represent the external rotational transform. As the STARWALL couplingassumes that there is no current in the vacuum region outside the simulation domain,the rotational transform in this region will not match the modeled stellarator. Thislimitation is not considered to be significant, because the simulation domain is alreadyextended beyond the plasma boundary, as discussed in Section 3. Therefore, the correctexternal rotational transform is captured in the region of interest. (a) (b) (c)(d) (e) (f) Figure 2.
Equilibrium profiles for the l = 2 stellarator (a-c) and QA device (d-f)test cases, and their axisymmetric approximations. The normalised toroidal flux, ˆΦ,is used as a radial coordinate. The tokamak approximation (blue) has a larger totaltoroidal current, I t , and vacuum magnetic shear than the stellarator (black). Thevirtual current approximation (red) has the same pressure, and safety factor profile asthe stellarator. I t continues to increase in the vacuum region, beyond ˆΦ = 1. The changes outlined above are simple to implement, and can just as easily beapplied to more advanced models, allowing for a fast, albeit qualitative assessment ofstellarator nonlinear dynamics that are intended to be compared against more accuratesimulations including the full geometry, once such models are available.
3. Equilibria of interest
Two stellarator equilibria are modeled in this study. A simple l = 2 stellarator isconsidered in the linear validation, alongside the main QA stellarator that is continuedinto the nonlinear phase. Their profiles are shown in Figure 2, alongside the profiles oftheir corresponding axisymmetric approximations.The l = 2 stellarator was designed as a simple linear test case for external kinks,to validate the virtual current model in a high aspect ratio, low beta stellarator. A flatcurrent density profile is used, such that the case would be strongly unstable to externalkinks. The QA configuration is a slightly modified version of the two field periodequilibrium previously modeled in [9]. In particular, the pressure profile was modifiedto be parabolic, such that the gradients of density and temperature at the plasma edgesmoothly transition to the vacuum conditions. The final safety factor profile has also Figure 3.
Workflow for linear validation between JOREK and CASTOR3D. been modified such that the edge safety factor is closer to the q = 2 rational surface. Toobtain a reasonable temperature of approximately 10 keV , the number density is set to3 . × m − .The workflow used to construct the axisymmetric approximations for the linearvalidation in Section 4 is outlined in Figure 3. The modifications are similar to thoseused in [9], but require additional steps to generate axisymmetric equilibria with thesame rotational transform outside the plasma. This is important as the magnetic shearof the plasma has a stabilising effect on the external kink. For these preliminarystudies, the original fixed boundary computations are extended smoothly into thevacuum region, linearly scaling the edge toroidal flux with the new area of the extendeddomain. This provides an approximate description of the stellarator equilibrium, whichis then transformed to a new axisymmetric approximation in the same way as hasbeen previously described. Numerical difficulties were encountered when solutions ofextended equilibria into the vacuum region were used in CASTOR3D. For this reason,the 3D equilibria are only modeled in CASTOR3D up to the plasma boundary, couplingthe vacuum and MHD equations at this interface.Comparing the extended equilibrium stellarator profiles with their tokamak andvirtual current approximations, as shown in Figure 2, the pressure profile is consistentbetween the three cases. The safety factor profile is the same for the stellarator andits approximations within the plasma, except at the plasma edge. In this region, itcan be seen that the magnetic shear in the tokamak case is larger than the stellaratorand virtual current approximation. The virtual current approximation matches thestellarator profile by having virtual current in the vacuum region to represent the fieldgenerated from external coils.A comparison of the flux surfaces of the generated equilibria for poloidal crosssections at φ = 0 . π, . π and 0 . π are shown in Figure 4. This shows that afterthe extension of the plasma, the QA equilibrium has not changed considerably from[9]. In addition the n = 0 Fourier harmonics of the QA equilibrium are plotted againstcorresponding surfaces from the axisymmetric approximation. Here, it can be seen thatthe Shafranov shift of the axisymmetric approximation is larger than the 3D case. For (a) (b) (c) (d)(e) (f) (g) (h) Figure 4.
Flux surface contours for the l = 2 stellarator (a-c) and QA configuration(e-g) test cases. Surfaces of the stellarator equilibria (black) and their virtualcurrent approximations (red) are overlayed at different toroidal planes. The n = 0Fourier components of the l = 2 stellarator (d) are aligned with its axisymmetricapproximation. A similar comparison for the QA configuration (h) shows a deviationin the Shafranov shift. the l = 2 stellarator, the pressure has been set to an artificially low value, removingthis discrepancy. This shows one of the limitations of the axisymmetric approximation,which is that the magnetic well of a stellarator cannot be captured. As a result, theplasma can shift further towards the low field side of the device, modifying the magneticshear, and local current density gradient, which can have an effect on the MHD stability.
4. Linear validation
The results of the linear growth rate comparison are shown in Figure 5. As in previousstudies, the mode number of the instabilities from CASTOR3D corresponds to thedominant toroidal mode number of the respective eigenfunction, and the results aresplit into odd and even mode families. The sinusoidal and cosinusoidal pairs for each n of the stellarator instabilities are also plotted. These orthogonal solutions are onlynon-degenerate at low toroidal mode numbers in the QA case. For the l = 2 stellarator,with a simpler 3D geometry, the two mode families do not diverge significantly.External kinks are challenging to validate because of the sensitivity of these0modes to the edge q , and the differing boundary conditions between JOREK andCASTOR3D. JOREK applies a Dirichlet boundary condition on the normal plasmavelocity component, such that the vacuum region needs to be modeled to observeexternal kinks correctly, without an artificial influence from the boundary conditions.CASTOR3D assumes an infinitely resistive, perfect vacuum outside the simulationboundary. To get a good match between the codes, this demanding vacuum condition,which is harder to match than realistic configurations, needs to be captured in JOREK.Comparing the tokamak approximation in JOREK and CASTOR3D, the externalkink growth rates show very good agreement for both the l = 2 and QA configuration.In the QA case, the dip in the growth rates close to n = 5 and n = 6 is attributed tomode coupling between ( m ±
1, n) pairs inside and outside of the plasma near ˆΦ = 1that provide a stabilising effect. The (9, 5) and (11, 6) modes are very close to theplasma edge. Small deviations in the q profile, and resistivity profile, which defines theplasma boundary, can lead to larger deviations between the two codes. In the l = 2stellarator case, the flat q profile was intentionally chosen to be above the closest lowersidebands of toroidal modes up to n = 7, such that the stabilising effect is reduced. Ingeneral, the comparison of the tokamak approximations show that the reduced MHDmodel used in this study is appropriate for modeling these external kinks.Comparing the virtual current model with the stellarator results from CASTOR3D,we see that the general trend of the instabilities is captured in both cases. Particularlyin the QA configuration, the deviations are largest where the dominant (m, n) mode isstabilised by its ( m ±
1, n) poloidal sidebands. In particular, the n = 4 mode found inCASTOR3D for the QA configuration was not observed in JOREK. The reason for themodification of the mode structure cannot be identified from the growth rates alone,but is expected to be due to the modification of the toroidally averaged flux surfacesbetween the two cases, as shown in Figure 4.The eigenfunctions for the radial displacement in the QA case are shown in Figure6. The results in JOREK and CASTOR3D are generated using Boozer coordinates[34]. The radial structure of the instabilities is remarkably similar for the tokamakapproximation. Comparing the virtual current model with 3D results, the structure ofthe captured contributions from the leading toroidal mode number of the instabilitycompare very well. The grey dashed lines that grow in amplitude towards the plasmaedge are from toroidal sidebands of the stellarator instability observed in CASTOR3D.An axisymmetric model does not capture the toroidal mode coupling of a stellarator,because the equilibrium magnetic field is not a function of the toroidal coordinate [35].Given the complexity of the case being modeled, the linear mode structure ofthe virtual current approximation, and the true stellarator agree surprisingly well. Inparticular, the growth rate of the n = 1 instability is captured. As it will be shown inthe next section, the n = 1 mode is dominant in the nonlinear dynamics. The virtualcurrent model is therefore considered a reasonable tool for a first study of these modesto provide an understanding of the nonlinear dynamics.1 (a)(b) Figure 5.
Linear growth rates of external kinks observed in the l = 2 (a) and QA(b) stellarators. To compare the virtual current model, with the stellarator modes, theaverage of the sinusoidal (triangles) and cosinusoidal (squares) orthogonal solutions ofthe odd (red) and even (blue) mode families are used.
5. Nonlinear Dynamics of External Kinks
In this Section, the tokamak and virtual current approximations are continued intothe nonlinear phase. The initial stellarator considered for this study, has a relativelysmall external rotational transform, such that the plasma current, I p , of its virtualcurrent approximation remains high at 8 . M A . For this reason, several additionalsimulations were run at reduced plasma currents corresponding to an increased externalrotational transform. The same equilibrium is used, reducing I p while maintaining thesame normalised plasma current density profile as the original stellarator. In such a2 (a) (b) (c)(d) (e) (f) Figure 6.
Radial displacement eigenfunction comparisons of JOREK (solid) withCASTOR3D (dashed) in Boozer coordinates. The tokamak approximations (a-c)show very good agreement. The virtual current model (d-f) instabilities have similarcontributions from the leading toroidal harmonic of the equivalent QA instability inCASTOR3D. way, the full parameter space can be explored to give an intuition of when the externalkink becomes nonlinearly stable for this case.Table 1 shows the numerical and diffusive parameters used in the following nonlinearsimulations. The simulations carried out in the previous section were set up to matchthe boundary conditions of CASTOR3D, modeling a highly resistive vacuum regionoutside the plasma. For nonlinear simulations, these conditions are not realistic andnumerically demanding, and so some numerical parameters were modified. A Spitzer-like profile is used for resistivity, viscosity and their equivalent hyper-diffusive terms, inorder to dissipate sub-grid resolution current sheets that form during the penetration ofthe kink into the hot core. A realistic, temperature dependent, parallel heat conductivityis used. The parallel conductivity is set based on the Spitzer-Haerm conductivity, suchthat χ (cid:107) = 2 . × m s − at the magnetic axis. Because the simulation domain does notinclude a limiter, or X-point, there are closed flux surfaces outside the plasma, which donot replicate the transport physics of the vacuum region expected in a real experiment.For this reason, a relatively large perpendicular diffusion is used in the vacuum region3 Table 1.
Diffusive parameters used in nonlinear external kink simulations. The ψ dependent profiles sharply transition from the lower to upper limit of their ranges atthe plasma boundary.Parameter Range Dependence χ (cid:107) [ m s − ] 68 . − . × Spitzer-Haerm χ ⊥ [ m s − ] 1 . − . × Step function D (cid:107) [ m s − ] 847 . D ⊥ [ m s − ] 0 . − . η [Ω m ] 1 . × − − . × − Spitzer-like η num . × − − . × − Spitzer-like µ [ kg m − s − ] 9 . × − − . × − Spitzer-like µ num . × − − . × − Spitzer-like
Table 2.
Growth rates and saturation times for different plasma currents.Case γ n =1 t sat [ ms ] t offset [ms]Tokamak 4 . × I p = 8 . M A . × I p = 6 . M A . × I p = 5 . M A . × I p = 4 . M A . × I p = 2 . M A - - - to approximate the physics of a limiter around the plasma where transport out of thesimulation domain should be high. In the plasma region, the perpendicular particle andthermal diffusion coefficients are spatially constant. Because the parallel momentumequation is neglected in this study for simplicity and to improve numerical stability, aparallel particle diffusivity is used as a proxy for the parallel mass transport that wouldotherwise be neglected. Heat and particle sources are used to maintain the density andtemperature profiles for the timescale of the simulations.Six cases have been simulated at plasma currents ranging from 11 . M A , in thetokamak approximation, to 8 . M A , representing the original QA configuration, andfour conceptual cases with lower current at 6 .
75, 5 .
25, 4 .
00 and 2 . M A respectively.The instabilities are modeled with a polar grid, discretised using 150 poloidal and 120radial elements. Free boundary conditions are applied to all toroidal harmonics, exceptfor the n = 0 component. A toroidal resolution scan considering the evolution of themagnetic energies in the tokamak approximation is shown in Figure 7. The n = 1toroidal harmonic which leads the instability is sufficiently resolved using 10 toroidalharmonics.The observed linear growth rates, and times to saturation are shown in Table2. With the modifications described above, the structure of the growth rates in thelinear validation is damped, but the external kink modes remain violently unstable.The growth rates decrease significantly with increasing virtual current, until for the4 Figure 7.
Toroidal resolution scan for the tokamak approximation, since it has themost violent dynamics. The n = 1 mode, which dominates the instability, is sufficientlyconverged using 10 toroidal harmonics. I p = 2 . M A case, the virtual current approximation becomes linearly stable to externalkinks. Some of the plots in this section use a normalised time, ˆ t , allowing all cases to becompared, despite their varying time scales. ˆ t is normalised by the linear growth rateof the n = 1 kink instability, γ n =1 , such thatˆ t = 1 . × − γ n =1 ( t − t offset ) (13)where the prefactor is the time normalisation to JOREK units, and t offset is usedso that, at t = 0, all simulated cases have the same magnetic energy in the n = 1 mode,which has become linearly unstable.The Ohmic decay of the plasma current, η ( j t − j v ) in equation 2, is applied onlyafter the initial kink instability has begun to saturate, in order to preserve the currentprofile during the linear phase. A more realistic bootstrap current profile cannot beevolved from the initial conditions because the initial pressure and bootstrap currentdensity profiles have not been prescribed self-consistently. Allowing the current densityprofile to decay exacerbates resistive instabilities, and as a result can be considereda worst case scenario for the dynamics, as enhanced resistivity is likely to lower thestability threshold for nonlinearly triggered modes [36]. In this section, the overall dynamics of the simulated cases are summarised usingthe averaged toroidal magnetic energy spectrum during the nonlinear phase, and theevolution of the total plasma current and thermal energy in time. In the early nonlinearphase, the external kink should lead to similar current dynamics as in a constrainedrelaxation process [37, 38]. The current profile flattens, redistributing into the vacuumregion. As a result, the inductance of the plasma decreases, leading to a current spike.5 (a)(b)
Figure 8.
Normalised plasma current (a) during the evolution of the external kink.The relative and absolute amplitude of the current spike (b) are plotted as a metricfor the violence of the initial kink.
The magnitude of the spike is a reasonable metric for the violence of the initial instability.Figure 8 shows the evolution of the normalised current profiles as a function of time. Asexpected, it can be seen that the relative spike amplitude increases with the total plasmacurrent. Eventually for the I p = 4 . M A case, the initial kink instability becomes somild that a current spike is not observed, and the current resistively decays.Figure 9 shows that most of the simulated cases have a very broad magnetic energyspectrum during the nonlinear phase. This is to be expected given that the edge q beingmodeled is below 2, and is therefore highly unstable. In the tokamak case especially,all toroidal harmonics are linearly unstable. As I p decreases, more modes are stabilised,such that the n = 1 component becomes the driving mode in the nonlinear saturation.While the energies of the n = 1, 2 and 3 modes reduce with the plasma current,the higher toroidal harmonics remain relatively constant during the nonlinear phase,indicating that these modes are not just driven by the lower toroidal harmonics, and6 Figure 9.
Magnetic energy spectrum during the nonlinear phase in which the externalkink has saturated. the stability threshold of internal modes has been crossed. This is in particular truefor the n = 4 and 5 modes, as discussed in Section 5.4. The I p = 4 . M A case showsa significantly lower saturated energy than the higher current cases. This seems tobe because the initial kink instability is not violent enough to trigger significant MHDactivity at internal surfaces.The change in thermal energy stored in the plasma volume over time is plotted forall cases in Figure 10. In most cases the thermal energy does not change significantlyin the simulated time frame. This is because the dynamics are strongest in the outerregion of the plasma, while most of the thermal energy is stored in the core, because ofthe assumed parabolic pressure profile. Somewhat surprisingly, the I p = 6 . M A caseloses the most thermal energy. The increase in thermal losses are too large to be due tothe longer timescale of the simulated instability alone. The behaviour at I p = 6 . M A is expected to be because of enhanced transport inside the plasma due to overlappinginternal modes as discussed in Section 5.4. At I p = 4 . M A , the thermal energy doesnot decay significantly. It is shown in Section 5.5 that confinement is maintained overmost of the plasma region in this case.
The general dynamics of the external kink are shown in Figure 11 using time traces ofthe energy, and pseudocolour plots of the temperature, in normalised units. Contoursof density are also overlayed on the temperature plots. The low toroidal harmonicslead the instability, resulting in an initial saturation of the external kink, while highertoroidal harmonics continue to rise after the saturation of the dominant modes. Asthese higher harmonics saturate, the outer region of the plasma becomes increasinglychaotic, leading to a loss of thermal energy.In the tokamak approximation, the n = 2 mode leads the instability, before being7 Figure 10.
Thermal energy during the evolution of the external kink. The I p =6 . M A case has a significantly increased transport of thermal energy out of theplasma. suppressed by the n = 1 mode. The saturation of the n = 1 mode is followed by arelatively stationary phase where the higher toroidal harmonics continue to grow, andsaturate. The combination of multiple toroidal harmonics that are all linearly unstableleads to significant loss of confinement, and heat transport through perpendicularconvection, and parallel conduction in the outer region of the plasma.As I p decreases, higher toroidal harmonics become increasingly stable, such that the n = 1 mode initially drives the other harmonics. This can be seen to lead to similar, butmilder dynamics in the virtual current case with I p = 8 . M A . The I p = 6 . M A and5 . M A cases follow the same initial trend, such that the n = 1 kink is milder. Howeverthe n = 4 and n = 5 modes play a larger role after the initial saturation, especially inthe I p = 6 . M A case. Here, it can be seen that there is a transient peaking of the n = 4 and n = 5 modes around t = 2 . ms . This indicates that the kink triggers furtherinternal instabilities led by these modes. A combination of this effect and the longertimescale of the instability lead to larger thermal losses by the end of the simulation,as can be seen from the constricted temperature plot in the I p = 6 . M A case at t = 3 . ms . Despite the longer time scale, the thermal losses in the I p = 5 . M A case are smaller. This indicates that the combined internal and external instabilitiesare stronger at I p = 6 . M A .The dynamics at I p = 4 . M A are much milder. The initial (2, 1) kink saturateswith a relatively small deformation of the plasma. The internal perturbations led bythe n = 5 mode close to the plasma edge do not penetrate significantly into the plasmaregion, or lead to a secondary deformations of the plasma edge. The time scale ofthe I p = 4 . M A case is orders of magnitude longer than at I p = 6 . M A , withmuch milder thermal losses that are unlikely to lead to a disruption without additionalcontributions to the destabilisation of the plasma that have not been modeled.8 (a) Tokamak Approximation (b) t = 0 .
471 ms (c) t = 0 .
520 ms (d) t = 0 .
633 ms(e) I p = 8 . MA (f) t = 1 .
236 ms (g) t = 1 .
295 ms (h) t = 1 .
530 ms(i) I p = 6 . MA (j) t = 2 .
801 ms (k) t = 3 .
001 ms (l) t = 3 .
178 ms(m) I p = 5 . MA (n) t = 7 .
309 ms (o) t = 7 .
899 ms (p) t = 8 .
238 ms(q) I p = 4 . MA (r) t = 66 .
495 ms (s) t = 72 .
980 ms (s) t = 80 .
029 ms
Figure 11.
Energy evolution of the toroidal harmonics, alongside pseudocolourplots of normalised temperature, during the dynamics. To show the approximatedeformation of the plasma, contour lines of the number density at n = 0 .
1, 0 .
3, 0 . .
7, and 0 . (a) ˆ t = 21 . t = 21 . t = 23 . t = 23 . Figure 12. q profile (a and c) and temperature profiles (b and d) averaged on then=0 normalised poloidal flux, ˆΨ, surfaces at ˆ t = 21 . t = 23 . q = 1 . q = 1 .
75 surfaces.
The internal modes that are triggered can be observed in Figure 12 as island structuresin the surface averaged radial plots of temperature. Magnetic islands lead to largeparallel heat transport, which leads to a local flattening of the temperature profilearound rational surfaces [39]. Two time points are shown during and after the strongMHD activity observed for the I p = 6 . M A case. The surfaces that are averagedover are defined by the n = 0 flux surfaces. The q profile in particular is therefore onlyapproximate, but gives an indication of where the rational surfaces are located. At thetime slices shown, island structures can be inferred from the temperature plots aroundpoints where the profiles are partially flattened. The q profile suggests that the (9, 5)and (7, 4) rational surfaces are important in the dynamics.To explain why the transport of thermal energy out of the plasma is largest at I p = 6 . M A , it is observed that the q profile is flattest in this case both during andafter the dynamics around the q = 1 . .
75 rational surfaces. The external kinkredistributes the plasma current, moving the peak in local current density at the plasmaedge both into the plasma and out into the vacuum. Overall, this leads to a flatteningof the q profile. In the tokamak case, redistribution occurs up to the region of strongnegative shear in the plasma core, where the q profile remains unchanged. The radial0extent of this process reduces with I p , as the kink dynamics become milder.For the I p = 6 . M A case, the redistributed current acts to mostly flatten the outerregion of negative shear (0 . < ˆΨ < .
8) in the q profile. This allows the n = 4 and5 modes to dominate in the mid-region of the plasma, such that their island structuresoverlap. As shown in the temperature plots from a later point in time, the temperatureprofile relaxes significantly after this point, while in the other cases, the islands are toofar apart, and so less thermal energy is transported out of the plasma.The significance of the (9, 5) and (7, 4) modes can also be shown using pseudocolourplots of the poloidal flux, as in Figure 13. The contours show only the n = 4 and n = 5components of the flux for the tokamak, I p = 6 . M A , and I p = 5 . M A cases.All three cases show (7, 4) and (9, 5) structures, indicating that internal modes havebeen triggered. In the tokamak case, (8, 4) and (10, 5) structures can also be observedsuggesting that the (2, 1) external kink still interferes with the internal modes.The modes that are triggered are likely to be either resistive double tearing modes,or infernal modes. There is typically a transition between the two modes depending onhow close the plasma is to the ideal stability threshold [40]. The modes are normallydistinguished by differing dependences of the linear growth rate on resistivity, whichcannot be analysed for the nonlinearly triggered modes in this paper. It can be seen inFigure 11 that at I p = 6 . M A , the n = 4 and n = 5 energies fall during the periodwhere significant thermal energy is lost, indicating that the modes might be pressuredriven. This observation and the apparent occurrence of the (7, 4) mode in plasmas thatdo not seem to have a q = 1 .
75 rational surface give some indication that the nature ofthese modes might be infernal.
The connection length of magnetic field lines to the simulation boundary can be usedas a metric for the loss of confinement. Field lines are traced for a maximum of 100toroidal turns around the torus, or until they connect with the simulation boundary. Theconnection length is determined for 1000 poloidal points on each of 51 radial surfaces,defined by the n = 0 component of the poloidal flux. The normalised harmonic mean iscalculated for each surface, and plotted over time in Figure 14.The dynamics is split into two phases. The initial ergodisation of the kink ischaracterised by a sudden loss in confinement near the plasma edge. As I p increases,this initial stochastisation driven by the external kink penetrates further into the coreregion. Later in time, there is a gradual loss of confinement further inside the plasma asMHD activity grows on the internal rational surfaces. After the strong internal modesin the I p = 6 . M A case, it can be seen that there is a partial recovery of the internalsurfaces. It is expected that if the other simulations were continued further, a similarrecovery may be observed once a comparable amount of thermal energy is lost, andthe drive for the internal modes diminishes. At I p = 4 . M A the initial kink onlypenetrates up to ˆ ψ ≈ .
9, and does not lead to a complete loss of confinement in this1 (a) (b) (c)(d) (e) (f)
Figure 13.
Pseudocolour plots of the n = 4 (a-c) and n = 5 (d-f) componentsof the m (cid:54) = 0 poloidal flux perturbation for the tokamak approximation (a and d), I p = 6 . M A (b and e) and I p = 5 . M A (c and f) cases, at ˆ t = 22 .
48. (7, 4) and(9, 5) internal modes have been triggered in all cases. region. The subsequent internal modes break up the flux surfaces close to the plasmaedge, but confinement is maintained over most of the plasma. Later in time, the edgecurrent denisity gradient and total plasma current reduce due to resistive effects. Thiscauses the edge safety factor to rise, and the drive for the instability to relax. As such,the energy in the toroidal perturbations decays, and the associated magnetic islandsshrink, there is a partial of recovery of confinement in the simulation domain, such thatlosses are mostly in the vacuum region.Poincar´e plots are shown in Figure 15 at the end of the simulated time. A dominant(2, 1) island structure is observed in all cases. In the tokamak approximation, the islandstructure is deformed by the remnants of the sub-dominant (4, 2) mode. Relatively littlestructure can be observed in the region where confinement is lost. This is because fieldlines very quickly connect to the simulation boundary, via the large (2, 1) magneticislands at the edge. This can be seen by plotting the strike points of the field lineson the simulation boundary as shown in Figure 16. The strike points trace out theintersection of the external islands with the simulation boundary, showing that themagnetic field has a very clear structure in the vacuum region.In the tokamak approximation, a small (2, 1) island has formed in the core, showing2 (a) Tokamak Approximation (b) I p = 6 . MA (c) I p = 4 . MA Figure 14.
Plots of normalised connection length to the simulation domain boundaryat the end of the simulated timeframe. (a) Tokamak (b) I p = 6 . MA (c) I p = 4 . MA Figure 15.
Poincar´e plots taken at ˆ t = 23 . I p = 6 . M A (b) and I p = 4 . M A case. the radial extent of the (2, 1) instability is larger than in cases with lower plasma current.This can also be seen from the visible kink in the core flux surfaces compared to thevirtual current cases. The (2, 1) deformation of the plasma reduces with I p . Similar tothe connection length plots, the Poincar´e plot of the I p = 4 . M A case show relativelylittle deformation of the plasma. Large islands are still apparent, but they are much lessergodic, which leads to better confinement of the field lines within the external (2, 1)island structure. Smaller (7, 4) and (9, 5) island chains have been identified just insidethe plasma region, but they remain at low enough energies to only mildly degrade theconfinement in the outer region of the plasma.
6. Conclusion
External kinks have been modeled with JOREK nonlinearly in axisymmetricapproximations of a QA configuration with varying fractions of external rotational3
Figure 16.
Strike points of field lines on the simulation boundary for the I p =6 . M A case. θ and φ are the poloidal and toroidal angle respectively. transform. This approach has been validated in the linear phase against CASTOR3D,showing reasonable agreement in the growth rates and eigenfunctions when comparedwith a 3D code, which supports the use of this model as a preliminary tool inunderstanding the dynamics of such configurations. In the nonlinear phase, the violenceof the initial instability is reduced with increased contribution of virtual currents tothe rotational transform, as measured by the current spike observed in the plasma.Nonlinearly triggered modes, in particular for this case the (7, 4) and (9, 5), mean thatthe loss in confinement can remain large, even after an appreciable external rotationaltransform has been introduced. This means that the variation of thermal losses withincreasing external rotational transform can be non-monotonic. There are cases at largeexternal rotational transform where internal modes are mild, indicating that the initialexternal modes could lead to a nonlinearly stable state.Further extensions to this work should consider the dependence of the nonlineardynamics on the edge safety factor and plasma pressure, as well as further destabilisingeffects such as vertical modes and impurities introduced from the plasma touching wallstructures. Eventually a 3D code will be necessary to look at the detailed dynamics,especially near the stability threshold for these kinks. The results of this paper havebeen carried out to inform the direction of future studies with a stellarator capable 3Dnonlinear code. When such a tool becomes available, the limitations of this work canbe overcome. Appendix
As part of the validation of the models used in this paper, resistive instabilities observedexperimentally in W7-X have been modeled and compared using CASTOR3D, JOREKand TM1 [28]. These instabilities are induced by Electron Cyclotron Current Drive(ECCD) near the core of the device, which leads to the formation of two q = 1 rational4surfaces. Resistive kink and double tearing modes can be linearly unstable at thesesurfaces [41].Axisymmetric approximations of a low beta W7-X equilibrium have beenconstructed, using both a circular boundary, and one defined by the n = 0 Fourierharmonics of the last closed flux surface of the equilibrium. The linear growth ratesof resistive modes are compared against full 3D results from CASTOR3D, as shownin Figure A1. In this instability, the tokamak approximation is more stable than thevirtual current model, despite having a larger plasma current. The reason for thisis because the current that represents the external rotational transform of a net zerocurrent stellarator is monotonically increasing towards the plasma edge. When thesecurrents are taken up by the plasma, as in the tokamak approximation, they reduce thecurrent density gradients around the rational surfaces, thereby reducing the drive for theinstabilities. The sharp change in the n = 1 growth rate of the tokamak approximationin Figure A1 (c) is because the (1, 1) mode can be either a resistive kink or doubletearing mode. The results shown for CASTOR3D are for double tearing modes. In thetokamak approximation, a resistive kink is observed in JOREK with a higher growthrate.The comparison shows good agreement between the three codes. The axisymmetricapproximations and the full 3D calculations match well even at relatively high toroidalmode numbers. While geometric effects play only a minor role in the linear dynamics,using the n = 0 Fourier components leads to a better match with the full 3D calculationsthan cylindrical approximations. Further studies were carried out for these resistivemodes, exploring the full 3D linear stability in more detail, and the nonlinear dynamicsin cylindrical geometry [28, 41]. Results have shown that axisymmetric approximationswork well, even for advanced stellarators, in the low beta, high aspect ratio limit. Acknowledgements
The authors would like to thank Qingquan Yu, for providing linear growth ratesused in the validation of the models in this paper, and Florian Hindenlang, NikitaNikulsin, Sophia Henneberg, and Carolin N¨uhrenberg for their contributions to thiswork through valuable insights and fruitful discussions. This work was supported bythe Max-Planck/Princeton Center for Plasma Physics. The authors acknowledge accessto the EUROfusion High Performance Computer (Marconi-Fusion) and the JFRS-1supercomputer in Japan to perform the presented simulations.
References [1] Henneberg S, Drevlak M, N¨uhrenberg C, Beidler C, Turkin Y, Loizu J and Helander P 2019
Nuclear Fusion Nuclear Fusion Physics of Plasmas (a)(b) Circular(c) Axisymmetric and full 3D Figure A1.
Equilibrium q profiles (a) of W7-X equilibria with and without ECCD,and linear growth rates of W7-X resistive modes calculated with JOREK, CASTOR3D,and TM1, using a circular approximation (b) and the n = 0 fourier coefficients ofthe W7-X boundary (c). The full 3D growth rates are compared with the n = 0approximations of W7-X.[4] ArchMiller M, Cianciosa M, Ennis D, Hanson J, Hartwell G, Hebert J, Herfindal J, Knowlton S,Ma X, Maurer D et al. Physics of Plasmas et al. Fusion science and technology Physics of Plasmas Nuclear Fusion S128–S202[8] Strumberger E and G¨unter S 2016
Nuclear Fusion Nuclear Fusion et al. Nuclear Fusion [11] Todo Y, Seki R, Spong D, Wang H, Suzuki Y, Yamamoto S, Nakajima N and Osakabe M 2017 Physics of plasmas [13] Hoelzl M, Huijsmans G, Pamela S, Becoulet M, Nardon E, Artola F, Nkonga B, Atanasiu C,Bandaru V, Bhole A et al. submitted to Nuclear Fusion. preprint at arxiv:2011.09120 [14] Nikulsin N, Hoelzl M, Zocco A, Lackner K and G¨unter S 2019 Physics of Plasmas submitted toJournal of Plasma Physics [16] Greene J M and Johnson J L 1961 The Physics of Fluids Ideal magnetohydrodynamics (New York: Plenum Press)[18] Greene J M, Johnson J L, Kruskal M D and Wilets L 1962
The Physics of Fluids Nuclear Fusion Nuclear Fusion NuclearFusion [22] Fredrickson E, Zarnstorff M and Lazarus E 2007 Fusion science and technology Journal of computational physics
SIAM Journal on scientific and statistical computing Parallel Computing pp 897–906[28] Yu Q, Strumberger E, Igochine V, Lackner K, Laqua H, Zanini M, Braune H, Hirsch M, H¨ofel U,Marsen S et al. Nuclear Fusion [29] Hirshman S P and Whitson J 1983
The Physics of fluids ESAIM: Mathematical Modelling andNumerical Analysis arXiv preprint arXiv:1508.04911 [32] Hoelzl M, Merkel P, Huysmans G, Nardon E, Strumberger E, McAdams R, Chapman I, G¨unterS and Lackner K 2012 Coupling JOREK and STARWALL codes for non-linear resistive-wallsimulations Journal of Physics: Conference Series [33] Artola F J, Beyer P, Huijsmans G, Loarte A and Hoelzl M 2018
Free-boundary simulationsof MHD plasma instabilities in tokamaks
Ph.D. thesis Aix-Marseille Universit´e URL https://hal-amu.archives-ouvertes.fr/tel-02012234v1 [34] Boozer A H 1981
The Physics of Fluids Physics of Fluids B: Plasma Physics Nuclear fusion Reviews of Modern Physics The Physics of Fluids Plasma physics and controlled fusion Physics of Fluids B: Plasma Physics et al. Nuclear Fusion60