Simulations of COMPASS Vertical Displacement Events with a self-consistent model for halo currents including neutrals and sheath boundary conditions
F.J. Artola, A. Loarte, E. Matveeva, J. Havlicek, T. Markovic, J. Adamek, J. Cavalier, L. Kripner, G.T.A. Huijsmans, M. Lehnen, M. Hoelzl, R. Panek, COMPASS team, JOREK team
SSimulations of COMPASS Vertical Displacement Events with aself-consistent model for halo currents including neutrals and sheathboundary conditions
F.J. Artola , A. Loarte , E. Matveeva , J. Havlicek , T. Markovic , J. Adamek , J.Cavalier , L. Kripner , G.T.A. Huijsmans , M. Lehnen , M. Hoelzl , R. Panek , theCOMPASS team , and the JOREK team ITER Organization, Route de Vinon sur Verdon, 13067 St Paul Lez Durance Cedex, France Max Planck Institute for Plasmaphysics, Boltzmannstr. 2, 85748 Garching, Germany CEA, IRFM, F-13108 St. Paul-lez-Durance cedex, France Eindhoven University of Technology, 5612 AP Eindhoven, The Netherlands Charles University, Faculty of Mathematics and Physics, Prague, Czech Republic Institute of Plasma Physics of CAS, Za Slovankou 3, 182 00 Prague, Czech Republic please refer to [M Hoelzl, G T A Huijsmans, S J P Pamela, M Becoulet, E Nardon, F J Artola, B Nkonga etal, Nuclear Fusion, in preparation] Corresponding author:
A. Loarte
Email: [email protected]
Abstract
The understanding of the halo current properties during disruptions is key to design and operatelarge scale tokamaks in view of the large thermal and electromagnetic loads that they entail. Forthe first time, we present a fully self-consistent model for halo current simulations including neutralparticles and sheath boundary conditions. The model is used to simulate Vertical DisplacementEvents (VDEs) occurring in the COMPASS tokamak. Recent COMPASS experiments have shownthat the parallel halo current density at the plasma-wall interface is limited by the ion saturationcurrent during VDE-induced disruptions. We show that usual MHD boundary conditions can leadto the violation of this physical limit and we implement this current density limitation through aboundary condition for the electrostatic potential. Sheath boundary conditions for the density, theheat flux, the parallel velocity and a realistic parameter choice (e.g. Spitzer η and Spitzer-Härm χ k values) extend present VDE simulations beyond the state of the art. Experimental measurementsof the current density, temperature and heat flux profiles at the COMPASS divertor are comparedwith the results obtained from axisymmetric simulations. Since the ion saturation current density( J sat ) is shown to be essential to determine the halo current profile, parametric scans are performedto study its dependence on different quantities such as the plasma resistivity and the particle andheat diffusion coefficients. In this respect, the plasma resistivity in the halo region broadenssignificantly the J sat profile, increasing the halo width at a similar total halo current. The operation of large scale tokamaks such as ITER and DEMO requires a Disruption MitigationSystem (DMS) to minimize the impact of disruption loads on the plasma facing components andon the vacuum vessel [14]. The design and optimization of such a system must be guided byexperiments and validated numerical codes. In this respect, high fidelity simulations should beable to predict the worst case scenarios in terms of heat and electromagnetic loads in order toquantify the load reduction provided by the different mitigation schemes.During disruptions, currents flowing along the open field lines (halo currents) can producelarge forces on the vacuum vessel and also large heat loads on the plasma facing components byconverting the plasma magnetic energy into thermal energy via Ohmic heating. The parametric a r X i v : . [ phy s i c s . p l a s m - ph ] J a n ependencies of the current density and the heat flux profiles during disruptions are not yet wellunderstood. MHD simulations for disruptions are numerically challenging and therefore previoussimulations have typically oversimplified the boundary conditions for different plasma quantities(e.g. Dirichlet conditions for temperature and density).In this paper we present, for the first time, a fully self-consistent model for halo current sim-ulations including neutral particles and a set of advanced boundary conditions. Sheath boundaryconditions for the density, the heat flux and the parallel velocity together with a realistic param-eter choice (e.g. Spitzer η and Spitzer-Härm χ k values) extend non-linear MHD VDE simulationsbeyond the state of the art. Additionally we introduce a boundary condition for the electrostaticpotential in order to limit the halo current density to the ion saturation current. This limitationcomes from basic sheath plasma physics and has recently been observed in COMPASS disruptionexperiments [1]. Neutral particles are required to calculate the evolution of the boundary density,which is key to determine the ion saturation current density. As a first step, the neutral particledensity is modeled with a fluid continuity equation and an effective diffusion coefficient. For thiswork we use the JOREK code to perform axisymetric MHD simulations of COMPASS VDEs andcompare the modeled results with experimental measurements. Since the ion saturation currentdensity ( J sat ) will be shown to be essential to determine the halo current profile, parametric scansare performed to study its dependence on different quantities such as the plasma resistivity andthe particle and heat diffusion coefficients. Recent simulations with the NIMROD code have ex-plored the effect of different boundary conditions for the plasma density, temperature and velocityon VDEs [4]. However the important effects of neutral particles and the ion saturation currentlimitation were not included.This paper is organized as follows: In section 2 the MHD model and the boundary conditionsare explained in detail. The numerical setup is presented in section 3 together with the parameterchoice used to simulate the COMPASS VDE experiments. The large effect of the sheath boundaryconditions on the evolution of the halo current is discussed in 4. In that section, the influence ofplasma viscosity, neutral particle reflection and the implemented minimum heat and particle fluxesare also discussed. The comparison of the simulations with the experimental results is given insection 5. The dependencies of the halo current profile on different parameters is researched insection 6. Finally the conclusions are summarized in section 7. The implicit non-linear code JOREK solves the extended magneto-hydrodynamic (MHD) equationsin realistic tokamak geometries including open field line regions [11, 8]. In JOREK the poloidalplane is discretized with quadrilateral bicubic Bezier elements using an isoparametric mapping [5].Fourier harmonics are used to decompose the toroidal direction and for the time stepping typicallythe Crank-Nicholson or the BDF2 Gears scheme are used. The code contains a hierarchy of differentMHD models with various extensions. The MHD model used for this work includes an equationfor fluid neutral particles ∂ A ∂t = v × B − η J − ∇ Φ , (1) ρ ∂ v ∂t = − ρ v · ∇ v − ∇ p + J × B + ∇ · τττ − ( ρρ n S ion − ρ α rec ) v , (2) ∂ρ∂t = −∇ · ( ρ v ) + ∇ · (DDD ∇ ρ ) + ( ρρ n S ion − ρ α rec ) , (3) ∂ρ n ∂t = ∇ · (DDD n ∇ ρ n ) − ( ρρ n S ion − ρ α rec ) , (4) ∂p∂t = − v · ∇ p − γp ∇ · v + ∇ · ( κκκ ∇ T ) + ( γ − ηJ − ξ ion ρρ n S ion − ρρ n L lines − ρ L brem (5)with the following reduced MHD ansatz for the magnetic field ( B ) and the mean plasma velocity( v ) B = ∇ ψ × ∇ φ + F ∇ φ, (6) v = − R F ∇ Φ × ∇ φ + v k , (7) here ψ is the poloidal magnetic flux and F = RB φ is a constant representing the main reducedMHD assumption, which is that the toroidal field is fixed in time. Note however that poloidal cur-rents are not fixed in time and evolve according to the current conservation and momentum balanceequations, but their contribution to the toroidal field is neglected due to the large vacuum field.A recent benchmark with full MHD codes has demonstrated that halo currents can be describedaccurately in this way [13]. The quantities shown in equations (1)-(7) are the magnetic potential( A ), the ion density ( ρ ), the neutral density ( ρ n ), the total pressure ( p ), the total temperature( T ≡ T e + T i ), the electrostatic potential (Φ) and the current density ( J ). The parameters relatedto the neutral particle fluid are the ionization and recombination rates ( S ion , α rec ), the Deuteriumionization energy ( ξ ion ) and the line ( L lines ) and bremsstrahlung ( L brem ) radiation coefficients.However the charge exchange process is not included, which leads to momentum losses at lowtemperatures ( T ∼ −
10 eV) and will be included in future works. Other parameters in equations(1)-(7) are the plasma resistivity ( η ), the stress tensor( τττ ), the thermal conductivity and the par-ticle diffusion coefficients ( κκκ, DDD , DDD n ) and the ratio of specific heats ( γ ). The thermal conductivitycoefficients κκκ present a high anisotropy (i.e. κ k (cid:29) κ ⊥ ) while the particle diffusion coefficients arenormally isotropic. Note that the reduction of the equations is ansatz-based, conserves energy [6]and does not result from geometrical ordering assumptions. The effect of impurity radiation playsan important role in the disruption dynamics, however this effect is not included in the presentsimulations and it is left for future work.The resistive wall and the free-boundary conditions are included by coupling JOREK to theSTARWALL code [16, 9, 3]. The coupling is performed by solving the full vacuum domain witha Green’s function method and therefore the vacuum does not need to be meshed. STARWALLuses the thin wall assumption and discretises the wall with linear triangular elements. Althoughthe employed wall formalism is implemented for 3D thin walls including holes, the setup used hereis restricted to an axisymmetric wall. Boundary conditions
In this subsection, the boundary conditions used throughout this paper are explained. The free-boundary conditions for the poloidal magnetic flux are explained in detail in the references [9, 3].
Parallel velocity ( v k ) The employed boundary condition for the parallel velocity to the magnetic field is a special caseof the Chodura-Riemann condition [18] at the magnetic pre-sheath v k = g ( b n ) c s (8)where c s ≡ p γk B T /m i is the ion sound-speed, and g ( b n ) is a function of the normal projectionof the magnetic field into the wall ( b n ≡ B · n /B ). Here n is a unit vector which is perpendicularto the boundary and points outside the plasma domain. The function g is needed to ensure thations always flow towards the wall ( v k · n >
0) and that the variable v k has a smooth transition inthe regions where b n changes sign. Otherwise the condition ( v k · n >
0) leads to discontinuities in v k along the boundary that cannot be resolved (e.g. v k would jump from c s to − c s ). The chosenform of the g function is g ( b n ) = sign( b n ) (cid:20)
14 (1 + tanh(( | b n | − g ) /g ) − g (cid:21) (9)and for this work the chosen coefficients are ( g , g , g ) = (0 . , . , . g (0) = 0and g ≈ ◦ . Heat flux
The boundary condition for the heat flux is based on reference [19] q · n = γ sh n k B T e v k · n + q min (10)where n is the particle density and γ sh is the sheath transmission coefficient. In this work thevalue γ sh = 11 is chosen as it has been found in COMPASS through experimental observations forsteady state plasmas [10]. Note that the evolution of γ sh during disruptions is presently unknown.Since the boundary conditions for the parallel flow are such that v k = 0 when the field lines re tangential to the wall (e.g. limiter point), a minimum heat flux ( q min ) has been introduced.Otherwise, if q · n = 0 at the tangency points, the thermal energy would accumulate artificially atthe boundary. We choose the minimum heat flux as q min = γ sh n e k B T e c s sin( ϑ min ), where unlessexplicitly stated, ϑ min = 2 ◦ . The choice of a minimum value is also motivated by the reference[15], in which significant particle and heat fluxes are found due to significant radial diffusion ofenergy and particles from the plasma at grazing angles. Ion flux
The employed boundary condition for the ion flux implies that ions leave the plasma domain atthe parallel velocity Γ · n = n v k · n + Γ min (11)Similar to the heat flux boundary condition, we introduce a minimum particle flux to avoid artificialaccumulation of particles, Γ min = nc s sin( ϑ min ). Neutral particle flux
In order to ensure conservation of particles through the simulation the neutral particle flux is Γ n · n = − Γ · n (12)which implies that all ions leaving the boundary come back as neutral particles. Electrostatic potential
Basic physics of the plasma implies that, in the presence of a sheath between the plasma and thematerial surfaces in contact with it, the maximum plasma current density that can flow alongthe field line is limited by the ion saturation current. Therefore in the presence of sufficientlylarge voltages induced by the disruption dynamics, the halo current could be directly determinedby the ion saturation current. This hypothesis has been confirmed in recent COMPASS VDEexperiments [1]. This was achieved by measuring independently the current flowing into the divertorthrough grounded Langmuir probes (to the vacuum vessel) and the ion saturation current at nearbylocations with biased Langmuir probes. During disruptions, the parallel current density ( J k ) wassimilar to the measured ion saturation current ( J sat ). As we will show later, MHD models donot normally include this limitation of the current density and can violate it. For that reason weexplore a suitable boundary condition for the MHD code JOREK including that limitation. Theparallel current integrated over a field line has the following evolution in reduced MHD Z ba ηJ k dl = − (Φ b − Φ a ) − Z ba B φ R | B | ∂ψ∂t dl (13)where a and b are respectively the starting and ending points of the field line and dl is the lengthdifferential along the field line. The latter relation is found by integrating the parallel electric fieldalong the field line. For a given flux variation (e.g. caused by the decay of the plasma current) thecurrent density along the field line is limited by the potential at its end points (Φ b , Φ a ). Thereforethe current density can be controlled by applying boundary conditions for Φ. At the plasma sheath,analytical expressions relating Φ and J k that feature such a limit are well known from Langmuirprobe theory [19]. The implementation of a boundary condition based on the Langmuir formulationis numerically very challenging and thus we have used a simplified boundary condition to limit thecurrent to the ion saturation currentΦ = ( Φ w + α ( J k − J sat ) · n , if J k · n ≥ J sat · n Φ w , if J k · n < J sat · n (14)where α is a constant. The boundary condition is such that when the parallel current ( J k ) exceedsthe ion saturation current ( J sat ), a voltage is driven with respect to the wall in order to reduce J k .By choosing a large α it is ensured that the parallel current remains close to J sat . If J k · n < J sat · n the potential is set to the wall potential, which is chosen to be Φ w = 0 (the boundary acts as aperfect conductor in the poloidal direction).The boundary condition for the potential can be expressed in terms of the JOREK variables.Assuming force-free balance in the open field- line region ( J × B = 0), the parallel current density s J k = − j B /F with j ≡ − RJ φ being the JOREK variable. The ion saturation current is J sat = sign( b n ) e n c s b and therefore the situation J k · n = J sat · n corresponds to j = j sat ≡− sign( b n ) e n c s F / | B | . Finally the JOREK boundary condition isΦ = ( − αF ( j − j sat ) B · n , if J k · n ≥ J sat · n , if J k · n < J sat · n (15) Vorticity and current density
In order to avoid second order derivatives in the code, auxiliary variables such as the toroidalvorticity w φ = F ∆Φ and the toroidal current density j = ∆ ∗ ψ are introduced. The boundarycondition for the vorticity is simply w φ = 0 and j evolves according to Ampere’s law ( j = ∆ ∗ ψ ). Implementation
The heat and particle fluxes are implemented as natural boundary conditions of the finite elementmethod. When performing the weak form and the partial integration of the DDD and the κκκ terms,the latter expressions for the fluxes are replaced in the arising boundary integrals. When Dirichletboundary conditions are applied, a penalization method is used. The method can be interpretedas adding a boundary integral term to the governing equations. For example the following term isadded to the parallel flow equation I Zv ∗ ( v k − g ( b n ) c s ) dS (16)where v ∗ is a test function of the FEM, dS is the boundary surface and Z = 10 is an arbitrarilylarge value used to enforce the boundary condition. In this subsection, we describe the plasma state prior to the VDE simulation. JOREK solves theGrad-Shafranov equation in order to obtain an initial condition that is accurate for the employedfinite element method. In its free-boundary version, the computation of the equilibrium requiresto specify the pressure and the toroidal current profiles as well as the geometry and the currentsof the PF coils. An initial equilibrium was computed based on an EFIT reconstruction of theCOMPASS discharge t = 1090 . q axis ∼ I p = 300 kA). The currents flowing in the 5 independent circuits of the COMPASS PF coilsare ( I BR , I BV , I EFPS , I
MFPS , I
SFPS ) = (0 . , . , − . , − . , − .
42) kA/turn. The I BR circuit isused for feedback control on the vertical position and its current has been increased to 0.9 kA/turnin the simulation to obtain a more accurate matching of the vertical position of the magnetic axis( Z axis ).For the plasma computational domain a polar grid has been chosen (see figure 1 (c)). Theboundary of the grid matches the COMPASS first wall in the regions of main plasma-wall inter-action. Its simplicity facilitates the implementation of the boundary conditions but some regionsare not represented accurately. The grid features radial mesh accumulation at the computationalboundary in order to resolve large temperature gradients that may arise. The vacuum vessel (bluecurve in figure 1 (c)) is an axisymmetric version of the COMPASS wall and it is discretized inSTARWALL with triangular elements using the thin wall approximation. The toroidal resistanceof the wall is 40% smaller than the value given by COMPASS specifications ( R w = 0 .
63 mΩ) inorder to have a better matching of the vertical position evolution with the experiment (see nextsection).Once the initial condition is given, an axisymmetric time evolution simulation (of ∼ η and κ k withtheir corresponding Spitzer and Spitzer-Härm values. A current source term ( ηj ) is added to themagnetic flux equation in order to fix the initial current profile over time ( ∂ t ψ = ... + η ( j − j ). )b) c)d) Vacuum vesselFirst wallSeparatrixLP arrayPolar grid
Figure 1:
Plasma configuration obtained after the steady state run in JOREK for the COMPASS shot t = 1090 . igure 2: Electron temperature distribution in logarithmic scale (left) and electron density distribution inlinear scale (right) at the end of the steady state run. The separatrix is indicated by the black curve.In the core, the steady state is reached when the Ohmic heating and the conductive perpendiculartransport are balanced. The chosen κ ⊥ is such that the final thermal energy ( W th = 3 . n e , T e ) measurements intothe EFIT reconstruction assuming T i = 0 . T e . The neutral diffusion coefficient was taken fromJOREK simulations without charge exchange that were compared to SOLPS-ITER simulations[17]. The ion and neutral particle densities are in equilibrium when the diffusive transport and theionization and recombination processes are balanced. The final core profiles are shown in figure 1(a) and (b). With the chosen parameters, the electron density accumulates at the High Field Side(HFS) near the lower X-point (see Figure 2 right), which explains the increasing plasma densitytowards the separatrix seen in Figure 1 (a). The final plasma temperature profile is also shown inFigure 2 (left) and several parameters describing the final steady state equilibrium are presentedin Figure 1 (d). Parameter Value Description D .
57 m / s Isotropic particle diffusion coefficient D n
228 m / s Isotropic neutral particle diffusion coefficient κ ⊥ . × (m s) − Perpendicular thermal conductivity κ k = κ ( T e /T e ) / κ = 3 . × (m s) − Parallel thermal conductivity η = η ( T e /T e ) − / η = 4 . × − Ω m Resistivity( µ k , µ ⊥ ) (3 . , . × − kg / (m s) Parallel/perpendicular dynamic viscosity γ sh
11 Sheath transmission coefficient
Table 1:
Parameters used during the steady state run. For the temperature dependent parameters thereference temperature is T e = 1 keV. Note that except for η and κ k the coefficients are spatially constant. I p0 ) Figure 3:
Comparison of two JOREK runs with different wall resistances (blue and green lines) and theexperimental time traces (red). (Top) Current change in the I BR circuit. (Middle) Vertical position of themagnetic axis. (Bottom) Total plasma current inside the plasma domain. The vertical dashed lines indicatethe two times ( t = 1096 . t = 1096 . ime (ms) Z a x i s ( m ) Time (ms) I p ( M A ) Experiment ( J sat = enc s ) J sat = ( = 0) Time (ms) I , h a l o ( M A ) Figure 4:
Comparison of two simulations with a different J sat value used for the boundary condition given byexpression 15. (Top) Vertical position of the magnetic axis. (Middle) Total plasma current. (Bottom) Toroidalplasma current in the halo region. The figure shows the stabilizing effect of the halo currents and the largeinfluence of the boundary condition for Φ for the dynamics of the VDE. During the steady state run, vertical position control is achieved with a PID controller acting onthe I BR PF coil circuit. In the experiments a vertical kick is performed by applying a currentwave-form to the I BR circuit in order to force a downwards VDE. In the simulation the samePF coil current variation (∆ I BR ) as in the experiment is specified and leads to a similar verticalposition evolution (see figure 3). We remind the reader that we have reduced the wall resistanceby 40% in order to have a better matching with the experiment. The case with the COMPASSreal wall resistance is also present in figure 3 showing a faster plasma displacement modelled byJOREK than measured in the experiment. For the VDE phase, the employed parameters were notmodified with respect to the steady state phase. The main difference is that the current sourceterm is switched-off and the current is allowed to decay. Also the code is run in axisymmetricmode. The boundary condition given by expression 15 is tested in this subsection and compared to theusual boundary condition applied in MHD simulations. When the usual MHD boundary conditionfor the electrostatic potential Φ = 0 is applied, large halo currents are induced in the SOL andthe decay time of the plasma current is much longer than observed in the experiment (see greenand red lines of Figure 4). In this case the halo current has a very large width and presents anon-monotonic profile as indicated by Figure 5 b), in which, the toroidal current density is shownat a late time of the VDE. A broad halo width ( w halo ) decreases the total resistance of the haloregion ( R halo ∼ η halo /w halo ) and consequently the total plasma current decays at a slower rate.Since in our simulations the core and the halo regions are treated with the same formalism, the fullplasma domain acts as a conductor with resistivity η ( T e ). Due to the electric field created by themoving plasma core, currents appearing in the far SOL (Figure 5 b) can be self-sustained through ) b) Figure 5:
Toroidal current density at Z axis = − .
27 m for the same two simulations as in Figure 4 withdifferent boundary conditions for the electrostatic potential. (a) Simulation with J sat = enc s and Φ employedin expression 15. (b) Simulation with Φ = 0 as boundary condition (equivalent to use J sat = ∞ in expression15). The figure shows the large influence of the boundary condition for Φ on the distribution of the halocurrents.Ohmic heating (which increases T e and makes the region more conductive).The case with the boundary condition Φ = 0 leads to non-physical results since the normalcurrent density ( J · n ) largely exceeds the ion saturation current ( J sat · n ) at different locationsalong the plasma-wall interface (see Figure 6). In order to visualize these quantities along theboundary of the plasma domain, the polar coordinate θ w is introduced. The polar coordinatesystem is centered at ( R, Z ) = (0 . ,
0) m such that θ w = 0 corresponds to the outer midplane( Z = 0). The results shown in Figure 6 imply that the usual boundary condition (Φ = 0) allowslarge currents to flow at locations with very low plasma densities (e.g. at θ w / π = 0 . J sat ≡ enc s ,the halo current is mainly constrained by the particle density profile ( n ) which evolves accordingto ionization, recombination, diffusion and convection. In practice, this implies that the halocurrent cannot be induced in the far SOL regions where the particle density is low. The lattercan be observed in Figure 5 a) as well as in Figure 7. As a consequence, the total induced halocurrent is much smaller than in the case without current limitation (see blue and green lines ofFigure 4). The halo current width is also severely reduced (see Figure 5 a)) and thus the totalplasma current decays at a faster rate, which is closer to experimental observations. Therefore theperformed simulations show that the boundary condition for the electrostatic potential (Φ) playsa very significant role for both plasma dynamics and distribution of the current profile in the haloregion (see Figure 7).The boundary condition shown in equation 15 provides an upper limit only for the positivecurrent ( J · n > p m i /m e ≈
61 larger thanthe ion saturation current. The limit to negative current in the LFS is not driven by the electronsaturation current but by the limit to the ion saturation current on the HFS and electric charge .5 0.6 0.7 0.8 0.9 1.0 w /2 J n ( M A / m ) J nJ sat n Figure 6:
Normal current and ion saturation current density at the wall-plasma interface as function of thepolar angle ( θ w ) for the simulation with Φ = 0 at the time when Z axis = − .
27 m. The figure shows that theparallel current density becomes larger than the ion saturation current for the usual MHD boundary condition,violating the experimentally found physical limit.conservation of the plasma, which leads to a similar negative current to J sat on the LFS despitethe Φ = 0 boundary condition there.Note that in the green shaded region of Figure 7 the halo current is somewhat larger than theion saturation current in the simulations. Unfortunately, the current could not be limited nearthe contact point or limiter point (green region) for the following reasons. Around the contactpoint (defined by B · n = 0), large gradients of the electrostatic potential are formed along theboundary ( ∇ tan Φ). The latter is due to the fact that at exactly the contact point Φ = 0, and inits vicinity, large voltages are applied in order to limit the parallel current (since B · n = 0). Sincethe parallel current density can locally increase due to these gradients via the local Ohm’s law( J k ∝ − B · ∇ Φ /η ), an unstable situation can occur around this point. In this case, an increaseof J k near the limiter point raises B · ∇ Φ, which in turn increases J k even further. To avoid thisissue, the following boundary condition is imposed to control the gradient in Φ along the boundaryΦ = − D lim ∂ Φ ∂l (17)where D lim is the distance in real space to the contact point and l is the length along the boundary.The latter boundary condition is applied up to 6 cm away from the limiter point, which is theminimum stable distance that has been found by a trial and error analysis. The effect of thisboundary condition on the electrostatic potential can be observed in the green region of figure7. With this setup, it has been found that for the latter phases of the VDE (e.g. Figure 7),the maximum normal current density ( J · n | max ) can exceed the maximum ion saturation current( J sat · n | max ) in this region by only up to ∼
30% of its value.Finally we investigate the origin of the non-monotonic current profile for the case with Φ = 0.Assuming spatially constant T e , J k and n and that the Ohmic heating and the parallel conductionare dominant, the evolution of the temperature in a SOL flux tube is ∂T e ∂t ≈ ηJ k n − c sh T / e A tube V tube = η T − / E k n − c sh A tube V tube ! T / e (18)where V tube is the volume of the flux tube, A tube is the cross-sectional area of the flux tube at itsend points, c sh ≡ γ sh p γk B /m i and we have used that E k = ηJ k . The latter equation showsthat for a given E k the flux tubes with smaller density and with smaller A tube /V tube ratios are themost efficient ones for increasing T e . The flux tubes in the upper part of the plasma domain arethe most efficient driving temperature because of the low density and the low A tube /V tube ratio.That is the reason why the current is preferentially induced on those regions and a non-monotonicprofile is observed in Figure 5 b). igure 7: (Top) Normal current and ion saturation current density at the wall-plasma interface as functionof the polar angle ( θ w ) for the "Base" VDE simulation at the time when Z axis = − .
27 m. The limiter point isindicated by the green line and the green area is the region where linear boundary conditions are applied forΦ. The purple dashed line is − J sat · n at the Low Field Side (see text for the physical limitations to negativecurrent in the LFS). (Bottom) Electrostatic potential as a function of θ w . /2 T e ( e V ) min = 2 min = 1 w /2 n e ( m ) Figure 8: T e and n e profiles along the wall for two different ϑ min values used in equations (10) and (11) at I p = 0 . I p . q min and Γ min In order to test the influence of q min and Γ min , two different simulations were run with ϑ min = 1 ◦ and ϑ min = 2 ◦ . The influence of ϑ min on the vertical position and current decay is not significant,but the T e near the contact point can be increased by 25% when ϑ min = 1 ◦ (see figure 8). The choice of the plasma perpendicular viscosity is found to play a very important role in thelimitation of J k . Near Alfvenic equilibrium and for cold plasmas, the current density is parallelto the magnetic field lines ( J × B = 0). In the absence of charge accumulation ( ∇ · J = 0), thecurrent density also satisfies B · ∇ ( J k /B ) ≈ F B · ∇ j = 0 (19)which implies that the quantity j ≡ RJ φ is approximately constant along the magnetic field-lines.When the boundary condition for Φ is applied, the current density at the boundary is locallylimited by J sat . This creates a gradient of current along field-lines leading to the excitation ofAlfven waves which restore B · ∇ j = 0. For time-scales longer than the Alfven time, gradientsalong the field-lines must vanish implying that J k ≤ J sat | wall . In the presence of unphysicallylarge fluid viscosities, a balance between the viscous force and the current gradient force can beestablished ( B · ∇ j ≈ µ ⊥ ∆ w ). Such a balance can prevent the excitation of Alfven waves and thehomogenization of the current along the field-lines leading to J k > J sat | wall . A case with largeviscosity is presented in figure 9 (left) showing the strong effect of the viscosity in the SOL currentprofile and the existence of large current density gradients along field-lines. Therefore specialcare must be taken with the viscosity choice when applying this type of boundary conditions. Toprevent unphysical results in our simulations, the value of µ ⊥ has been chosen small enough toavoid gradients of current density and such that a smaller value does not influence the plasmadynamics. In this subsection we discuss the large effect of neutral particles on the electron density and on theion saturation current. The VDE simulation was restarted from the end of the steady state run andthe neutral particle flux was to zero ( Γ n · n = 0). For that case, the electron density and the ionsaturation current is compared to the Base simulation at a latter phase of the VDE in Figure 10.When the neutral particle flux is set to zero, the plasma ions leave the computational domain atthe sound speed and do not come back as neutral particles. Therefore the boundary acts as a large igure 9: Current density in the far SOL in (left) simulation with µ ⊥ = 1 . × − kg m/s and (right)simulation with µ ⊥ = 1 . × − kg m/s. The viscosity coefficients are spatially constant. The white contoursindicate three different flux surfaces. When a large unphysical viscosity is present, current gradients along thefield lines can be sustained. w /2 J s a t n ( M A / m ) w/ neutral fluxw/o neutral flux w /2 n e ( m ) Figure 10:
Normal ion saturation current density (top) and electron density (bottom) at the wall-plasmainterface as function of the polar angle ( θ w ) for the time when Z axis = − .
27 m. The blue line corresponds tothe Base VDE simulation and the red line to a VDE case in which the neutral particle flux has been set to zero( Γ n · n = 0). The figure shows the importance of including neutral particles to calculate the electron particledensity and the ion saturation current at the plasma-wall interface. igure 11: Normal current density, electron temperature and normal heat flux on the divertor at (left) t = 1096 . I p = 0 . I p = 0 .
21 MA ( t = 1096 .
40 ms in the experiment and t = 1096 . q k ≈ T e | J k | and γ sh = 11. Theshaded region in green corresponds to the area where the J sat limitation of the current density is not effectiveand a linear boundary condition is applied for the electrostatic potential as explained in section 4.1.ion sink that leads to plasma densities lower than the Base simulation by two orders of magnitude(see Figure 10). As a consequence, the ion saturation current density is strongly reduced and halocurrents cannot be induced during the VDE. The latter demonstrates that neutral particles mustbe included when using sheath boundary conditions in order to a obtain a self-consistent evolutionof the plasma density and the ion saturation current at the plasma-wall interface. For the COMPASS experimental campaign on VDEs and current flows, two arrays of rooftop-shaped Langmuir probes (LPs) and one array of Ball-pen probes (BPPs) [2] were used to measurethe parallel current density and the electron temperature at the lower divertor. The location of theLPs array along the divertor is indicated in figure 1 (c). Since during COMPASS disruptions it isfound that J k ≈ J sat , the experimental heat-flux can be estimated by q k = γ sh T e J sat ≈ γ sh T e | J k | .Note however that this estimation is given for reference and that γ sh is assumed to be 11 while itcould vary significantly during the disruption since strong halo currents circulate in the SOL. Infigure 11, the JOREK values are compared to the experimental results along the boundary at twodifferent times.After a downwards motion of 10 cm and before the decay of the plasma current ( t = 1096 . R lim = 0 . J k = 0 ( R lim = 0 .
499 m). TheJOREK current profile is broader in the HFS than in the experiment, but due to the short widthof the LP array it is not possible to conclude the same at the LFS. The large power flux in theexperiment compared to JOREK is most likely due to the fact that the thermal quench is takingplace at 1096.0 ms in the experiment (as indicated by H α measurements) while such a quench isnot included in this axisymmetric JOREK simulations. ime (ms) Z a x i s ( m ) halo × 1 halo × 3 halo × 6 Time (ms) I p ( M A ) Time (ms) I , h a l o ( M A ) Figure 12: (Left) Vertical position of the magnetic axis, total toroidal plasma current and toroidal halocurrent in the plasma as function of time. (Right) Normal ion saturation current density, electron temperatureand normal heat flux on the divertor at I p = 0 . I p = 0 .
21 MA as function of the major radius. The lines withdifferent colors correspond to simulations with different factors multiplying the plasma resistivity in the haloregion. Note that in the upper right figure the absolute value of the experimental current density (red line) isplotted for reference.In figure 11 (right) the same quantities are compared at the time where the total current reaches70% of its pre-disruptive value in the experiment and in the simulation. Excellent agreement isfound on the location of the limiter position ( R lim = 0 .
50 m). This particular time point is chosenbecause the current densities reach their maximum values ∼ . The experiments indicatethat at this point, the temperatures decay to ∼
10 eV and that the heat-flux is significantly reduced.For the JOREK Base simulation, a thermal quench was not triggered artificially and cannot occurself-consistently due to the axisymmetric simulation setup, thus at this point, the plasma core stillcontains 28% of its pre-disruptive thermal energy. Due to the large pressures around the limiterpoint, the JOREK results show large heat fluxes compared to the experiments. However, where J k is limited by J sat in the simulation, similar current densities and temperatures are found, althoughthe JOREK simulations show that T e falls with a smaller decay length. The influence of severalsimulation parameters onto the results is shown in the next section. In this section we study how the different parameters shown in table 1 influence the profiles at thedivertor region. During all simulated cases, the parallel current was always found to be limited bythe ion saturation current (except for R ∈ [ R lim − . , R lim + 0 . J sat profile. The plasma resistivity dictates the diffusion of the current density and influences the energy balancethrough Ohmic heating. In figure 12 we study the influence of the resistivity in the halo region byincreasing it by a given factor only in that domain. Such an exercise can be viewed as a scan ofthe effective charge ( Z eff ) in the halo region. In this case the resistivity is modified in both the ime (ms) Z a x i s ( m ) ohmichalo × 0 ohmichalo × 1 ohmichalo × 3 ohmichalo × 6 Time (ms) I p ( M A ) Time (ms) I , h a l o ( M A ) Figure 13: (Left) Vertical position of the magnetic axis, total toroidal plasma current and toroidal halocurrent in the plasma as function of time. (Right) Normal ion saturation current density, electron temperatureand normal heat flux on the divertor at I p = 0 . I p = 0 .
21 MA as function of the major radius. The lineswith different colors correspond to simulations with different factors multiplying the Ohmic heating term in thehalo region without modifying current diffusion. Note that in the upper right figure the absolute value of theexperimental current density (red line) is plotted for reference. ime (ms) Z a x i s ( m ) P ohmiccore × 0 P ohmiccore × 1 Time (ms) I p ( M A ) Time (ms) I , h a l o ( M A ) Figure 14: (Left) Vertical position of the magnetic axis, total toroidal plasma current and toroidal halocurrent in the plasma as function of time. (Right) Normal ion saturation current density, electron temperatureand normal heat flux on the divertor at I p = 0 . I p = 0 .
21 MA as function of the major radius. The blackline corresponds to the base simulation and the yellow line to a case with the core Ohmic heating switched off.Note that in the upper right figure the absolute value of the experimental current density (red line) is plottedfor reference.induction equation ( ηJ k term) and in the pressure equation ( ηJ k term). The time traces in Figure12 (left) indicate that increasing η halo causes the vertical plasma motion and the current decay totake place earlier in time, although the decay rates do not change significantly. Also the toroidalhalo current has a weak dependence on η halo . The profiles at the time when I p = 0 . I p (Figure12 right) show that the increase in η halo flattens the J sat profile as well as the temperature andheat flux profiles.In order to distinguish between the effect of the resistive diffusion and the effect Ohmic heating,a set of cases were run where only the Ohmic heating term was multiplied by different factors inthe halo region. These cases are presented in Figure 13, where it can be seen that the Ohmicheating alone can effectively increase T e and q in the halo region. The increased Ohmic poweracts as an additional heating source that is able to increase the plasma temperature and due tothis also the pressure. For the factor range 1-6, the effect of the Ohmic power on the J sat profileis very weak, which means that the halo pressure is due mostly through a temperature increaseand thus at an almost constant particle flux (note that J sat ∼ p/ √ T ). However when neglectingcompletely the Ohmic heating in the halo region (blue line), the J sat magnitude is approximatelyreduced by a factor of 2. Note that in this case, the total toroidal halo current is not rigidly limitedto J sat in the ± J sat profile for Z eff ≥ ime (ms) Z a x i s ( m ) base D n × 0.5 D n × 3 D × 3 Time (ms) I p ( M A ) Time (ms) I , h a l o ( M A ) Figure 15: (Left) Vertical position of the magnetic axis, total toroidal plasma current and toroidal halocurrent in the plasma as function of time. (Right) Normal ion saturation current density, electron temperatureand normal heat flux on the divertor at I p = 0 . I p = 0 .
21 MA as function of the major radius. The legendscorrespond to a case with D n increased by a factor of 3 in all the plasma (green line), a case with D n decreasedby a factor of 2 (yellow line), a case with the ion diffusion increased by a factor of 3 (blue line), the basesimulation (black line) and the experimental results (red line). Note that in the upper right figure the absolutevalue of the experimental current density (red line) is plotted for reference.is much smaller than the power directly deposited in the halo by Ohmic heating. For example forthe Base simulation at t = 1096 .
38 ms (70% of I p ), the power deposited by Ohmic heating in thecore is 1.1 MW and the power directly deposited in the halo region by Ohmic heating is 7.8 MW. The particle diffusion coefficients can directly affect the distribution of the plasma density andtherefore can modify the ion saturation current density. In Figure 15 the standard simulation isrepeated with different ion and neutral particle diffusion coefficients. When increasing the neutralparticle diffusion coefficient by a factor of three (yellow line), the plasma density in the far SOL isreduced and the temperature increases, although the overall J sat profile is not strongly modified.A larger D n tends to decrease the plasma density at the wall due to the increase of the neutralmean free path (i.e. the ion particle source locally decreases as the neutrals are ionized at a furtherdistance from the wall). Decreasing this coefficient by a factor of two shows a small effect in the J sat profile, with a weak tendency to be increased in the far SOL. The increase of the ion diffusioncoefficient by a factor of three (green line) has a larger effect on the J sat profile. In this case theaugmented D can diffuse the plasma density at the wall and decrease J sat by about a factor of 2.Note that temperature does not increase accordingly and that finally the heat flux is reduced awayfrom the contact point (e.g. at R = 0 .
45 m). The reason is that a smaller J sat limiting the currentflow away from the contact point, obliges the halo current to be re-induced near the LCFS and todeposit its associated Ohmic heating in that region instead.For steady state plasmas, the ratio between the perpendicular and the parallel heat conductioncoefficients typically determines the width of the heat flux profile in the SOL. However for theconsidered VDE case, increasing the κ ⊥ /κ k ratio by a factor of three (either by modifying κ ⊥ or κ k ) does not increase the width or the magnitude of the heat flux and the J sat profiles (see blueand green lines of Figure 16). This indicates that during VDEs the SOL width is determined by ime (ms) Z a x i s ( m ) base/3× 3× 3 Time (ms) I p ( M A ) Time (ms) I , h a l o ( M A ) Figure 16: (Left) Vertical position of the magnetic axis, total toroidal plasma current and toroidal halocurrent in the plasma as function of time. (Right) Normal on saturation current density, electron temperatureand parallel heat flux on the divertor at I p = 0 . I p = 0 .
21 MA as function of the major radius. The legendscorrespond to cases with κ k decreased and increased by a factor of 3 (yellow and green curves), a case with κ ⊥ increased by a factor of 3 (blue line) the base simulation (black line) and the experimental results (red line).Note that in the upper right figure the absolute value of the experimental current density (red line) is plottedfor reference. ther processes such as the resistive current diffusion shown in the previous section. However theincrease of the parallel heat conduction by a factor of three (yellow line of Figure 16) is able toreduce the J sat profile more significantly. In this case the divertor temperature near the limiterpoint raises as the upstream temperature is conducted faster along the field lines. As the LCFSregion is hotter and therefore more electrically conducting, a significant fraction of the toroidalhalo current is preferentially induced near the LCFS and deposits its associated Ohmic heating inthat region, leading to the decrease of the J sat magnitude at a further distance in the SOL (e.g.at R = 0 .
45 m).
In this paper a fully self-consistent model for halo currents has been presented and explored numer-ically with the JOREK code. The reduced MHD model includes an equation for neutral particlesand a set of advanced boundary conditions. A boundary condition for the electrostatic potential(Φ) has been implemented in order to limit the halo current density ( J ) to the ion saturation cur-rent ( J sat ∼ n √ T ) as expected from basic sheath plasma physics and observed in recent COMPASSexperiments [1]. The boundary condition for Φ requires the implementation of sheath boundaryconditions for the parallel velocity, the ion density, the neutral density and the temperature in or-der to obtain a self-consistent evolution of J sat at the plasma-wall interface. The choice of realisticplasma parameters (e.g. Spitzer η ( T ) and Spitzer-Härm κ k ( T )) together with this model extendsthe status of halo current simulations beyond the present state of the art.The limitation of the halo current density to J sat is found to play a key role. The presentedaxisymmetric simulations of COMPASS VDE experiments show that with typical MHD boundaryconditions (Φ = 0), the current density becomes much larger than the ion saturation currentleading to non-physical results. When the boundary condition for the electrostatic potential isapplied during the VDE, the halo current density is always limited by the ion saturation current J ≈ J sat as observed in experiments, implying that the halo current density is mainly given by theparticle density. In this case both the halo width and total halo current are reduced in comparisonto the case in which the current density is not limited. The implementation of such boundarycondition requires a special treatment at the plasma-wall tangency point for numerical stability. Inaddition, the MHD simulations must be performed with a sufficiently low plasma viscosity in orderto minimize current density gradients forming along the open field lines. The implementation ofneutral particles is also found to be necessary when using sheath boundary conditions, otherwisethe plasma-wall interface acts as a large ion sink that leads to very low particle densities and halocurrent densities.Axisymmetric simulations of the COMPASS VDE discharge t = 1096 . J sat has been studied. In this respect the resistivity in the haloregion flattens the J sat profile through resistive diffusion, leading to an increased halo width. Thedependencies of J sat on the neutral particle diffusion and on the heat conduction coefficients areweaker than for the resistivity. However J sat is also sensitive to the ion diffusion coefficient, whichcan reduce its magnitude considerably by reducing the plasma density at the boundary. In allthe performed parametric studies the total halo current is only weakly affected by the specificparameter choice whether it affects or not J sat . This implies that the current flow limit by J sat forthe halo current density has a stronger impact on the halo width than on the total halo currentwith the latter being determined by VDE dynamics. ssuming the same J sat limitation for halo currents during ITER VDEs, the ITER halo widthwill scale like I p / ( q J sat ) since present experiments indicate that the total halo current scaleswith I p /q [12, 7]. Future simulations will be conducted with the presented model in order todetermine whether the latter scaling is to be expected and how J sat evolves for ITER disruptions. Acknowledgements
This work was supported by the ITER Monaco Fellowship Program. ITER is the Nuclear Fa-cility INB no. 174. This paper explores physics processes during the plasma operation of thetokamak when disruptions take place; nevertheless the nuclear operator is not constrained bythe results presented here. The views and opinions expressed herein do not necessarily reflectthose of the ITER Organization. The simulations presented here have been performed usingthe Marconi-Fusion supercomputer. This work has been carried out within the framework ofthe EUROfusion Consortium and has received funding from the Euratom research and train-ing program 2014-2018 and 2019-2020 under grant agreement No 633053. The views and opin-ions expressed herein do not necessarily reflect those of the European Commission. The workwas co-funded by MEYS projects number 8D15001 and LM2018117. This work has been carriedout within the framework of the project COMPASS-U: Tokamak for cutting-edge fusion research(No. CZ.02 . . / . / . / / References [1] J. Adamek, F. Artola, A. Loarte, J. Cavalier, R. Pitts, J. Havlicek, E. Matveeva, M. Hron, andR. Panek. Halo current density limitation during VDE-induced disruptions in the COMPASStokamak. to be submitted , 2020.[2] J. Adamek, J. Seidl, J. Horacek, M. Komm, T. Eich, R. Panek, J. Cavalier, A. De-vitre, M. Peterka, P. Vondracek, J. Stöckel, D. Sestak, O. Grover, P. Bilkova, P. Böhm,J. Varju, A. Havranek, V. Weinzettl, J. Lovell, M. Dimitrova, K. Mitosinkova, R. De-jarnac, M. Hron, and and. Electron temperature and heat load measurements in the COM-PASS divertor using the new system of probes.
Nuclear Fusion , 57(11):116017, aug 2017. doi:10.1088/1741-4326/aa7e09 .[3] F. Artola Such.
Free-boundary simulations of MHD plasma instabilities in tokamaks . The-ses, Université Aix Marseille, Nov. 2018. URL: https://tel.archives-ouvertes.fr/tel-02012234 .[4] K. J. Bunkers and C. R. Sovinec. The influence of boundary and edge-plasma modeling incomputations of axisymmetric vertical displacement.
Physics of Plasmas , 27(11):112505, 2020. arXiv:https://doi.org/10.1063/5.0023604 , doi:10.1063/5.0023604 .[5] O. Czarny and G. Huysmans. Bézier surfaces and finite elements for MHD simulations. Journalof computational physics , 227(16):7423–7445, 2008. doi:10.1016/j.jcp.2008.04.001 .[6] E. Franck, M. Hölzl, A. Lessig, and E. Sonnendrücker. Energy conservation and numerical sta-bility for the reduced mhd models of the non-linear jorek code. arXiv preprint arXiv:1408.2099 ,2014.[7] R. Granetz, I. Hutchinson, J. Sorci, J. Irby, B. LaBombard, and D. Gwinn. Disruptions andhalo currents in Alcator C-Mod.
Nuclear fusion , 36(5):545, 1996. doi:10.1088/0029-5515/36/5/i02 .[8] M. Hoelzl, G. Huijsmans, S. Pamela, M. Becoulet, E. Nardon, F. Artola, and B. N. et al. TheJOREK non-linear extended MHD code and applications to large-scale instabilities and theircontrol in magnetically confined fusion plasmas.
Nuclear Fusion (in preparation) , 2020.[9] M. Hölzl, P. Merkel, G. Huysmans, E. Nardon, E. Strumberger, R. McAdams, I. Chapman,S. Günter, and K. Lackner. Coupling JOREK and STARWALL codes for non-linear resistive-wall simulations. In
Journal of Physics: Conference Series , volume 401, page 012010. IOPPublishing, 2012. doi:10.1088/1742-6596/401/1/012010 .[10] J. Horacek, J. Adamek, M. Komm, J. Seidl, P. Vondracek, A. Jardin, C. Guillemaut, S. Elmore,A. Thornton, and K. J. et al. Scaling of L-mode heat flux for ITER and COMPASS-U divertors,based on five tokamaks.
Nuclear Fusion , 60(6):066016, 2020. doi:10.1088/1741-4326/ab7e47 .
11] G. Huysmans and O. Czarny. MHD stability in X-point geometry: simulation of ELMs.
Nuclear fusion , 47(7):659, 2007. doi:10.1088/0029-5515/47/7/016 .[12] P. Knight, G. Castle, A. Morris, A. Caloutsis, and C. Gimblett. Analysis of verticaldisplacement events and halo currents in COMPASS-D.
Nuclear Fusion , 40(3):325–337,mar 2000. URL: https://doi.org/10.1088%2F0029-5515%2F40%2F3%2F304 , doi:10.1088/0029-5515/40/3/304 .[13] I. Krebs, F. Artola, C. Sovinec, S. Jardin, K. Bunkers, M. Hoelzl, and N. Ferraro. Axisym-metric simulations of vertical displacement events in tokamaks: A benchmark of M3D-C1,NIMROD, and JOREK. Physics of Plasmas , 27(2):022505, 2020. doi:10.1063/1.5127664 .[14] M. Lehnen, K. Aleynikova, P. Aleynikov, D. Campbell, P. Drewelow, N. Eidietis, Y. Gas-paryan, R. Granetz, Y. Gribov, N. Hartmann, et al. Disruptions in ITER and strategiesfor their control and mitigation.
Journal of Nuclear Materials , 463:39–48, 2015. doi:10.1016/j.jnucmat.2014.10.075 .[15] G. Matthews, S. Fielding, G. McCracken, C. Pitcher, P. Stangeby, and M. Ulrickson. Investiga-tion of the fluxes to a surface at grazing angles of incidence in the tokamak boundary.
PlasmaPhysics and Controlled Fusion , 32(14):1301, 1990. doi:10.1088/0741-3335/32/14/004 .[16] P. Merkel and E. Strumberger. Linear MHD stability studies with the STARWALL code. arXiv preprint , 2015. URL: https://arxiv.org/abs/1508.04911 .[17] S. Smith, S. Pamela, A. Fil, M. Hölzl, G. Huijsmans, A. Kirk, D. Moulton, O. Myatra,A. Thornton, and H. W. and. Simulations of edge localised mode instabilities in MAST-u super-x tokamak plasmas.
Nuclear Fusion , 60(6):066021, may 2020. doi:10.1088/1741-4326/ab826a .[18] P. C. Stangeby. The Bohm–Chodura plasma sheath criterion.
Physics of Plasmas , 2(3):702–706, 1995. doi:10.1063/1.871483 .[19] P. C. Stangeby et al.
The plasma boundary of magnetic fusion devices , volume 224. Instituteof Physics Pub. Philadelphia, Pennsylvania, 2000., volume 224. Instituteof Physics Pub. Philadelphia, Pennsylvania, 2000.