Beryllium melt instabilities and ejection during unmitigated current quenches in ITER

The dynamics of transient liquid beryllium flows induced on the ITER first wall during the current quench stage of unmitigated vertical displacement events are modelled by means of two-dimensional Navier–Stokes simulations. The study focuses on melt that is driven to the first wall panels’ chamfered edges, where free-surface instabilities are the most likely to be seeded. Beyond their impact on plasma-facing component damage, these instabilities potentially result in material ejection in the form of droplets, which may ultimately solidify into dust and accumulate in the vessel. Based on prior integrated numerical predictions of quenching magnetic equilibria, wall energy deposition and melt-related damage in a concrete worst-case disruption scenario, the simulations suggest that, although the liquid layer is significantly destabilized, only 5% of the total melt mass created on the wall surface is lost through ejection. This result can serve as a basis to refine the estimates of the real transient-induced beryllium dust inventory expected in ITER.


Introduction
Quantitative estimates of metallic dust generation during the operation of fusion reactors play a pivotal role in nuclear safety and licensing issues [1][2][3][4][5]. Two main dust sources are expected in ITER with its beryllium (Be) first wall and tungsten (W) divertor: the delamination of layered co-deposits growing on plasma-facing components (PFCs) over the course of normal operation, and the destabilization of surface melt pools created during off-normal high-energy transients. The * Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.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. latter mechanism is particularly appropriate in the case of unmitigated disruptions, whose occurrence is more likely in the early phases of ITER operation, when experience is first being gained with the planned ITER disruption mitigation system.
In this regard, it is important to distinguish between the thermal quench (TQ) and current quench (CQ) stages of plasma disruptions. During TQs, whose characteristic duration on ITER is expected to be a few milliseconds, liquid pools may form by flash melting and be traversed by strong eddy currents induced by the rapid changes in poloidal magnetic field [6]. The resulting Lorentz forces are directed towards the plasma, and hence can seed Rayleigh-Taylor instabilities, which are expected to develop quickly enough to cause melt loss [7]. Conversely, under CQ conditions, molten layers are sustained by plasma heat fluxes persisting for up to hundreds of milliseconds, and mainly driven by surface-parallel Lorentz forces due e.g. to halo currents [8,9].
Evidence of transient macroscopic PFC melting has been observed repeatedly in contemporary tokamaks [10][11][12][13][14]. Among these experimental results, traces of melt ejection from specially prepared misaligned W components in ASDEX Upgrade [12] and from Be upper dump plates in JET [13] during deliberately induced vertical displacement events (VDEs) are particularly noteworthy. In both cases, fluid-dynamics instabilities appear to be triggered by the presence of edges or similar geometric obstacles on the PFC surface. Whereas considerable progress has been achieved in modelling the bulk molten layer dynamics thanks to simulation tools such as the MEMOS-U shallow-water melt code [15][16][17], separate dedicated studies must be conducted in parallel to further the current understanding of ejection processes. Recent efforts on this front have focused on numerical solutions of the Navier-Stokes equations on reduced computational domains, benchmarked against published data [18][19][20] and successfully validated against observations of Be splashing in JET [13,21]. These studies have provided new insights on the destabilization mechanisms at play in fusion environments and have allowed a variety of analytical estimates and dimensionless criteria to be tested [4,5].
The purpose of this paper is to present the first simulation of Be melt instabilities in a concrete ITER CQ scenario. Disturbances in the liquid layer are caused by the presence of a 45 • chamfer at the edge of ITER first wall panels (FWPs). The Fluent R computational fluid dynamics software is used to predict the evolution of a flowing molten layer injected with prescribed space-and time-dependent profiles echoing the variations of plasma-wall interactions during the CQ. The main dynamical features of the flow are analysed qualitatively and quantitatively, with a focus on free-surface breakups and material ejection events, whose characteristics are of great importance to dust migration, survival and accumulation studies [22].

Input from prior simulations
The complete response of an ITER FWP to CQ loads involves characteristic lengths up to 1 m over several hundreds of milliseconds, while the dynamic instabilities potentially responsible for melt ejection and splashing are typically associated with distances of the order of 100 μm and sub-millisecond durations [4,21]. This wide scale separation imposes a composite modelling strategy wherein different numerical tools are used to address specific aspects of the PFC response and damage in an integrated fashion. In particular, the study presented in this paper, which focuses on the dynamics of Be melt flows at FWP chamfers, requires external input regarding the characteristic thickness, velocity, and temperature of the melt layers created during disruptions. This information was extracted from MEMOS-U macroscopic melt dynamics simulations of an unmitigated 15 MA/5.3 T upward VDE striking FWP #8, which has been identified as a worst-case scenario in terms of wall damage [9]. The evolution of the disrupting Three-dimensional view of the thermal loading on the ITER first wall from SMITER calculations of an unmitigated upward VDE (see [9] for more details) at the time of first melt edge crossing (t = 0). The colour bar gives surface heat flux values in W m −2 . The particular instance of FWP #8 simulated in this paper is labelled, and the white arrow indicates the poloidal chamfer that is crossed by Be melt. plasma was simulated by the DINA code [23], whose output was fed into the field-line tracing software SMITER [24] to reconstruct the heat flux and halo current density impinging on the PFCs, as exemplified by figure 1. These loads were then projected onto the surface of the rectangular cuboid domain used in MEMOS-U simulations to predict the formation and subsequent motion of Be melt.
According to MEMOS-U calculations, as a result of its acceleration by Lorentz forces, the liquid Be pool reaches the chamfered edge of the panel towards the end of the CQ. This instant will be referred to as t = 0 in all the following, so that the entire interaction between the quenching plasma and FWP #8 takes place from t −340 ms to t 60 ms.   These values are given in table 1 at 20 ms intervals until t = 60 ms, at which point MEMOS-U predicts that melt motion has essentially ceased. It should be mentioned here that the shallow-water equations currently implemented in MEMOS-U [17] neglect the normal stresses due to surface tension, potentially resulting in the formation of an unphysically large shock-like ridge at the leading edge of the moving melt layer, as visible in figure 2(a). Indeed, the friction force-of order μu/h 2 , where μ is the dynamic viscosity-can become very large near the edge of the liquid pool where h → 0. This leads to a local thickening of the layer into a ridge as fluid parcels initially located away from the edge, and hence experiencing weaker viscous damping, catch up with it. Letting L be the characteristic width of this melt pile-up phenomenon, i.e. the length scale across which h varies significantly, surface tension results in an effective body force of order γh/L 3 , where γ denotes surface tension [25,26]. This prevents the growth of the ridge beyond h ∼ Ca 1/3 L [26], where the capillary number Ca = μu/γ encodes the balance between viscosity and surface tension. In the MEMOS-U results of interest here, i.e. those relating to the edge of the molten layer first reaching the FWP chamfer, Ca 1/3 = 0.25 and L ∼ 1 cm, hence the observed h values which locally exceed a few millimetres should be regarded as spurious. Nevertheless, these numerical artefacts remain localized to a thin band of computational MEMOS-U cells and are not expected to have any major effect on the bulk dynamics of the liquid pool, where h/L 10 −3 Ca 1/3 . Therefore, the anomalous melt thickness datum at t = 0 is discarded and h(t = 0) is extrapolated backwards from h(t = 20 ms), assuming by default that it remains constant during that interval.

Main equations and material properties
Similarly to [20,21], Fluent is used to solve the twodimensional Navier-Stokes system for an incompressible twophase flow with heat transfer using the volume-of-fluid method [27]: where v is the velocity, α the metal volume fraction, ρ the mass density, p the pressure, g the gravitational field, T the temperature, ε the specific enthalpy and λ the thermal conductivity. In (3), the effects of surface tension are accounted for according to the continuum surface force model [28] via a term involving the algebraic curvature κ of the free surface, which is computed natively by Fluent from higher-order spatial derivatives of α. The enthalpy-porosity method [29] is employed, wherein the enthalpy-temperature relationship allows for temperature-dependent values of the volumetric heat capacity c v and accounts for solidification and melting via the specific latent heat of fusion ΔH fus and the molten fraction  [33] to accommodate the constant-density assumption) and thermal conductivity (dashed curve [34]); (b) surface tension (solid curve [31,32]) and dynamic viscosity (dashed curve [35]). β, which is implemented as a piecewise linear approximation of H(T − T melt ) [21], where T melt is the melting point and H is the Heaviside step function. The momentum source term F s in (3) is active when β = 0 and serves to freeze the solid metal in place [20,29]. Furthermore, no external energy sources are considered in (4): the effects of plasma heating are already accounted for in the MEMOS-U temperature output, while the data given in table 1 can be used to show that surface cooling fluxes due to thermal radiation and vaporization do not exceed 1% of the conductive heat flux. Beryllium is modelled according to the same data compendium implemented in MEMOS-U [17,30], with  [31][32][33]; its remaining properties are allowed to vary with temperature as plotted in figure 3. Although mass density is assumed to be constant, which implies that the mechanical effects of thermal expansion are neglected, the values used for the specific heat capacity c p are corrected so that c v (T) = ρc p (T) is consistent with the available literature [33]. The properties of the background gaseous phase-which is meant to represent the dilute gas or weakly ionized plasma in the vicinity of the FWP-are set to negligible values to ensure that it behaves as a ghost fluid that has no noticeable effect on the dynamics of the metal.

Initial and boundary conditions
The computational domain used in the simulation is sketched in figure 4. It is split into two main regions: a closed region bounded by impermeable, no-slip walls, representing the Be FWP volume, and an open region in which multiphase flow  is modelled. These two regions are separated by a thin wall with zero thermal resistance and a static contact angle of 10 • to emulate the good wetting of liquid Be on its own solid phase [36]. The geometry of the thin wall reproduces the 45 • chamfer at the FWP edge, located at x = y = 0. The positive x and y directions in Fluent correspond respectively to the negative Z and positive X directions in the MEMOS-U coordinates used in figure 2.
At the start of the simulation, the closed region is filled with solid Be, with a temperature profile as shown in figure 5, where z(x, y) is the distance to the initial plasma-facing (x 0, y = 0) surface: where the five main parameters are obtained by linear time interpolation from the data listed in table 1. Equation (9) corresponds to a Poiseuille flow in the melt layer, while (10) imposes a linear temperature profile in the metal, bounded above by T surf . The injected gaseous ghost fluid (for y > y m + h) has a uniform velocity matching that of the liquid surface and a uniform temperature T bulk . Plots of the inlet profiles are shown in figure 6 at t = 30 ms for easier visualisation. The total 179.25 mm 2 area of the simulation domain is meshed by approximately 160 000 cells, with a characteristic resolution of 25 μm in the region located within 2 mm of the thin wall, and 50 μm elsewhere. Navier-Stokes calculations are carried out using the PISO pressure-velocity coupling algorithm [37], with adaptive time stepping in order to maintain the Courant number of the flow below 0.1. This results in time steps of the order of 0.1 μs and a typical computational cost of 125 core-hours per millisecond of simulated time (using Intel R Xeon R Gold 6134 CPUs).

Evolution of melt instabilities
Due to the large values of the injected melt thickness and velocity, which translate into Weber numbers, We = ρhu 2 /γ, above 20 at the beginning of the simulation, free-surface instabilities start developing within a few milliseconds. In contrast to [21], where a straight PFC edge was shown to facilitate the complete jetting of an incoming melt layer at such We values, the milder chamfer angle allows the downstream flow to remain partially attached to the wall since the upstream momentum does not need to be transferred entirely onto the y direction as the liquid crosses the corner.
The increasingly irregular features of the flow during the first 15 ms are illustrated by figure 7. Panels (a)-(c) are representative of the initial geometry of the downstream layer, with a main bulge near the wall angle followed by smaller and shorter-wavelength free-surface oscillations whose amplitude grows beyond the linear regime, leading to wave breaking phenomena akin to the spilling (panel (b)) and plunging (panel (c)) breakers often observed in water waves [38][39][40]. The growing instabilities downstream of the main bulge lead to the formation of small sheet-like structures and first melt ejection events as shown in panels (d) and (e). Moreover, the strongly non-uniform layer thickness along the chamfer couples with heat conduction into the solid FWP, effectively splitting the liquid into nearly independent parcels, as was noted in [21]. This promotes further instabilities in the upstream flow until the main bulge breaks (panel (e)), which results in high-vorticity dynamics as exemplified by the formation of numerous gas-filled cavities in panels ( f )-(h).
The aftermath of this first instability growth stage is followed by a relative recovery of the molten layer as the main bulge reforms and the liquid thickness becomes more uniform, as shown by figure 8(a), which is qualitatively similar to figure 7(b). However, unstable downstream waves start developing again as demonstrated in figure 8(b). This global flow behaviour can be compared to the alternating bulge filling and relaxation dynamics observed in [21]. As the main corner bulge breaks a second time, it evolves into an elongated fluttering structure (figures 9(a)-(d)), whose morphology evokes the sinuous mode known in linear fluid sheet instability studies [41]. Since most of the injected Be melt becomes 'airborne', no thermal energy is supplied to the sloped FWP surface,  Under the action of surface tension forces, working to reduce the total area of the liquid's free surface, the aforementioned cavity ruptures into an elongated ligament. This event, displayed in figure 10, constitutes the last occurrence of material ejection in the simulation. Whereas Be injection into the computational domain continues for approximately 30 ms, the lower values of h, u and T surf at the inlet (see table 1) do not lead to additional free-surface break-ups. Instead of the main bulge reforming, figure 11(a) shows that a train of smallamplitude millimetric capillary waves starts propagating from the corner, sharing some qualitative features with the spilling dynamics observed experimentally in [39]. As the simulation progresses further, conduction cooling becomes dominant and the re-solidified Be layer on the chamfer grows to a thickness of nearly 1 mm ( figure 11(b)).

Material ejection
Approximately 10 melt ejection events are observed in the simulation output. Whereas a couple of the ligaments so created are eventually reintegrated into the main layer, most remain separated from it until they exit the computational domain. The equivalent radii of the ejecta range from 50 μm to 1.25 mm, leading to a total 9 mm 2 surface area of detached material, which can be compared to the 100 mm 2 of liquid Be injected at the inlet through the entire simulation. More than 80% of the The center-of-mass velocity of the ejecta typically lies between 3 m s −1 and 5 m s −1 , i.e. values comparable to those of the attached downstream flow. Although significant speed gradients-corresponding to variations of the order of 1 m s −1 across the ligaments' dimensions-can be observed, the velocity field is found to maintain a uniform direction, parallel to the sloped wall. Such internal shear configurations produce elastic oscillations on the free surface of the ejecta. These are apparently rapidly dissipated by viscosity in the case of small ligaments, which then adopt a circular cross-sectional shape, while larger ejecta are seen to retain an oscillating elliptic cross-section.

Discussion and conclusions
Given the chaotic nature of the simulated melt flow, the fine details of the observed dynamics are expected to depend heavily on initial and boundary conditions, and hence not to be readily extrapolated into predictions for ITER. Rather, the global qualitative features of the flow behaviour, as well as rough assessments of the main spatio-temporal scales involved, are likely to constitute robust results. A first remarkable feature of the numerical output is that analytical estimates from linear analyses of Rayleigh-Taylor instabilities are consistent with the observed flow dynamics. Indeed, the dispersion relation [42,43] ν 2 = γk ρh 2 We − k 2 h 2 tanh(kh), (11) where ν is the growth rate and k the wavenumber, predicts that the most unstable modes have a wavelength of a few millimetres and a characteristic growth time of the order of 0.5 ms, with larger-wavelength modes growing more slowly. This is reflected by the initial destabilization of the melt layer in figure 7, wherein short-wavelength downstream modes break before the larger main bulge. Hence, such simple linear estimates remain useful for order-of-magnitude predictions of melt flow behaviour. Despite the highly unstable dynamics of the melt flow, material ejection is found to represent only 10% of the incoming liquid. This is especially noteworthy in the light of previous numerical findings that Weber numbers as small as 5 can lead to complete jetting in right-angle corner configurations [4,21], and that two-dimensional simulations appear to overestimate ejecta sizes [21]. It should, however, be mentioned that the geometry of actual ITER FWPs includes a second 45 • angle, located roughly 4 cm away from the first one. The results presented in this paper suggest that some melt will reach this second corner in an already unstable state, likely leading to more ejection events. However, the velocities of both the attached melt flow and the detached ligaments are found to be parallel to the FWP chamfer, which implies that, at least for the particular FWP considered here, most of the ejected material moves towards remote, plasma-shadowed surfaces such as inter-panel gaps. There, liquid ejecta will likely splash and ultimately cool into frozen spatter traces similarly to what has been observed in JET [13]. Such splashing impacts may cause small secondary droplets to be released towards the main plasma chamber, but their mass is probably negligible with respect to the impinging liquid projectiles [44].
Overall, the results presented here can serve as a basis for a downward re-evaluation of the original upper-limit estimates of melt ejection during unmitigated VDEs [1], which were based on the conservative assumption that 25% of the entire melt produced by plasma-wall interaction is lost, and ultimately converted into dust. According to the MEMOS-U results [9] used in this paper, up to half of the liquid Be volume present on the surface of FWP #8 reaches the chamfer. Given the aforementioned simulation results showing that 10% of the liquid crossing the edge detaches from the wall, only 5% of the total induced melt is expected to be ejected in the worstcase scenario. Moreover, if such lost material does deposit on plasma-shadowed areas, its conversion into solid dust is unlikely [5]. Nevertheless, global predictive assessments of melt creation and destabilization during ITER disruptions will not be complete for as long as the TQ stage remains an open issue. This will be the subject of future work. IO/20/CT/4300002230). ITER is the Nuclear Facility INB No. 174. This publication is provided for scientific purposes only. Its contents should not be considered as commitments from the ITER Organization as a nuclear operator in the frame of the licensing process. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.