Beryllium melt instabilities and ejection during unmitigated current quenches in ITER

Quantitative estimates of metallic dust generation during the operation of fusion reactors play a pivotal role in nuclear safety and licensing issues [1–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 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–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–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–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^® 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 break-ups and material ejection events, whose characteristics are of great importance to dust migration, survival and accumulation studies [22].

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 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.

Zoom In Zoom Out Reset image size

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. The spatially varying properties of the entire melt layer at t = 0 are shown in more details in figure 2. The point $\left(Y=50\enspace \text{cm},Z=0.5\enspace \text{cm}\right)$ in the MEMOS-U domain is chosen as the reference location from which the melt characteristics are extracted. More precisely, the instantaneous state of the molten layer is summarized by five main figures of merit:

• The position y_m of the solidification front with respect to the nominal FWP surface,
• The melt thickness h,
• The depth-averaged velocity u,
• The surface temperature T_surf,
• The temperature gradient ∇T_m at the solidification front, in the direction normal to the nominal FWP surface.

Zoom In Zoom Out Reset image size

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.

Table 1. Melt layer characteristics extracted from MEMOS-U simulations and used in equations (8)–(10) to construct inlet boundary conditions.

t	$\left[\text{ms}\right]$	0	20	40	60
ym	$\left[\mu \text{m}\right]$	0	0	124	319
h	$\left[\mu \text{m}\right]$	531	531	265	70
u	$\left[\text{m}\enspace {\text{s}}^{-1}\right]$	5.36	4.58	2.87	0
Tsurf	$\left[\text{K}\right]$	1795	1706	1566	1566
∇Tm	$\left[\text{K}\enspace \mu {\text{m}}^{-1}\right]$	1.403	1.375	0.730	0.500

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², 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³, 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⁻³ ≪ Ca^1/3. Therefore, the anomalous melt thickness datum at t = 0 is discarded and $h\left(t=0\right)$ is extrapolated backwards from $h\left(t=20\enspace \text{ms}\right)$, assuming by default that it remains constant during that interval.

3.1. Main equations and material properties

Similarly to [20, 21], Fluent is used to solve the two-dimensional Navier–Stokes system for an incompressible two-phase flow with heat transfer using the volume-of-fluid method [27]:

$$\begin{equation}\nabla \cdot \boldsymbol{v}=0,\qquad \qquad\end{equation} \tag{ 1 }$$

$$\begin{equation}\frac{\partial \alpha }{\partial t}+\nabla \cdot \left(\alpha \boldsymbol{v}\right)=0,\qquad \qquad\qquad \qquad\end{equation} \tag{ 2 }$$

$$\begin{align}\hfill \frac{\partial }{\partial t}\left(\rho \boldsymbol{v}\right)+\nabla \cdot \left(\rho \boldsymbol{v}\boldsymbol{v}\right)& =\nabla \cdot \left[\mu \left(\nabla \boldsymbol{v}+\nabla {\boldsymbol{v}}^{\top }\right)\right]\hfill \\ \hfill & \quad -\nabla p+\rho \boldsymbol{g}+\gamma \kappa \nabla \alpha +{\boldsymbol{F}}_{\text{s}},\hfill \end{align} \tag{ 3 }$$

$$\begin{equation}\frac{\partial }{\partial t}\left(\rho \varepsilon \right)+\nabla \cdot \left(\rho \varepsilon \boldsymbol{v}\right)=\nabla \cdot \left(\lambda \nabla T\right),\qquad \qquad\quad \end{equation} \tag{ 4 }$$

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

$$\begin{equation}\rho \varepsilon ={\int }^{T}{c}_{\text{v}}\left({T}^{\prime }\right)\mathrm{d}{T}^{\prime }+\rho \beta \left(T\right){\Delta}{H}_{\text{fus}}\end{equation} \tag{ 5 }$$

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 β, which is implemented as a piecewise linear approximation of $\text{H}\left(T-{T}_{\text{melt}}\right)$ [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 ρ = 1690 kg m⁻³, T_melt = 1560 K and ΔH_fus = 1.64 MJ kg⁻¹ [31–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.

Zoom In Zoom Out Reset image size

3.2. 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.

Zoom In Zoom Out Reset image size

At the start of the simulation, the closed region is filled with solid Be, with a temperature profile as shown in figure 5,

$$\begin{align}\hfill {T}_{\text{wall}}^{0}\left(z\right)& ={T}_{\text{bulk}}+\left({T}_{\text{melt}}-{T}_{\text{bulk}}\right)\hfill \\ \hfill & \quad \times \left[\mathrm{exp}\left(-\frac{{z}^{2}}{{\delta }^{2}}\right)-\sqrt{\pi }\frac{z}{\delta }\mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{z}{\delta }\right)\right],\hfill \end{align} \tag{ 6 }$$

where $z\left(x,y\right)$ is the distance to the initial plasma-facing $\left(x\leqslant 0,y=0\right)$ surface:

$$\begin{equation}z\left(x,y\right)=\begin{cases}\vert y\vert \quad \hfill & \quad \text{if}\;x\leqslant 0,\hfill \\ \sqrt{{x}^{2}+{y}^{2}}\quad \hfill & \quad \text{if}\;x > 0.\hfill \end{cases}\end{equation} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

Equation (6) corresponds to one-dimensional heat diffusion in a uniform material under constant surface heat flux. The bulk FWP temperature T_bulk = 373 K and heat diffusion length δ = 1.5 mm are chosen to approximate MEMOS-U output at t = 0. The open region is initially filled entirely by immobile gas at T_bulk to minimize spurious heat transfer across the thin wall during early simulation time steps, before it is covered by liquid.

The outer computational boundaries of the closed region are adiabatic walls, whereas the open region is externally bounded by constant-pressure outlets and a velocity inlet (x = −2 mm, y ⩾ 0) through which gas and liquid Be are injected. Inlet profiles for the volume fraction, velocity along x and temperature are defined by

$$\begin{equation}{\alpha }_{\text{in}}=\text{H}\left({y}_{\text{m}}+h-y\right),\end{equation} \tag{ 8 }$$

$$\begin{equation}{v}_{\text{in}}=\frac{3u}{2{h}^{2}}\text{H}\left(y-{y}_{\text{m}}\right)\left[{h}^{2}-{\left({y}_{\text{m}}+h-y\right)}^{2}{\alpha }_{\text{in}}\right],\end{equation} \tag{ 9 }$$

$$\begin{align}\hfill {T}_{\text{in}}& ={\alpha }_{\text{in}}\mathrm{min}\left[{T}_{\text{surf}},{T}_{\text{melt}}+\left(y-{y}_{\text{m}}\right)\nabla {T}_{\text{m}}\right]\hfill \\ \hfill & \quad +\left(1-{\alpha }_{\text{in}}\right){T}_{\text{bulk}},\hfill \end{align} \tag{ 10 }$$

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.

Zoom In Zoom Out Reset image size

The total 179.25 mm² 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^® Xeon^® Gold 6134 CPUs).

4.1. Evolution of melt instabilities

Due to the large values of the injected melt thickness and velocity, which translate into Weber numbers, We = ρhu²/γ, 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–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).

Zoom In Zoom Out Reset image size

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, allowing the solidification front to progress. Ultimately, the sheet collapses by ejecting several fluid ligaments and forming a large cavity as detailed in figures 9(e)–(h).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 small-amplitude 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)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.2. Material ejection

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]

$$\begin{equation}{\nu }^{2}=\frac{\gamma k}{\rho {h}^{2}}\left(\text{We}-{k}^{2}{h}^{2}\right)\mathrm{tanh}\left(kh\right),\end{equation} \tag{ 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 worst-case 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.

This work has been carried out as part of a comprehensive study of disruption transient-induced beryllium melt splashing funded by the ITER Organization (Service Contract 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.