Resolidification-controlled melt dynamics under fast transient tokamak plasma loads

Studies of macroscopic melt motion induced by fast transient power loads and the ensuing damage to plasma-facing components are critical for ITER. Based on state-of-the-art experiments, heat transfer is argued to be strongly entangled with fluid dynamics. The physics model of the MEMOS-U code is introduced. Simulations are reported of recent tokamak experiments concerning deliberate transient melting of beryllium main chamber tiles (JET) and tungsten divertor components (ASDEX Upgrade, JET). MEMOS-U is demonstrated to capture the main physics responsible for melt dynamics and to reproduce the observed surface deformation. The intricate role of resolidification is elucidated.

be overcome in the development of fusion reactors. ITER will feature a beryllium (Be) first wall and a full tungsten (W) divertor [1], whose design is mainly driven by the steady-state power handling requirement [2]. However, the PFC integrity is mostly threatened by fast transient power loading on millisecond timescales, most commonly owing to magnetohydrodynamic instabilities such as edge-localized modes (ELMs), vertical displacement events (VDEs) or major disruptions. In fact, key aims of the Pre-Fusion Power Operation phase of ITER include the demonstration of effective disruption mitigation and reliable ELM suppression [3]. Due to the high plasma energies and power flux densities, even when mitigation schemes are considered, surface melting is essentially inevitable in the non-active phases of ITER operation [2,4]. Metallic melt is subject to external forces which displace the material and may cause large-scale surface deformation, severely compromising power handling [5]. Hence, reliable assessments of the possible consequences of transient melt events are crucial for ITER and beyond.
Self-consistent simulations of such problems are highly challenging, since simultaneous treatment of the severely deformed free surface and moving solid-liquid boundary is complicated by the intrinsic 3D geometry and the prolonged plasma exposures, ranging from several tens of milliseconds in disruption scenarios to several seconds in repetitive transient loading. Simplifications have been sought that are based on timescale separation assumptions between fluid dynamics and heat transfer. For heat diffusion much slower than fluid instability growth and liquid splashing, thermal analysis can be discarded and the problem reduces to the evolution of a prescribed melt pool [6][7][8]. For heat diffusion much faster than macroscopic fluid motion, thermal modelling can be performed in a static geometry, circumventing the need to introduce deformed domains [9,10]. Both approaches are challenged by a careful examination of recent experimental evidence which reveals that the timescales of heat transfer and fluid motion are comparable. This refers in particular to dedicated experiments on repetitive transient melting in JET [11] and ASDEX-Upgrade (AUG) [12] with deliberately misaligned W components, along with observations of disruption-induced Be melting in JET [13].
During interaction with repetitive ELMs or disrupting plasmas, the wetted area of PFCs is typically much less than their overall area [11][12][13]. Hence, newly formed liquid metal pools are surrounded by a progressively colder solid surface. Under the effect of plasma-induced forces, molten material is accelerated out of the pool and onto the solid, where it is subject to strong temperature gradients that facilitate heat conduction and prompt resolidification. We reveal that in such conditions, melt motion falls under an as-yet unexplored regime in which resolidification governs dynamics by determining the melt layer thickness and thus controlling viscous damping. The proposed new model, implemented in the MEMOS-U code, a significantly upgraded version of MEMOS-3D [9][10][11], is tested against tokamak experiments that involve different plasma scenarios, as well as PFC materials and exposure geometries. The successful reproduction of key experimental results demonstrates that the model adequately accounts for the essential physical processes at play.

Model equations of MEMOS-U
The macroscopic motion of liquid metals is described selfconsistently by the incompressible resistive thermoelectric magnetohydrodynamic equations and the heat convectiondiffusion equation [14]. Post-mortem observations on PFCs subject to tokamak-induced melting [11][12][13] have revealed that the melt extent greatly exceeds its depth and that a predominant direction exists in the frozen flow without signatures of liquid disintegration, justifying the shallow water approximation. The shallow water equations are obtained by integrating the Navier-Stokes equations across the liquid depth and reduce the fluid equation dimensionality by one [15], drastically decreasing computational cost. Combined with the magnetostatic limit [16] in the case of uniform material composition, the equations become Here (h, b 1 ,ẋ vap , U, J, B, T, T s , P, f d ) are the melt layer thickness, liquid-solid interface position, free surface position rate of change due to vaporization, depth-averaged fluid velocity, current density, magnetic flux density, bulk and surface temperature, ambient pressure, and plasma drag force; (ρ m , c p , σ, k, µ, ρ e , S) denote the temperature-dependent mass density, specific isobaric heat capacity, surface tension, thermal conductivity, dynamic viscosity, electrical resistivity, and absolute thermoelectric power; ∇ t refers to the gradient along the sample surface and ⟨⟩ t denotes depth-averaging. Note that the plasma pressure and plasma drag force constitute external input. MEMOS-U employs the finite difference method on a rectilinear grid. The fluid equations, equations (1), (2), are solved on a 2D geometry along the sample surface imposing a free flow boundary condition at the edge of the melt pool. The heat equation, equation (3), is solved on the full 3D domain that corresponds to the actual sample with surface fluxes due to plasma heating, vaporization, thermal radiation and thermionic emission [17,18]. Latent heat transfer is evaluated with the heat integration method, wherein the phase change enthalpy is counted in each grid point [19]. Evaporation mass loss follows the Knudsen formula [20]. The current density distribution should be solved on the 3D domain imposing the surface current. In ELM-driven W melt motion applications, this is the local escaping thermionic current dictated by the surface temperature and plasma fluxes [16]. In disruption-driven Be melt motion scenarios, it is strongly plasma-dependent and estimated from measurements with simplifying assumptions.

Transient W melting
Very similar experiments have been performed in JET [11] and AUG [12] by deliberately inserting a 1 mm W leading edge of short poloidal extent in the outer strike point region such that the full parallel-field heat flux is intercepted. Exposures to repetitive, unmitigated Type-I ELMs in high power H-modes provided quasi-periodic transient heat flux bursts superimposed on the much weaker inter-ELM stationary loads leading to repeated cycles of melting, motion and resolidification. Near the W melting point (3695 K), thermionic emission becomes important. Assuming plasma ambipolarity, the escaping thermionic current J esc th triggers a replacement current (with a ⟨J r ⟩ ∥ ≃ J esc th /4 [16] depth-averaged radial component) that flows through the molten layer and drives the dominant J × B force density. Due to space-charge limitation [21,22], J esc th is suppressed to a constant plasma-dependent value varying from 15 MAm −2 near the ELM peak to 5 MAm −2 in inter-ELM periods [23,24].
Furthermore, thermionic emission introduces a complex dependence on the plasma conditions and an intricate coupling between momentum and energy exchange. More specifically, the resulting current density dictates the melt acceleration and simultaneously affects surface cooling via the flux (J esc th /e)W f , with W f the work function. This crucial role has been recognized in MEMOS-U, where the description of the escaping thermionic current is based on dedicated particle-in-cell simulations [24,25].
To elucidate the interplay between fluid dynamics and heat transfer, we first consider a single ELM heat pulse and then the combined effect of consecutive ELM pulses. It will be instructive to discuss the ordering between the fluid motion and the heat transfer timescales. Since the liquid pools under question are strongly non-uniform and transient, we shall pay attention to the characteristic motion timescale responsible for the final surface deformation profile. This can be obtained from the pool width L (a few cm) and the dominant Lorentz acceleration (hundreds of m s−1 2 ) leading to several ms, noticeably longer than the < 3 ms ELM duration. Concerning heat transfer, the effect of convection and conduction is comparable for most of the melting event (extent as well as duration). Conduction being mainly responsible for resolidification, one can use the heat diffusion timescale δ 2 /α with δ the diffusion length (ca. pool depth, few hundred microns), α the thermal diffusivity (∼ 20 mm 2 /s within the W liquid phase [26]) for order of magnitude estimates.

Single transient W melt event
A representative ELM heat flux can be well captured by a Gaussian pulse with ∼ cm width and ∼ ms duration, see the caption of figure 1.
In the first regime of slow heat diffusion compared to fluid motion, dynamics are independent of heat transfer so that thermal analysis can be discarded. Phase transitions cannot be treated, thus the amount of melt should be prescribed. This simplification makes the treatment of non-linear dynamics beyond the shallow water approximation (including droplet ejection) computationally feasible. However, in this case, the melt depth and melt extent have to be introduced as external parameters and, more importantly, the expected long-term surface evolution strongly contradicts experimental evidence. In particular, the melt will be accelerated out of the pool by the Lorentz force, which is partially counteracted by viscous damping. Small-scale perturbations due to capillary forces will be of no practical importance in light of the large displaced volumes. Recalling that liquid metals wet their own solid excellently due to the strong interfacial bonding [27,28], the resulting macroscopic dynamics are rather intuitive: the liquid will move out of the pool and stretch over the adjacent solid until a thin film is formed; no mass accumulation near the pool edge is possible.
The other two regimes require the combination of the depth-averaged Navier-Stokes equations with heat transfer. In the case of fast heat diffusion compared to fluid motion; the convective term of equation (3) can be dropped, the heat conduction simulation domain remains undeformed and melt column height continuity (see equation (1)) can also be linearized in the limit of small surface perturbations compared to the melt thickness. The emerging equations are decoupled in the sense that the heat equation provides input to the continuity and momentum equations (interface position b 1 , melt thickness h), but fluid motion no longer affects heat conduction. Finally, the case of comparable heat diffusion and liquid motion timescales corresponds to the fully coupled MEMOS-U model described earlier. Figure 1 presents the temperature response to a Gaussian heat pulse and the associated surface deformation, for the cases of heat transfer independent of fluid motion (a) and of fully coupled heat transfer-fluid motion (b). The heat pulse generates a melt pool of depth ∼ 180 µm and 2 cm horizontal extent. After ELM incidence, cooling is mainly due to heat conduction, with contributions from thermionic emission, vaporization and thermal radiation. At early exposure stages, J × Binduced melt displacement is small and deformed surface profiles coincide. As displacement progresses, results begin to significantly deviate. In (a,b), strong heat conduction occurs only in the region where melting was originally induced. Due to the heat pulse spatial profile, the melt layer thins out towards the edges, leading to increased viscous damping near pool boundaries. Combined with the larger acceleration, that is also induced for a longer time, in the hotter and deeper central region, this results in a fluid velocity that decreases from its maximum value at the centre to zero at the pool edges. The final surface deformation corresponds to the black dashed line in (b). On the other hand, in (c,d), the melt moves onto the colder adjacent solid which provides an additional conduction channel for the deposited energy. In (d), liquid motion on the solid surface is arrested by resolidification, leading to a profile featuring a crater and mass accumulation near the edge.

ELM-driven W melt motion in AUG
We model recent W leading edge exposures to 2.5 T ELMy H-mode pulses (#33504, #33508, #33509) in AUG [12]. All pulses were similar in terms of input power and ELM properties. The absorbed heat flux, extracted from IR measurements in #33511, is used as an input. The average ELM frequency was 72 Hz, the ELM duration < 3 ms, the peak ELM heat flux within 0.5 − 1 GWm −2 , the inter-ELM heat flux within 30 − 100 MWm −2 . Strong spatio-temporal variations in the experimental ELM heat flux leave an imprint in melt production; each transient melt pool has a unique depth and position along the exposed edge (figure 2).
The evolution of the simulated melt layer thickness during pulse #33509 is illustrated in figure 2, without (a) and with (b) full coupling between heat transfer and fluid motion. The differences can be understood on the basis of the single event comparison in figure 1. The ELMs repeatedly strike roughly at the sample's centre. In figure 2(a), each new ELM hits a previously exposed region that remained hot or even liquid, eventually causing sustained melting near the end of the exposure. Overall, this leads to a significant overestimation of the predicted surface deformation and material pile-up at a location inconsistent with observations, see figure 3(a). In figure 2(b), incident ELMs strike a colder surface, preventing sustained melting. Hot material moves progressively along the sample towards weaker heat flux locations. At the end of the exposure, material reaches the sample's edge and flows out of the domain (indicated by black arrows), replicating the molten material ejection observed in the experiment [12].
We proceed with a detailed comparison between the fully coupled output and post-mortem observations. The final surface deformation is determined by non-additive accumulation of melt displaced from the individual transient pools depicted in figure 2(b). Figure 3 confirms that the main features of the final MEMOS-U surface profile are in close agreement with the experiment. As in the observations, the excavation region lies at the sample's centre and liquid metal is accelerated out of the transient melt pool towards the adjacent solid surface, where it resolidifies. The material pileup reaches all the way up to the sample edge, where some ejection occurs. The quantities of interest are the excavated volume (6.8-7.8 mm 3 in the experiment [12], 5.4 mm 3 in the simulations) and the ejected material (0.7-1.74 mm 3 [12] vs 1.5 mm 3 ).
Turning to dynamics (equation (2)), melt motion is driven by the Lorentz force density (J × B) and is mainly counteracted by depth averaged shear stresses (−3µU/h 2 ) which strongly depend on the melt thickness that is regulated by resolidification. Melt velocities are determined by the transient balance between the above, since contributions from the remaining non-inertial momentum terms are insignificant. This leads to typical velocities and accelerations of 1 − 2 m s −1 and up to 700ms −2 , respectively.

ELM-driven W melt motion in JET
The misaligned W lamella was subject to 9 consecutive pulses with surface deformation mainly inflicted during pulse #84779 [11], modelled here. The incident heat flux was reconstructed from IR thermography data [29]. The typical ELM duration was < 3 ms, the average ELM frequency 30 Hz, the intra-ELM heat flux 0.5 − 1.5 GWm −2 , the inter-ELM heat flux 50-200 MWm −2 . Experimentally extracted heat loads induce superficial melting incompatible with observations [11]. A 10% up-scaling, that lies well within the estimated uncertainties, results in the best agreement. As illustrated in figure 4, the MEMOS-U deformed surface profile is in good agreement with observations [11] concerning the crater position (2.0 cm from the edge), its depth (0.6 mm) and extension (2 cm) as well as the melt pile-up location (0.7 cm from the edge), its height (up to 1.7 mm) and its extension (1 cm). Since the heat load and sample temperature evolution are very similar to the AUG exposures, the physics discussions remain valid.

Disruption-driven Be melt motion
Upward VDEs during the first two JET ITER-Like Wall (ILW) campaigns led to melting of Be upper dump plates (UDP), each of which comprises 8 individual tiles [30]. Most of the damage was inflicted during ILW2, with 29 pulses identified as the main culprits, largely featuring deliberate VDEs [13]. Surface deformation was mostly observed on the low-field side DP-8 plates [13]. Melt is produced on the wetted side (Face 1), moves over the corner onto the shadowed side (Face 2) and re-solidifies. The emerging pile-up of frozen Be liquid resembles an inverted 'waterfall' [13,30].
Experimental input involves integrated quantities; the total VDE induced current and total UDP conducted energy. Thermocouple data lead to 20-30 kJ for the energy conducted by a DP-8 plate over the pulse duration [13]. The observed 30-60 cm 2 damaged area can be used for heat flux estimates. For square uniform heat pulses, assumption of energy deposition during the thermal quench (∼ few ms) translates to heat fluxes of a few GW m −2 that would lead to intense vaporization, but shallow melt layers (≲ 100 µm), contradicting observations [13]. In contrast, energy deposition over the current quench (CQ) timescale of few tens of ms translates to ∼ 0.1 GW m −2 which suffices to reproduce the observed melt pools. Typical CQ halo current densities are ∼ 0.1 MA m −2 [30].
Due to the shallow water assumption, MEMOS-U simulations are set up in simplified geometry, with the shadowed UDP face rotated by 90 • ( figure 5). This does not affect heat transfer, which is insensitive to such geometrical effects on these timescales given the ∼ 1 mm heat diffusion length. The Be melt passage across the sharp edge is not captured but its subsequent motion on the shadowed face is correctly described. Figure 5 shows the post-exposure pattern after the deposition of 0.04 MJ over 35 ms according to a non-uniform profile detailed in the caption. The 50 kA m −2 halo current density is assumed uniform, normal to Face 1 (due to the large scale difference between wetted area and melt depth) and constant over the pulse. Thermionic emission from Be is negligible owing to its low melting point and high work function. The waterfall extension (several cm) and total volume (∼ 1 cm 3 ) can be matched to observations by tuning the current density (responsible for melt acceleration) Contrary to ELM-induced events where motion mainly occurs after each pulse termination, longer, continuous VDE heat fluxes (resulting in deeper pools up to 500 µm) lead to uninterrupted melt production while Lorentz acceleration removes the hotter upper layers. This strong energymomentum transfer coupling enables a more efficient heat conduction and melt generation compared to a stagnated pool, which would exhibit strongly increasing surface temperature and enhanced vaporization. Balance between the resolidification rate and the uninterrupted melt supply (dictated by heat transfer in the wetted area) determines the moving layer thickness and final deformation profile. Similar to ELMinduced events, the J × B force (here due to the halo current flowing through the liquid pool) drives motion, depth averaged shear stresses damp it and other non-inertial terms are negligible. The melt velocity, dictated by transient momentum balance, is within 1 − 4 m s −1 and acceleration within 100-400 m s −2 .

Summary
Resolidification-controlled liquid dynamics are ubiquitous in contemporary tokamak melting events due to the small ratio of the plasma-wetted area to the overall plasma-facing unit area, which results in localized melt pools bounded by progressively colder solid surfaces. In such scenarios, heat transfer is strongly entangled with fluid motion, forbidding simplified descriptions based on timescale separation arguments. Comparison with experiments has demonstrated that the MEMOS-U code adequately accounts for the key physical processes involved in fast transient melt events and has lent confidence to its use as a predictive tool for ITER. The modelling reveals that despite significant differences between ELM-induced W melting in the divertor and disruption-induced Be melting in the main chamber (in terms of heat loading, bulk current source, material composition), metallic melt motion is governed by similar physical mechanisms.