Code JOREK Research Articles

MHD simulations of small ELMs at low triangularity in ASDEX Upgrade

The development of small and no-ELM regimes for ITER is a high priority topic due to the risks associated with type-I ELMs. By considering non-linear extended magnetohydrodynamic (MHD) simulations of the ASDEX Upgrade tokamak with the JOREK code, we probe a regime that avoids type-I ELMs completely, provided that the separatrix density is high enough. The dynamics of the pedestal in this regime are observed to be qualitatively similar to the so-called quasi-continuous exhaust regime in several ways. Repetitive type-I ELMs are substituted by roughly constant levels of outward transport, caused by peeling-ballooning modes (with dominant ballooning characteristics) which are localised in the last 5% of the confined region (in normalised poloidal flux). The simulated low triangularity plasma transitions to a type-I ELMy H-mode if the separatrix density is sufficiently reduced or if the input heating power is sufficiently increased. The stabilising factors that play a role in the suppression of the small ELMs are also investigated by analysing the simulations, and the importance of including diamagnetic effects in the simulations is highlighted. By considering a scan in the pedestal resistivity and by comparing the poloidal velocity of the modes to theoretical estimates for ideal and resistive modes, we identify the underlying instabilities as resistive peeling-ballooning modes. Decreasing the resistivity below experimentally-relevant conditions (i.e. going towards ideal MHD), the peeling-ballooning modes that constrain the pedestal below the type-I ELM stability boundary display sharply decreasing growth rates.

A generalised formulation of G-continuous Bezier elements applied to non-linear MHD simulations

The international tokamak ITER is progressing towards assembly completion and first-plasma operation, which will be a physics and engineering challenge for the fusion community. In the preparation for ITER experimental scenarios, non-linear MHD simulations are playing an essential role to actively understand and predict the behaviour and stability of tokamak plasmas in future fusion power plant. The development of MHD codes like JOREK is a key aspect of this research effort, and provides invaluable insight into the plasma stability and the control of global and localised plasma events, like Edge-Localised-Mode and disruptions. In this paper, we present an operational implementation of a new, generalised formulation of Bezier finite-elements applied to the JOREK code, a significant advancement from the previously G1-continuous bi-cubic Bezier elements. This new mathematical method enables any polynomial order of Bezier elements, with a guarantee of G-continuity at the level of (n−1)/2, for any odd n, where n is the order of the Bezier polynomials. The generalised method is defined, and a rigorous mathematical proof is provided for the G-continuity requirement. Key details on the code implementation are mentioned, together with a suite of tests to demonstrate the mathematical reliability of the finite-element method, as well as the practical usability for typical non-linear tokamak MHD simulations. A demonstration for a state-of-the-art simulation of an Edge-Localised-Mode instability in the future ITER tokamak, with realistic grid geometry, finalises the study.

On the role of preexisting MHD activity for the plasma response to massive deuterium injection

As part of a reliable disruption mitigation system (SPI) for ITER, pure deuterium shattered pellet injection (SPI) has been proposed as a way of avoiding hot tail runaway electron generation. It offers the possibility of diluting the plasma and, thereby, cooling it down by a large factor without immediately triggering a thermal quench (TQ). However, the reliability of this and similar SPI approaches could be reduced by preexisting MHD modes, which are usually present during the pre-TQ phase, when the disruption mitigation scheme is being triggered. To address this question, this theoretical study investigates massive deuterium injection into an MHD active ASDEX Upgrade plasma using the non-linear MHD code JOREK. Cases with and without preexisting 2/1 islands are studied. Scans are performed in the preexisting island size, the number of atoms injected, and the relative phase of the injection location with respect to the island. Realistic values of resistivity and heat diffusion anisotropy are considered. This provides insights into the physical mechanisms at play and the relevant time scales involved. Results largely indicate that plasma dilution by deuterium also seems to work reliably in the presence of preexisting MHD activity. Nevertheless, when injecting in phase with the X-point of a large preexisting island, the TQ can occur earlier than without. Altogether, simulations increase confidence in the reliability of plasma dilution by deuterium injection and its applicability to ITER.

A simulation chain for reflectometry and non-linear MHD: type-I ELM case

A complete chain from a non-linear MHD plasma model simulation through full-wave code simulations implementing synthetic conventional reflectometry is established. For this purpose, the two-dimensional full-wave code REFMUL is employed together with MHD descriptions obtained from the JOREK code. First results of the integrated modeling are presented here where a type-I ELM crash, leading to a fast collapse of the H-mode pedestal, was taken as case-study. The REFMUL simulations were customized to implement synthetic reflectometry in conventional set-up using fixed frequency probing with O-mode waves. Posing challenging conditions for reflectometry, the type-I ELM crash reveals some of the merits and caveats of the diagnostic technique. This work also opens up the possibility to extend modeling to other MHD or ELM studies and provide support to experimental observations with reflectometry.

Thermal quench and current profile relaxation dynamics in massive-material-injection-triggered tokamak disruptions

3D non-linear magnetohydrodynamic simulations of a disruption triggered by a massive injection of argon gas in JET are performed with the JOREK code. The key role of the thermal drive of the m = 2, n = 1 tearing mode (i.e. the drive from helical cooling inside the island) in the disruption process is highlighted by varying the amplitude and position of the argon source across simulations, and also during a simulation. In cases where this drive persists in spite of the development of magnetic stochasticity, which is favoured by moving the argon source in an ad hoc way from the plasma edge into the 2/1 island at some point in the simulation, a relaxation in the region (roughly) takes place. This relaxation generates a plasma current spike comparable to the experimental one. Simulations are compared in detail to measurements via synthetic diagnostics, validating the model to some degree.

Comparing spontaneous and pellet-triggered ELMs via non-linear extended MHD simulations

Injecting frozen deuterium pellets into an ELMy H-mode plasma is a well established scheme for triggering edge localized modes (ELMs) before they naturally occur. This paper presents non-linear simulations of spontaneous type-I ELMs and pellet-triggered ELMs in ASDEX Upgrade performed with the extended MHD code JOREK. A thorough comparison of the non-linear dynamics of these events is provided. In particular, pellet-triggered ELMs are simulated by injecting deuterium pellets into different time points during the pedestal build-up described in A Cathey et al (2020 Nuclear Fusion 60 124007). Realistic ExB and diamagnetic background plasma flows as well as the time dependent bootstrap current evolution are included during the build-up to accurately capture the balance between stabilising and destabilising terms for the edge instabilities. Dependencies on the pellet size and injection times are studied. The spatio-temporal structures of the modes and the resulting divertor heat fluxes are compared in detail between spontaneous and triggered ELMs. We observe that the premature excitation of ELMs by means of pellet injection is caused by a helical perturbation described by a toroidal mode number of n = 1. In accordance with experimental observations, the pellet-triggered ELMs show reduced thermal energy losses and a narrower divertor wetted area with respect to spontaneous ELMs. The peak divertor energy fluence is seen to decrease when ELMs are triggered by pellets injected earlier during the pedestal build-up.

Experimental study of ELM induced fast-ion transport using passive FIDA spectroscopy at the ASDEX Upgrade tokamak

Measurements using a recently installed edge fast-ion D-alpha (FIDA) diagnostic at the ASDEX Upgrade tokamak show a clear effect of edge localised modes (ELMs) on the passive FIDA signals. While a reduction in the passive FIDA emission is observed in the scrape-off layer (SOL) region, measurements close to the last closed flux surface show an increase in signals shortly after ELMs, followed by a decrease. The decrease provides a clear sign of fast-ion losses in the SOL, while the increase can be explained by an enhanced neutral density during ELMs inside the plasma. In addition, small ELMs are observed, which barely change the neutral density and plasma position but still cause significant changes in the passive FIDA signals. A comparison of the measurements with forward modelling shows that 60% to 80% of the fast ions are lost by ELMs outside the last closed flux surface. In addition, a 20% decrease of the fast-ion density in a range up to 4 cm within the last closed flux surface can be inferred. This range agrees well with the latest modelling results of ELMs using the non-linear MHD code JOREK and shows that less than 0.3% of all fast ions are lost by ELMs.

Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection

The radiation response and the MHD destabilization during the thermal quench after a mixed species shattered pellet injection with impurity species neon and argon are investigated via 3D non-linear MHD simulation using the JOREK code. Both the n = 0 global current profile contraction and the local helical cooling at each rational surface caused by the pellet fragments are found to be responsible for MHD destabilization after the injection. Significant current driven mode growth is observed as the fragments cross low order rational surfaces, resulting in rapidly inward propagating stochastic magnetic field, ultimately causing the core temperature collapse. The thermal quench (TQ) is triggered as the fragments arrive on the q = 1 or q = 2 surface depending on the exact q profile and thus mode structure. When injecting from a single toroidal location, strong radiation asymmetry is found before and during the TQ as a result of the unrelaxed impurity density profile along the field line and asymmetric outward heat flux. Such asymmetry gradually relaxes over the course of the TQ, and is entirely eliminated by the end of it. Simulation results indicate that the aforementioned asymmetric radiation behavior could be significantly mitigated by injection from toroidally opposite locations, provided that the time delay between the two injectors is shorter than 1 ms. It is also found that the MHD response are sensitive to the relative timing and injection configuration in these multiple injection cases.

Parametric decay instabilities near the second-harmonic upper hybrid resonance in fusion plasmas

Parametric decay instabilities (PDIs) lead to the generation of strong frequency-shifted radiation when powerful X-mode polarized microwaves, injected for electron cyclotron resonance heating (ECRH) of fusion plasmas, cross a region of non-monotonic electron density near the second-harmonic upper hybrid resonance (UHR). For the standard second-harmonic X-mode ECRH scenarios used at the ASDEX Upgrade tokamak, the second-harmonic UHR occurs near the plasma edge, meaning that PDIs occur in connection with phenomena leading to non-monotonic edge electron density profiles, e.g. blobs, edge-localized modes (ELMs), and inter-ELM modes. We present the first study of strong signals near half the ECRH frequency in fusion-relevant plasmas and demonstrate their PDI-like features, such as an ECRH power threshold required for their occurrence, through analog modulations of the ECRH power. The signals near half the ECRH frequency are found to be far more prevalent than expected based on previous theories, and we hence present a qualitative explanation of their origin, in terms of cross-polarized microwaves generated through non-linear interactions of the waves involved in the PDIs, which captures the basic features of the experimental observations. We additionally present examples of the PDI-related microwave signals just below the ECRH frequency which indicate the location of the microwave source. Furthermore, based on a non-linear magnetohydrodynamics simulation of an ELM crash in ASDEX Upgrade, performed using the JOREK code, the duration of the PDI-like microwave spikes observed in connection with ELMs is shown to match the transit time of an ELM filament though an ECRH beam. Finally, the PDI power threshold expected theoretically, calculated using the electron density and temperature profiles from the JOREK simulation, is shown to match the experimentally observed values.

Scattering of ion cyclotron range of frequency waves by filaments and ELMs

The scattering of ion cyclotron range of frequency (ICRF) waves by filaments and ELMs has been comprehensively studied via modeling. The influence of filaments on the propagation of ICRF waves is studied by building a 2D COMSOL model which solves the wave equation. It is benchmarked against the Mie scattering theoretical model. Cases either with perfect matching layer boundary condition or with poloidal periodic boundary condition are considered. Parameter scans, including the filament density, the filament radius, the number of filaments as well as the distance between filaments, are performed. Then, the influence of ELMs on the propagation of ICRF waves in realistic geometry is studied by importing density distributions calculated by 3D non-linear MHD code JOREK into the 3D antenna code RAPLICASOL. The evolution of density and wave fields during an ELM is analyzed for cases with a low and high scrape-off layer density, respectively. Both the COMSOL simulation and the realistic simulation with JOREK and RAPLICASOL show that the wave fields and Poynting flux can be significantly and globally modified by the ELM filaments. Due to the mutual influence of many density filaments, radially elongated and poloidally distributed stripe structures of the wave fields develop, leading to poloidally inhomogeneous wave fields. Furthermore, the influence of ELMs on the antenna power and antenna coupling resistance is investigated with experimental results in ASDEX Upgrade.

Global ITG eigenmodes: From ballooning angle and radial shift to Reynolds stress and nonlinear saturation

We present global linear and nonlinear simulations of ion temperature gradient instabilities based on a fluid formulation, with an adapted version of the JOREK code. These simulations are performed in realistic global tokamak equilibria based on the solution of the Grad–Shafranov equation. Benchmarking of linear growth rates was successfully completed with respect to previously published data. We find two distinct types of eigenstructures, depending on the magnetic shear. For high shear, when the coupling of poloidal harmonics is strong, ballooning-type eigenmodes are formed, which are up-down asymmetric with a finite ballooning angle, θ0. The poloidal harmonics which form the global eigenmode are found to demonstrate a radial shift, being centered well outside of their corresponding rational surface. Stronger diamagnetic effects increase both θ0 and proportionately shift the m harmonics to larger radii (by as much as two rational surfaces). In the low shear regime, the unstable eigenmodes become narrowly localized between neighboring pairs of rational surfaces, and exhibit no up-down asymmetry. Our simulations also show the generation of finite Reynolds stress due to nonlocal/global profile effects. This stress possesses both poloidally symmetric (n=m=0) and asymmetric (finite-m) components. Turbulent saturation in nonlinear simulations is demonstrated for both shear regimes.

Nonlinear modeling of the effect of n = 2 resonant magnetic field perturbation on peeling-ballooning modes in KSTAR

Using the nonlinear 3D MHD code JOREK with reduced MHD equations (visco-resistive MHD), we have successfully simulated a recent n = 2 resonant magnetic perturbation (RMP)-driven edge localized mode (ELM) suppression in KSTAR. We have found that such ELM suppression has been attributable not only to the degraded pedestal but also to the direct coupling between the peeling-ballooning mode (PBM) and RMP-driven plasma response. Notably, the pedestal pressure gradient is reduced as the radial transport is enhanced because of the formation of the stochastic layer and increased convection fluxes due to tearing and the kink-peeling mode driven by RMPs. The increased transport in the stochastic layer is due to the parallel transport across the stochastic fields, described by the Braginskii model in the simulation. While the linear stability of the PBM improves owing to the degraded pedestal, it is not a sole contributor to ELM suppression, in that the nonlinear mode coupling plays a more critical role. This outcome is consistent with previous studies where mode coupling affects the ELM mitigation or suppression. In addition, PBM locking has been numerically achieved during the ELM suppression phase, which may support the relationship between at the pedestal and the onset of ELM suppression. We suggest that PBM locking can enhance the mode interactions between RMPs and PBMs, which is significant for ELM suppression.

Modeling of TAE mode excitation with an antenna in realistic X-point geometry

Experimentally, it is observed that Toroidal Alfvén Eigenmodes (TAEs) are difficult to excite with an external antenna when the plasma is in X-point geometry. Here, the effect of the X-point geometry on the efficiency of the TAE excitation with the external antenna is investigated. In the first part of this paper, the influence of the near-last closed flux surface layer from the core side on the damping of the TAE modes is calculated using the linear resistive eigenvalue MHD code CASTOR. The resistive damping is identified as the main cause of the TAE damping in the open gap in the Alfvén continuum. It is shown that one aspect of the difficulty of excitation of the TAE modes in X-point geometry is an increased damping from the region inside the separatrix. However, the increased damping with the plasma boundary approaching the separatrix is not general and depends on the density profile shape. The second part of this paper discusses the influence of the TAE behavior in the limiter and X-point geometries including the scrape-off layer (SOL) in the reduced MHD code JOREK. It is shown that the dominant effect on the damping of the original TAE mode observed in the limiter configuration is caused by the additional damping in the region of open field lines, i.e., the SOL.

Non-linear magnetohydrodynamic simulations of pellet triggered edge-localized modes in JET

Non-linear magnetohydrodynamic simulations of pellet-triggered edge-localized modes (ELMs) in JET plasma have been carried out with the JOREK code. The pellet particle fuelling efficiency and the power flux at the divertor target during the pellet-triggered ELM have been studied. The pellet injection in unstable plasma delivers the particle fuelling but the pellet fuelling rate is smaller than the rate of particle loss during the pellet triggered ELM. The JOREK simulations estimate the power flux at the divertor target and found good agreement with the experimental observation. The energy deposition of the pellet triggered ELM shows a toroidally asymmetric profile. However, and due to this toroidal asymmetry, this effect cannot be captured by the existing layout of the divertor infra-red (IR) cameras available in JET. This work highlights the benefit of having a larger number of IR cameras to analyse the heat flux for the experiments which are assumed to be toroidally asymmetric, such as the pellet and/or gas injection experiments.

Simulating the nonlinear interaction of relativistic electrons and tokamak plasma instabilities: Implementation and validation of a fluid model

For the simulation of disruptions in tokamak fusion plasmas, a fluid model describing the evolution of relativistic runaway electrons and their interaction with the background plasma is presented. The overall aim of the model is to self-consistently describe the nonlinear coupled evolution of runaway electrons (REs) and plasma instabilities during disruptions. In this model, the runaway electrons are considered as a separate fluid species in which the initial seed is generated through the Dreicer source, which eventually grows by the avalanche mechanism (further relevant source mechanisms can easily be added). Advection of the runaway electrons is considered primarily along field lines, but also taking into account the E×B drift. The model is implemented in the nonlinear magnetohydrodynamic (MHD) code jorek based on Bezier finite elements, with current coupling to the thermal plasma. Benchmarking of the code with the one-dimensional runaway electron code go is done using an artificial thermal quench on a circular plasma. As a first demonstration, the code is applied to the problem of an axisymmetric cold vertical displacement event in an ITER plasma, revealing significantly different dynamics between cases computed with and without runaway electrons. Though it is not yet feasible to achieve fully realistic runaway electron velocities close to the speed of light in complete simulations of slowly evolving plasma instabilities, the code is demonstrated to be suitable to study various kinds of MHD-RE interactions in MHD-active and disruption relevant plasmas.

A wall-aligned grid generator for non-linear simulations of MHD instabilities in tokamak plasmas

Block-structured mesh generation techniques have been well addressed in the CFD community for automobile and aerospace studies, and their applicability to magnetic fusion is highly relevant, due to the complexity of the plasma-facing wall structures inside a tokamak device. Typically applied to non-linear simulations of MHD instabilities relevant to magnetically confined fusion, the JOREK code was originally developed with a 2D grid composed of isoparametric bi-cubic Bézier finite elements, that are aligned to the magnetic equilibrium of tokamak plasmas (the third dimension being represented by Fourier harmonics). To improve the applicability of these simulations, the grid-generator has been generalised to provide a robust extension method, using a block-structured mesh approach, which allows the simulations of arbitrary domains of tokamak vacuum vessels. Such boundary-aligned grids require the adaptation of boundary conditions along the edge of the new domain. Demonstrative non-linear simulations of plasma edge instabilities are presented to validate the robustness of the new grid, and future potential physics applications for tokamak plasmas are discussed. The methods presented here may be of interest to the wider community, beyond tokamak physics, wherever imposing arbitrary boundaries to quadrilateral finite elements is required.

Non-linear modeling of the threshold between ELM mitigation and ELM suppression by resonant magnetic perturbations in ASDEX upgrade

The interaction between Edge-Localized Modes (ELMs) and Resonant Magnetic Perturbations (RMPs) is modeled with the magnetohydrodynamic code JOREK using experimental parameters from ASDEX Upgrade discharges. According to the modeling, the ELM mitigation or suppression is optimal when the amplification of both tearing and peeling-kink responses results in a better RMP penetration. The ELM mitigation or suppression is not only due to the reduction of the pressure gradient but predominantly arises from the toroidal coupling between the ELMs and the RMP-induced mode at the plasma edge, forcing the edge modes to saturate at a low level. The bifurcation from ELM mitigation to ELM suppression is observed when the RMP amplitude is increased. ELM mitigation is characterized by rotating modes at the edge, while the mode locking to RMPs is induced by the resonant braking of the electron perpendicular flow in the ELM suppression regime.

Kinetic modeling of ELM-induced tungsten transport in a tokamak plasma

Impurity accumulation in the core plasma leads to fuel dilution and higher radiative losses that can lead to loss of the H-mode, to thermal collapse of the plasma, and eventually even to a disruption in tokamaks. In present experiments, it has been shown that Edge Localized Modes (ELMs) at sufficiently high frequency are required to prevent W accumulation in the core, by expelling impurities from the edge plasma region, thus preventing their penetration into the plasma core. We present a full-orbit particle extension of the MHD code JOREK suitable for simulating impurity transport during ELMs. This model has been applied to the simulation of an ELM crash in ASDEX Upgrade, where we have quantified the displacement of W particles across flux surfaces. The transport mechanism is shown to be the particle E × B-drifts due to the electric field created by the MHD instability underlying the ELM. In- and outward transport is observed as a series of interchange motions, leading to a superdiffusive behavior. This causes not only the particles near the plasma pedestal to move outwards but also the particles outside of the pedestal to move inwards. This has important consequences for operation with W in ITER, where it is expected to be screened in the pedestal, and ELMs are shown here to increase the core W density. A comparison with existing diffusive modeling is made, showing a qualitative agreement and the limitations of this simplified modeling approach.

Numerical study of tearing mode seeding in tokamak X-point plasma

A detailed understanding of island seeding is crucial to avoid neoclassical tearing modes and their negative consequences like confinement degradation and disruptions. In the present work, we investigate the growth of 2/1 islands in response to magnetic perturbations. Although we use externally applied perturbations produced by resonant magnetic perturbation (RMP) coils for this study, the results are directly transferable to island seeding by other MHD instabilities creating a resonant magnetic field component at the rational surface. Experimental results for 2/1 island penetration from ASDEX Upgrade are presented extending previous studies. Simulations are based on an ASDEX Upgrade L-mode discharge with low collisionality and active RMP coils. Our numerical studies are performed with the 3D, two-fluid, nonlinear MHD code JOREK. All three phases of mode seeding observed in the experiment are also seen in the simulations: first, a weak response phase characterized by large perpendicular electron flow velocities followed by a fast growth of the magnetic island size accompanied by a reduction of the perpendicular electron velocity and finally the saturation to a fully formed island state with perpendicular electron velocity close to zero. Thresholds for mode penetration are observed in the plasma rotation as well as in the RMP coil current. A hysteresis of the island size and electron perpendicular velocity is observed between the ramping up and down of the RMP amplitude consistent with an analytically predicted bifurcation. The transition from dominant kink/bending to tearing parity during the penetration is investigated.

3D non-linear MHD simulation of the MHD response and density increase as a result of shattered pellet injection

The MHD response and the penetration of a deuterium shattered pellet into a JET plasma is investigated via the non-linear reduced MHD code JOREK with the neutral gas shielding (NGS) ablation model. The dominant MHD destabilizing mechanism by the injection is identified as the local helical cooling at each rational surface, as opposed to the global current profile contraction. Thus the injected fragments destabilize each rational surface as they pass through them. The injection penetration is found to be much better compared to MGI, with the convective transport caused by core MHD instabilities (e.g. 1/1 kink) contributing significantly to the core penetration. Moreover, the injection with realistic JET SPI system configurations is simulated in order to provide some insights into future operations, and the impact on the total assimilation and penetration depth of varying injection parameters such as the injection velocity or fineness of shattering is assessed. Further, the effect of changing the target equilibrium temperature or q profile on the assimilation and penetration is also investigated. Such analysis will form the basis of further investigation into a desirable configuration for the future SPI system in ITER.