Numerical investigation of disruption characteristics for the snowflake divertor configuration in HL-2M

Cold and hot vertical displacement events (VDEs) are frequently related to the disruption of vertically-elongated tokamaks. The weak poloidal magnetic field around the null-points of a snowflake divertor configuration may influence the vertical displacement process. In this paper, the major disruption with a cold VDE and the vertical disruption in the HL-2M tokamak are investigated by the DINA code. In order to better illustrate the effect from the weak poloidal field, a double-null snowflake configuration is compared with the standard divertor (SD) configuration under the same plasma parameters. Computational results show that the weak poloidal magnetic field can be partly beneficial for mitigating the vertical instability of the plasma under small perturbations. For major disruption, the peak poloidal halo current fraction is almost the same between the snowflake and the SD configurations. However, this fraction becomes much larger for the snowflake in the event of a hot VDE. Furthermore, during the disruption for a snowflake configuration, the distribution of electromagnetic force on a vacuum vessel gets more non-uniform during the current quench.


Introduction
Disruptions that might lead to great damage of the machine are a major safety issue in tokamak operation. During the disruption, large thermal loads on the plasma facing comp onents can cause an erosion of the materials. Meanwhile large electromagnetic (EM) forces on the vacuum vessel (VV) and other invessel conducting components (e.g. blanket modules, first wall (FW) panels and invessel coils) can lead to the displacement or distortion of these structures. Runaway electrons generated during the current quench (CQ) of disruptions are particularly damaging for the FW, poten tially resulting in the deep melting of main chamber armors. Therefore, assessing the magnitude of these loads and dam ages, via computational tools, is of paramount importance for the engineering design of a future machine. Predictive model ling is even more crucial for reactorscale devices that cannot afford a disruption.
For a verticallyelongated tokamak, disruptions are frequently related to the cold vertical displacement event (VDE) and the hot VDE [1]. The released thermal energy at the beginning of a major disruption (MD) causes the

Numerical investigation of disruption characteristics for the snowflake divertor configuration in HL-2M
Cold and hot vertical displacement events (VDEs) are frequently related to the disruption of verticallyelongated tokamaks. The weak poloidal magnetic field around the nullpoints of a snowflake divertor configuration may influence the vertical displacement process. In this paper, the major disruption with a cold VDE and the vertical disruption in the HL2M tokamak are investigated by the DINA code. In order to better illustrate the effect from the weak poloidal field, a doublenull snowflake configuration is compared with the standard divertor (SD) configuration under the same plasma parameters. Computational results show that the weak poloidal magnetic field can be partly beneficial for mitigating the vertical instability of the plasma under small perturbations. For major disruption, the peak poloidal halo current fraction is almost the same between the snowflake and the SD configurations. However, this fraction becomes much larger for the snowflake in the event of a hot VDE. Furthermore, during the disruption for a snowflake configuration, the distribution of electromagnetic force on a vacuum vessel gets more nonuniform during the current quench.
Keywords: disruption, VDE, snowflake divertor (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. plasma pressure to rapidly decay. Large changes of plasma para meters lead to further, fast development of the ver tical instability, followed by an uncontrolled vertical dis placement. Subsequently, the CQ occurs when the plasma contacts the wall. This VDE, occurring after the thermal quench (TQ) due to the major disruption, is often termed the cold VDE [2,3]. Another type of VDE often occurs due to the failure of the vertical position feedback control system, excessive elongations, or even the large magneto hydro dynamic (MHD) events such as the edge localized modes. In such cases, small perturbations normally grow to a large amplitude due to nonlinear motions of electrons and ions, resulting in uncontrolled vertical displacement of the plasma column. Subsequently, the TQ and the CQ occur almost simultaneously. This scenario is called the hot VDE [4,5].
The HL2M tokamak is a mediumsized coppercon ductor machine, currently under construction in China. The device has a major radius of 1.78 m, with a plasma minor radius of 0.55-0.65 m. The design uses the advantage of demountable toroidal field (TF) coils to increase the exper imental flexibility in studying highperformance plasmas, thus exploring technical and engineering potentials toward ITER and future fusion reactors [6]. The poloidal field (PF) coils are placed inside the TF coils; this makes it fea sible for the machine to generate vertical highlyelongated plasmas. The study of VDEs is thus of crucial importance for HL2M.
The snowflake divertor (SF) configuration is favorable for mitigating the deposition of the heat loads. This configuration has been realized in some of the current tokamak devices, and has been considered in future fusion devices including HL2M. The key feature of a SF configuration is the pres ence of a large zone with a very weak poloidal magnetic field around the nullpoint [7,8]. What is presently not well understood, however, is the effect of such a weak PF field on the VDE process. The VDE modelling of the SF config uration, as well as the assessment of EM loads during the VDE, is thus not only important for the subsequent upgrade of HL2M, but also for the design of future fusion reactors. In this paper, we shall utilize the predictive mode of the DINA code, to investigate the disruption in both cold and hot VDE regimes, for the HL2M plasmas in a SF configuration. The (up-down) symmetric PF coil system in HL2M has the capa bility of generating doublenull configurations. A doublenull SF configuration is thus chosen for the present investigation. In order to illustrate the disruption characteristics of the SF better, a comparative study is performed for a standard divertor (SD) configuration with the same plasma parameters. We shall focus on the temporal behaviors of the main plasma para meters, the peak poloidal halo current fraction, and the electro magnetic forces as well as the force distribution along the VV during the disruption. Section 2 introduces the DINA model as applied to the HL2M tokamak. Sections 3 and 4 report the details of numericallysimulated MD, and the hot VDEs using the DINA code, respectively. A summary is drawn in section 5.

DINA model applied to HL-2M
The DINA code [9] is a mature axisymmetric timedepen dent equilibrium code, which has been successfully applied to model DIIID [10], JT60U [11], ASDEXU [12], MAST [13] and TCV [14] plasmas, as well as used for a predictive study in ITER [15,18]. The code is capable of investigating disruptions well. When applied to HL2M, the structure model includes 16 PF coils, one central solenoid (CS) coil, and 88 magnetic probes as well as 52 loops. The VV is divided into 290 filaments. The divertor target plates, which are insulated along the toroidal direction, are generated by a limiter module of the code. The computational domain, cov ering a rectangular area from 0.9 to 2.6 m in the Rdirection and from −1.8 to 1.8 m in the Zdirection, is divided into 8 385 grid elements. The computation model, as shown in figure 1, is used to prepare the Green functions for analyzing the EM responses.
The predictive mode of the DINA code is used for the VDEs simulation. In order to demonstrate the prominent VDE characteristics, a highlyelongated plasma shape is considered in this work. Subject to the limitation of the vertical space of HL2M, it is easier to obtain such an elongated plasma in the simulation, by choosing a smaller minor radius (a = 0.55 m). Furthermore, in order to investigate the effect of the weak PF of the SF on the path of the VDE, a low triangularity is chosen (δ = 0.25 95 ) for the snowflake configuration. The plasma current is chosen as 1 MA, in favor of increasing the initial safety factor before the TQ. These parameters, as well as the common parameters designed for HL2M plasmas are shown in table 1. For the initial equilibrium, the plasma current density is chosen as follows where a k , b k are the fitting coefficients (k = 0-2), R is the major radius and R 0 is the geometrical center of the plasma. The eddy current, in each VV filament, is represented as the sum of up to five poloidal Fourier harmonics Here, θ i is the poloidal angle. The coefficients a k , b k , A m and B m , in the above equations (1)-(2), can be obtained by the singular value decomposition technique. Based on these representations, the initial equilibrium doublenull divertor configurations (the snowflake and the standard configura tions with the same global plasma parameters) are computed and are shown in figure 2. A twodimensional free boundary plasma equilibrium solver, in the presence of external poloidal magnetic fields, combined with circuit equations for the VV filaments and the poloidal flux diffusion equation, is utilized in order to obtain the converged magnetic flux sur faces. The onedimensional transport equations, including the ion and electron energy diffusion equations, as well as the particle diffusion equations, are then applied for com puting the plasmapressure profile for the next time step [9,18]. This how the evolution of the plasma configuration is implemented in the code.
In addition, a model of the halo current [19] is included in the DINA code. The physical model of the halo area expan sion during the plasma disruption is based on the assumption of either the conservation, or partial decrease, of the toroidal magnetic flux inside the poloidal surface S, combining the core plasma and the halo region. Because the toroidal magn etic flux is proportional to S, the scaling for the halo width w(t) in the DINA code can be obtained by the expression as shown in equations (3). Here S 0 and I p0 are the poloidal surface and the plasma current, respectively, before the start of the TQ. S(t, w) and ( ) I t w , p are the aforementioned quantities after the TQ. C is an adjustable coefficient for the decrease of the toroidal flux, which is equal to 6.0 in the present model. The poloidal and toroidal components of the halo current are determined by equations (4) and (5), respectively. This allows us to assess the magnitude of the halo current.

Investigation of the MD with a cold VDE
In the generic predictive mode of the DINA code, the one dimensional transport equations, and the poloidal flux Table 1. Main plasma parameters for the simulation: I p (the plasma current), R 0 (the geometric center), a (the minor radius), κ 95 (the elongation on the 95% magnetic flux surface), β p (the normalized poloidal pressure), l i (the internal inductance), δ 95 (the triangularity on the 95% magnetic flux surface), B t (the toroidal magnetic field on the geometric center).
Parameters  averaged on magnetic surfaces, are selfconsistently evolved. However, since the energy transport is highly abnormal during and after the TQ phase, an accurate prediction of the disruption onset is still difficult to achieve at present. Much experimental statistic work has been done on investigating the precursors of disruptions [16,17]. In the simulation, the trigger time of the TQ is artificially prescribed in the code based on the experience from experiments, without solving the transport equations [18].
In simulating the MD, we define an equilibrium state as the starting point (t = 0 ms). The onset time of the TQ is set before the vertical displacement appearing at 1.1 ms, and the dura tion time of the TQ is 1.0 ms (from 1.1 to 2.1 ms). During the TQ, the poloidal normalized pressure (β p ) quickly drops due to fast energy release. The internal inductance of the plasma increases to conserve the total poloidal magnetic flux. Large scale changes of plasma parameters, such as the β p and the internal inductance, gives rise to the fast increase of the ver tical instability. In the final stage of the TQ, an uncontrolled vertical displacement occurs and the plasma column moves downwards.
As shown in figure 3, accompanied by the downwards movement of the plasma, the edge safety factor q 95 and the edge elongation factor κ 95 drop quickly. The poloidal cross section shrinks rapidly. Meanwhile the internal inductance increases quickly; this is associated with the fact that the plasma current density profile becomes more peaked with time. By the time the plasma column contacts the divertor baf fles (at 2.1 ms), the scenario evolves into a limiter configu ration, with the plasma central current density reaching the maximum. Simultaneously the total vertical EM force on the VV also reaches the first peak value ( figure 4). Normally, the net EM force on the VV is negative with the downwards shift of the plasma. However the peak value of the EM force is positive in this simulation. This may be related to the peaked plasma current profile as shown in figure 4.
The CQ rate entirely depends on the plasma resistivity, which is determined by the electron temperature T e and the   effective charge Z eff . Here the T e in both the plasma core and the halo region are set to be 10 eV. The profile of Z eff is assumed to be flat, with a constant value of 1.7. The CQ can be further divided into two phases-the core plasma CQ phase and the halo CQ phase. As the core plasma current rapidly drops, the toroidal component of the halo current is driven in the scrape off layer (SOL). Meanwhile the fast shrinkage of the plasma cross section, determined by the force balance equation, causes the toroidal magnetic flux and the internal inductance to decay, further inducing the poloidal halo current [19][20][21] in the SOL region. As shown in figure 5, the toroidal and the poloidal comp onents of the halo current reach their maximum when the core plasma current entirely quenches. A comparison of the results during the MD shows that the temporal behaviors of the plasma parameters (figure 3), the maximum total vertical EM force on the VV (figure 4), and the poloidal halo current peak (figure 5), are almost the same between the SF and the SD configurations. The largescale perturbations during the TQ move the plasma so fast, that the weak PF of the snowflake configuration does not seem to have sufficient time to affect the temporal behavior before the plasma contacts the baffles.
The 290filament representation of the VV makes it feasible to investigate the spatial distribution of the vertical EM forces on the VV. The time corresponding to the two peaked values of the total vertical EM force are chosen as the reference here. Figures 6 and 7 do not show obvious differences between the SF and the SD configurations, in terms of the spatial distribu tion of the eddy currents. The magnetic configuration formed by the divertor coils dominates over the EM force distribution on the VV. The only visible effect is that the spatial distribution of the EM forces on the VV becomes much less uniform near the divertor baffles under the effect of the SF coils.

Investigation of the hot VDE
For the hot VDE simulation, the vertical feedback control system is turned off at t = 0 ms. Small perturbations from the plasma result in the development of the vertical instability. As  time goes onwards, the plasma column starts to move down wards. Compared with the standard configuration with the same κ 95 , the presence of the weak PF in the SF increases the elongation (κ x ) of the LCFS. Normally a higher κ x leads to a faster development of the vertical instability [22]. However, the simulated development seems to be relatively slow in the SF case, as shown in (a) of figure 8, for the time traces of the magnetic axis along the Z direction. Compared with the SD configuration, the uncontrolled vertical displacement defined by the criterion of Z > 0.05 m, occurs earlier in with the SF. On the other hand, it seems that the weak PF of the snowflake configuration may mitigate the development of the vertical instability before the VDE occurs. Figure 8(b) shows that, from the beginning of the VDE to the TQ onset, the VDE velocity of the plasma column gets relatively slower. This is also demonstrated in figure 9. The lower negative divertor coil current, which is used to provide the weaker poloidal field, partly restrains the downward movement of the plasma. The snowflake configuration thus seems to be more beneficial in mitigating the hot vertical displacement than the standard configuration.
During the hot vertical displacement, the edge safety factor q 95 , as well as the internal inductance (Li3), decreases with the shrinkage of the plasma cross section. The value of q 95 is thus used as the onset criterion for the TQ [11]. When the plasma column contacts the divertor baffles, a limiter configuration appears as shown in figure 10. The TQ is triggered when q 95 drops to less than 3.0 at the 16.9 ms. Here we still assume that the duration of the TQ is 1.0 ms. As shown in figure 11, the plasma cross section decreases during 16.9-17.1 ms quickly,  During this phase, we pay part icular attention to the peak poloidal halo current fraction and the EM force on the VV. These are critical factors for the engineering design of HL2M. As the safety factor decreases, the core plasma becomes increasingly unstable. Eventually the decay rate of the core plasma current, the shrinkage rate of the crosssectional area, and the vertical velocity of the plasma column reach their maximum almost simultane ously. This leads to a spike in the halo current. As shown in figure 12, the peak poloidal halo current fraction, which appeared at 18.7 ms in the SF case, reaches 57.2%, which is 24.4% higher than that which appeared at 15.6 ms in the SD case. The generation of the higher peak poloidal halo current in the SF seems to be mainly related to the faster decay rate of the core plasma current. It is known that the divertor coils, which are used to generate the weak PF in the SF, are applied to a negative toroidal current. Figure 13 shows that under the pushing effect of the negative current, the crosssectional area gets smaller (0.66 m 2 ) than that (0.72 m 2 ) in the SD VDE case when the CQ occurs (13.6 ms for the SD and 17.1 ms for the SF), which gives rise to a lower q 95 . The lower q 95 causes the core plasma to be vertically more unstable, leading to a faster quench of the core plasma, and finally to the larger halo spike.  The corresponding evolution of the configuration is shown in figure 14.
As shown in figure 15, two peaks form for the total down ward vertical EM force; one peak during the core plasma CQ phase and the other peak during the halo CQ phase. For the SD case, the maximal total vertical EM force occurs at the second peak. However, for the SF configuration, the maximal force appears at the first peak, with the whole VV suffering a net downward EM force up to 17.1 tons. The spatial dis tribution of the EM force on the VV is very similar to that simulated for the MD with a cold VDE, as shown in figures 6 and 7.

Summary
The disruption characteristics of a doublenull SF con figuration in the HL2M tokamak are numerically investi gated by utilizing the predictive mode of the DINA code.
In the modelled major disruption, a downward cold VDE is generated due to largescale perturbations during the TQ. Meanwhile the internal inductance increases quickly, corresponding to fast peaking of the plasma current pro file. The peaked plasma current profile may be responsible for the upward total vertical EM force in the final stage of the TQ. When the plasma touches the divertor baffles, the CQ occurs. The fastquenching core plasma current gives rise to a high halo current. The poloidal halo current frac tion reaches 33.5%. The maximum total vertical EM force generated in the CQ stage is downward, reaching 6.5 tons. Evolutions of the aforementioned parameters are almost the same as those of the MD event in the SD configuration, pro vided that the main plasma parameters are kept the same. It appears that, except for a more nonuniform spatial distribu tion of the EM force on the VV, no obvious effect from the weaker PF of the snowflake configuration is exerted during the MD process.

SD SF
However, during the hot VDE process, the weaker PF of the snowflake configuration seems to play a more important role. Compared to the standard configuration, the presence of a weaker PF can mitigate the vertical instability caused by small perturbations from the plasma. Moreover, the divertor coils with negative (opposite) currents, which are used to generate the weak poloidal field, can also slow down the vertical velocity of the plasma column even when the hot vertical displacement occurs. This may be beneficial for the vertical position control. We notice that the pushing effect of the negative divertor coil current causes a lower edge safety factor when the CQ occurs. This may make the plasma more unstable along the vertical direction, which eventually gives rise to a faster CQ of the core plasma (within 1.5 ms for the SF, and within 1.9 ms for the SD) and a higher halo current fraction. The peak poloidal halo current fraction reaches 57.2% in the hot VDE event, being 24.4% higher than that of the standard VDE case. Furthermore, the maximum total  vertical EM force reaches 17.1 tons. The spatial distribution of the EM force on the VV becomes more nonuniform near the divertor baffles, very similar to that observed during the MD case.