Simulations of liquid metal flows over plasma-facing component edges and application to beryllium melt events in JET

Navier–Stokes simulations of liquid beryllium (Be) flows over the straight edge of plasma-facing components are carried out in conditions emulating upper dump plate (UDP) melting observed experimentally in JET. The results demonstrate the existence of three main hydrodynamic regimes featuring various degrees of downstream flow attachment to the underlying solid surface. Transitions between these regimes are characterized by critical values of the Weber number, which quantifies the relative strength of fluid inertia and surface tension, thereby providing a general stability criterion that can be applied to any instance of transient melt events in fusion devices. The predictive capabilities of the model are tested by comparing numerical output with JET data regarding the morphology of the frozen melt layers and the location of Be droplets splashed onto nearby vacuum vessel surfaces as a result of disruption current quench plasmas interacting with the solid Be tiles protecting the upper main chamber regions. Simulations accounting for the coupling between fluid flow and heat transfer confirm the key role played by re-solidification as a stabilizing process, as previously found through macroscopic melt dynamics calculations performed with the MEMOS-U code. The favourable agreement found between the simulations and the general characteristics of the JET Be UDP melt splashing give confidence that the same approach can be applied to estimate the possibility of such mechanisms occurring during disruptions on ITER.


Introduction
Transient surface melting of plasma-facing components (PFCs) under high-energy plasma loads has been observed * 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.
repeatedly in tokamaks equipped with a metallic first wall [1][2][3][4][5]. The dynamics of the resulting liquid pools are potentially subject to a variety of instabilities that lead to material ejection in the form of droplets, which may then re-solidify into dust particles. In fact, transient melt events are foreseen as a major contributor to the beryllium (Be) dust production in ITER, along with the delamination of co-deposited impurity layers [6]. Reliable characterizations of melt splashing are therefore crucial to ensure that the safety limits imposed on the in-vessel dust inventory [6,7] are met. This issue has received renewed attention following recent studies of PFC thermal response and macroscopic surface melt dynamics in which the MEMOS-U code [8,9] is applied to the current quench phase of unmitigated vertical displacement events in ITER, predicting the formation of extended melt pools on the Be first wall panels [10,11].
JET operation with the ITER-like wall offers clear empirical evidence that liquid surface layers are produced on Be PFCs [2,4]. In particular, extensive melting is observed on the upper dump plates (UDPs) as a consequence of upward vertical displacement events [4]. Under the action of Lorentz forces due to halo currents, Be melt produced on the plasmawetted face of the UDPs is found to move over the plate edge towards the low-field side. While crossing the edge, a fraction of the liquid Be ejects as droplets onto the nearby vacuum vessel wall, whereas the remaining melt flows on the plasmashadowed face of the UDP and ultimately solidifies into a so-called 'waterfall' structure [4]. Splashing from PFC edges is in fact the most prominent example of macroscopic melt instability in existing tokamaks [3]. It involves physical mechanisms which are discarded in the shallow-water equations implemented in MEMOS-U and must therefore be the subject of separate dedicated studies. In ITER, the geometrical constraints for such processes to develop may be satisfied at the edge of first wall panels or in case of radial misalignment of poloidal PFC stacks.
This paper presents first modelling results on the stability of liquid Be flows over a right-angle corner geometry similar to that of the JET UDPs. Two-dimensional incompressible Navier-Stokes simulations are carried out using customized ANSYS Fluent set-ups that were benchmarked against previous publications on pressure-driven flows during unipolar arcing [12]. Repeated calculations are run to explore the full range of melt depths and velocities relevant to JET current quenches [8]. After a discussion of the various predicted hydrodynamic regimes and their stability, momentum damping due to gradual re-solidification on the plasma-shadowed area [8] is investigated by means of simulations accounting for the coupling between fluid dynamics and heat transfer. The results are then compared to the available experimental data regarding waterfall morphology and the presence of frozen Be splash traces on the JET vacuum vessel [4].

Hydrodynamic simulations
ANSYS Fluent simulations are set up to solve the twodimensional incompressible Navier-Stokes system using the volume-of-fluid method [13] to track the position of the free surface which separates the liquid metal from the background gaseous phase (which corresponds physically to dilute gas or a weakly-ionized plasma remnant). The model equations are written in the so-called one-fluid formulation, which means in particular that the pressure and velocity fields are shared by the liquid metal and the gas: ∂α ∂t where v is the velocity, α is the metal volume fraction, p is the pressure, ρ is the mass density and μ is the dynamic viscosity. The volumetric body force term F includes gravity, the Lorentz force due to disruption halo currents, and surface tension effects via the continuum surface force model [14]: where g denotes the gravitational acceleration, γ the surface tension and κ the algebraic curvature of the free surface. Note that gravity, which is oriented in the positive y direction (see figure 1) so as to mimic the experimental conditions at the UDPs in JET, was found to have no noticeable effect on the simulation results and could in fact have been neglected. Local material properties are obtained as volume-fraction averages, e.g.
in which the subscripts 'm' and 'g' respectively refer to the liquid metal and the gaseous phase. As in [12], the properties assigned to the gaseous phase are not meant to be physically accurate; they are chosen to be negligible with respect to their metal counterparts in order to ensure that the gas behaves as a so-called ghost fluid which does not affect the dynamics of the liquid. Liquid Be is modelled with ρ m = 1690 kg m −3 , μ m = 1.2 × 10 −3 Pa s −1 and γ = 1.3 N m −1 [15][16][17].
The simulations are run on an inverted L-shaped domain inscribed in an 8 × 10mm 2 rectangle, sketched in figure 1 as the 'open region'. It is meshed with approximately 75 000 quad cells whose characteristic size ranges from 25 μm in the part of the domain defined by x 3 mm and y 1 mm, where most of the liquid is located, to 50 μm in the remaining area. Moreover, 10 layers of thin mesh cells are stacked in the region within 35 μm of the no-slip walls in order to resolve velocity gradients at these locations. The equilibrium contact angle at the walls is set to 10 • in order to emulate the good wettability of solid metals by their own liquid phase [18]. Starting from an initial state in which the domain is entirely filled with immobile gas, a liquid metal layer of height h is injected at the inlet with a fully developed half-parabolic velocity profile whose average is v. Above the liquid, the gas injection velocity is uniform so as to enforce zero shear stress at the free surface, Once in the domain, the liquid layer is subjected to a horizontal Lorentz force whose value compensates viscous drag in order to propagate the inlet velocity profile up to the point where the melt reaches the corner: The simulations are run in laminar mode with the PISO pressure-velocity coupling algorithm [19]. Adaptive time stepping is controlled by a set-point Courant number of 0.1, leading to steps in the range 0.1-0.5 μs. Each time step consists of 20 internal solver iterations, which yields a representative computational cost of 30 core-hours per simulated millisecond of physical time (using Intel R Xeon R Gold 6134 CPUs).

Simulations including heat transfer
In order to test the impact of heat transfer, a second simulation set-up including thermal modelling has been prepared. With respect to the hydrodynamic simulations described in section 2.1, the heat convection-diffusion equation is added to the physical model, where ε is the specific enthalpy, T is the temperature and λ is the thermal conductivity. Metal solidification and melting are implemented as detailed in [12] through the introduction of the molten fraction β (T), which linearly increases from 0 to 1 in the range 1559-1561 K, and the addition of a springlike body force term to the momentum equation (3), which enforces zero velocity when β = 0. The thermal conductivity and specific heat capacity of Be are assumed to be constant at λ m = 80 W m −1 K −1 and c m = 3200 J kg −1 K −1 , respectively, while L melt = 1.64 × 10 6 J kg −1 is used for the latent heat of fusion [20,21]. The density and thermal properties of solid Be are set to be equal to their liquid-phase counterparts, hence thermal expansion effects, along with all other phenomena that can attributed to the temperature dependence of material properties, are neglected in the simulations. Therefore, differences with hydrodynamic simulations can be imputed exclusively to heat diffusion and phase transitions. Although non-uniformities associated to the temperature dependence of material properties can significantly affect the detailed features of the flow (such as the shape of free-surface oscillations), they are not expected to alter the main phenomena and characteristic estimates discussed in sections 3 and 4. The simulation domain is extended to the entire 8 × 10 mm 2 rectangle, including the 'wall region' sketched in figure 1, which is meshed with square cells of size 25 μm and is entirely filled with metal. All four boundaries of the wall region are impermeable no-slip walls, implying that neither mass nor momentum can be transferred between the two regions. The two walls in contact with the open region are infinitely thin and have zero thermal resistance, meaning that heat transfer between the two regions is unimpeded. Finally, the adiabatic boundary condition is imposed on the remaining two walls.
At the start of the simulation, the temperature in the wall region is assumed to follow the analytical solution for onedimensional heat diffusion under constant heat flux: where T bulk and T surf are the bulk and surface temperatures of the solid, respectively, and δ is the characteristic heat diffusion length. Here, T surf = 1500 K is chosen to be slightly below the melting point T melt = 1560 K of Be, while T bulk = 580 K and δ = λ m √ π (T surf − T bulk ) /q 620 μm have been selected in accordance with the initial plate temperature and average plasma heat flux q = 210 MW m −2 used in MEMOS-U simulations of the JET current quench induced Be melting [8]. The initial temperature of the gas is set at T bulk in order to avoid spurious temperature gradients that would result in artificial heat fluxes into the wall region. The temperature of the liquid metal injected at the inlet is T in = 1700 K. At such temperatures, the radiative cooling flux at the free surface remains below 0.5 MW m −2 even if perfect thermal emissivity is assumed. This is much smaller than the characteristic conduction heat flux λ m (T in − T melt ) /h > 50 MW m −2 , which justifies neglecting thermal radiation in equation (8). Due to the extra equation to solve and the larger computational domain, the typical cost of simulations with heat transfer is 45 core-hours per simulated millisecond of physical time.

Hydrodynamic flow regimes
A total of 26   dynamics modelling [8]. Calculations were run for a maximum physical time of 50 ms and the results can be broadly classified in three main categories: the attached regime, the quasi-steady jetting regime, and the intermediate regime.
The attached regime is realized for relatively small values of h and v and corresponds to simulations in which the flow remains completely attached to the wall during the entirety of the 50 ms of simulation time, never entering in contact with the right outlet boundary. The typical dynamics of liquid layers in the attached regime are shown in figure 2: in the vicinity of the corner, the competition between inertia and surface tension leads to the formation of a bulge, which gets filled by the incoming metal until the curvature of the free surface becomes small enough to allow the fluid's velocity vector to rotate by 90 • . Fluid parcels whose velocity completes the rotation are transferred along the vertical wall, which translates as a relaxation of the bulge into a free-surface wave propagating in the negative y direction. As new liquid is injected at the inlet, the process repeats and the downstream free surface oscillates globally as the content of the bulge transfers vertically in a pulse-like manner. As figure 2 demonstrates, the amplitude of the bulge oscillations can reach several millimeters, at which point downstream mass transfers occur in the form of a breaking wave, which may lead to the formation of gas-filled cavities inside the melt layer.
In the quasi-steady jetting regime, which corresponds to large values of h and v, the inertia of the liquid is too large to allow the fluid to accelerate into bulk vertical motion. As a consequence, the melt layer jets away from the corner as a single sheet, whose angle with respect to the horizontal wall echoes the incomplete rotation of the fluid velocity vector near the corner. After some time, the flow eventually enters a quasi-steady state in which the free surface remains essentially unchanged until the end of the simulation. Two examples of these quasi-steady states, with different jetting angles, are shown in figure 3.
Simulation parameters leading to results fitting in neither of the two aforementioned categories are classified as belonging to the intermediate regime. In this regime, free surface oscillations near the corner are large enough for the metal to interact with the right computational boundary, but the flow does not settle into a quasi-steady state within 50 ms. Contact with the right outlet boundary can occur via the ejection of large cylindrical ligaments from the corner bulge (see figure 4) or via the formation of elongated liquid sheets which eventually collapse onto the vertical wall (see figure 5). It can be noted that the time required for the oscillation amplitude to become large enough to induce material ejection or interaction with the right boundary is of the order of 20 ms, that is comparable with the characteristic current quench timescale in JET [8].

Classification by Weber numbers
Additional insight into the hydrodynamic regimes described in section 3.1 can be gained by exploiting a crude analogy between the flows of interest and the development of Rayleigh-Taylor instabilities at the surface of a planar fluid film. This analogy consists in considering the liquid layer as if it were lying over a flat surface and subjected to a normal body force F ∼ ρ m v 2 /h, whose magnitude emulates the centrifugal force that appears in a reference frame following the rotation of the actual fluid across the corner [22][23][24]. Under this assumption, the dispersion relation of linearly unstable Rayleigh-Taylor modes for an inviscid film of thickness h [25], where k is the wavenumber, ν the growth rate and Bo = Fh 2 /γ the bond number, can be applied, provided that Bo is replaced by the Weber number We = ρ m hv 2 /γ. Given that the characteristic time required by the melt to cross the corner is h/v, it is expected that an unstable mode can develop when hν/v is larger than some critical dimensionless value. Plugging in the largest growth rate predicted by equation (10)  Finally, it should be emphasized that dimensionless stability criteria can be directly applied to any metal. For example, transient surface melting accompanied by melt ejection has been observed during the exposure of a tungsten leading edge in ASDEX Upgrade [3]. That experiment is similar to the JET UDP case insofar as ejection occurs when liquid layers reach the straight edge of the sample. Detailed MEMOS-U simulations of the ASDEX Upgrade exposure predict that melt pools arrive at the corner at a velocity of 1-2 m s −1 with a typical thickness of 300 μm. Using ρ m = 1.6 × 10 4 kg m −3 and γ = 2.4 N m −1 as representative properties for liquid tungsten [26], the corresponding Weber numbers are in the range 2-8. These values fall mostly in the quasi-steady jetting regime, in agreement with the observations of re-solidified melt and debris documented in figure 9(d) of [3].

Effect of re-solidification
Owing to the increased computational cost of simulations coupling fluid dynamics and thermal modelling, the thorough parameter scan presented above has not been repeated. Instead, the particular scenario shown in figure 4, with h = 200 μm and v = 3 m s −1 is selected for simulations including heat transfer. While the flow belongs to the intermediate regime of hydrodynamic simulations, once heat transfer is taken into account, it falls into the more stable attached regime. This stabilization is due to the progressive re-solidification of the melt layer downstream of the corner, where the relatively cold solid surface conducts thermal energy away from the liquid. This effect, which has already been discussed in MEMOS-U shallow-water simulations [8], is displayed in figure 7, where it can be seen that the progression of the solidification front reduces the downstream melt depth, leading in particular to higher viscous drag.

Melt ejection angles and morphology of the re-solidified surface
Post-mortem measurements of splashed and re-solidified Be melt near the JET UDPs are reported in [4]. The published data consist mainly of images of the solid waterfall-like structure formed by melt that crossed the plate's edge onto the plasmashadowed side, as well as splash traces found on the nearby vacuum vessel wall. As such, most of the available information is qualitative, although some characteristic length scales can be quantified. In particular, the range of droplet ejection angles can be estimated from photographs as 0 • -40 • by using a nearby 'mushroom' limiter as a landmark (see figure 8(a)). These angles are consistent with those found in some simulation results in the intermediate regime when liquid ligaments are ejected, but also with low-We instances of the quasi-steady jetting regime in which the liquid sheet exits the simulation domain through the lower boundary (as in the left panel of figure 3). The presence of millimeter-range, roughly spherical features on the waterfall (see figure 8(b)) is also consistent with the liquid bulges and large free-surface oscillations observed in the simulations.
Furthermore, numerical results provide some insight on the complex multi-layer structure of the frozen Be waterfall (see figure 8(c)). Indeed, while the most natural interpretation of the experimental images is to link each layer with a distinct melt event, simulations including heat transfer open the possibility for some of the layers to originate from a single disruption. The proposed mechanism stems from the intermittent melt mass transfers from the corner to the vertical side of the plate and is illustrated in figure 7. Whereas the liquid filling the corner bulge remains well above the melting point due to the relatively large distances separating it from the cold wall, the morphology of the thinner melt further away from the corner enables much more efficient conduction cooling, possibly allowing for complete downstream re-solidification before the next mass transfer occurs. The hot liquid then transported from the corner effectively forms a new layer that propagates onto the recently solidified material. These mass transfers can even be accompanied by local re-melting of the underlying layer, contributing further to the non-uniformity of the downstream flow. Overall, although surface composition analyses presented in [4] indicate that the waterfall was not formed in a single disruption, simulations suggest that it is impossible to establish a oneto-one correspondence between the number of layers and the number of melt events.

Characteristic size and mass of the ejecta
Since the information available regarding the splash traces found on the vacuum vessel wall consists exclusively of photographs, cross-sectional area is the only quantity that can reliably be measured. Extra assumptions on the real threedimensional spatter shape are required in order to estimate the ejected mass and the size of the ejecta before their collision with the wall. Two extreme hypotheses are used in [4]: the measurable surface is interpreted as either the cross-section of a sphere or the top side of a flattened cylinder with height h sp = 0.8 μm. Taking into account a resolution of 6 px mm −1 [4], which corresponds to a smallest measurable surface S px = 2.78 × 10 −8 m 2 , the radius of the smallest spherical droplet whose trace can be detected on the images is estimated as either S px /π = 94 μm if spherical shape is assumed, or 3S px h sp /4π 1/3 = 17 μm for a cylindrical shape. Likewise, if one considers that the largest splash trace visible on the photographs is made of roughly 50 px, the corresponding droplet radius is either 665 μm or 64 μm. This large uncertainty is reflected in the two limiting estimates reported in [4] for the total mass splashed on the wall, which differ by a factor 660. The spherical shape assumption is ultimately the one retained by Jepu et al based on the similarity between the resulting splashed mass and their estimate of the total UDP mass loss. However, this conclusion is questionable because: (i) weight loss measurements are carried out on tiles located roughly half a meter away poloidally from the UDP edge where melt ejection occurs, in particular no weighing data is presented for the tile where the frozen waterfall is observed; (ii) the measured mass loss varies greatly between different tiles, to the point that the total mass loss estimate can be almost entirely accounted for by the contribution from a single tile (see figure 5 in [4]). Physical models of droplet-wall impacts can be exploited in an attempt to reduce the uncertainty of the splash mass estimates. In such models, the outcome of the collision between a hot liquid drop and a cold solid surface (i.e. well below the melting point) is determined by the competition between spreading dynamics driven by inertia and wetting, and damping mechanisms governed by viscosity and re-solidification [27]. The relative importance of each of these processes depends on the size and velocity of the impinging droplet, as well as the respective temperatures and thermo-physical properties of both materials in contact. Analytical estimates of the maximum droplet spread are available in the literature, based on considerations of the global energy budget of the impact. Namely, the initial energy of the droplet-comprised of kinetic and surface contributions-is dissipated into viscous damping, motion suppression at the re-solidification front, and the formation of new interfaces with the underlying solid and the background gaseous medium [28]. In the context of interest here, typical values of the maximum spread factor ξ max , i.e. the ratio of the re-solidified spatter radius to the initial droplet radius, are estimated to lie in the range ξ max = 3-5 [28]. This suggests that the assumption of cylindrical splash  shape, which corresponds to ξ max = 16S px /9πh 2 sp 1/6 = 5.4 for a single-pixel trace, is closer to reality. It should however be emphasized that these estimates are only indicative; as noted in [28], more refined models accounting for thermal contact resistance and the three-dimensional nature of heat conduction into the substrate may be required for a precise comparison.
Nevertheless, the ligaments ejected during simulations belonging to the intermediate regime, with equivalent radii in the range 600 μm-1.2 mm, are found to be significantly larger than experimental estimates. Several potential explanations may be brought forward to explain this discrepancy. These include: (i) the lack of a third spatial dimension in the simulations, which can modify the way surface tension acts to stabilize the flow and determine the shape and size of the ejecta; (ii) the effect of realistic time-varying melt thickness and velocity at the inlet, which influences free-surface dynamics near the corner; (iii) possible break-up phenomena fed by the kinetic energy stored in surface oscillations of the ligaments. It may also be noted that the fluid sheets present in the quasi-steady jetting regime are likely subject to Plateau-Rayleigh instabilities [29], causing them to split into droplets whose characteristic size is comparable to that of the original sheet, i.e. the incoming melt depth. Such sizes would then be in much closer agreement with the experimental data. A detailed study of the relative importance of each of the aforementioned phenomena requires three-dimensional simulations run on extended spatial domains. The computational cost of such calculations effectively forbids parametric surveys; more realistic simulations of concrete, fully specified transient melt events will be the object of future work.

Conclusions and outlook
Two-dimensional ANSYS Fluent Navier-Stokes simulations of liquid Be layers flowing over a right-angle corner were run in conditions emulating disruption current quench induced melting of Be UDP tiles of the JET ITER-like wall [4]. A parametric survey of the relevant ranges of melt depths and velocities [8] was conducted, demonstrating the possibility for three main flow regimes that differ by the degree of attachment of the liquid on the downstream surface. Transitions between these regimes are shown to be governed by critical values of the fluid Weber number, thereby providing a dimensionless stability criterion that can be applied to any material. Detailed comparisons between simulation results and surface analyses carried out at JET [4] display good agreement in terms of melt ejection angles and morphological features of the re-solidified metal. However, significant discrepancies with experimental data are found in regard to the characteristic ejecta sizes. These are presumably due in the most part to the absence of three-dimensional processes in the simulations, the use of time-independent inlet parameters and the possibility of liquid break-up phenomena occurring on larger spatio-temporal scales than those sampled in the calculations.
First simulation attempts coupling hydrodynamics and heat transfer confirm the stabilizing effect of melt re-solidification on plasma-shadowed surfaces, as previously found by macroscopic shallow-water calculations using the MEMOS-U code [8]. Juxtapositions of numerical results with and without thermal modelling elucidate the effective reduction of the melt thickness downstream of the corner and its role in suppressing free-surface oscillations.
The general understanding of melt flows near PFC edges gained from parametric surveys will serve as a basis for the analysis of future numerical studies focused on concrete current quench scenarios in ITER, to determine whether significant dust and droplet production can occur during melt events. Three-dimensional simulations are expected to yield better quantitative estimates of the characteristics of melt ejection processes and allow the investigation of transverse phenomena such as the possible formation of rivulets [30]. The spreading dynamics of hot droplets impinging on cold surfaces in fusionrelevant conditions will also be explored separately using dedicated calculations. This will not only facilitate the link between numerical results and post-mortem observations of splashing traces, but also open the way for more accurate descriptions of mechanical impacts in dust and droplet transport codes [31].