OpenFOAM modeling of beryllium melt motion and splashing from first wall in ITER

Beryllium (Be) is a material which will be used as a plasma facing component in ITER due to its unique properties of high thermal conductivity, low density, and high strength. However, under extreme conditions of high temperature and pressure, Be can melt at the surface of tiles and molten droplets can be ejected into the reactor leading to disruption of fusion plasma. The pressure, mass density, velocity of Be vapor, and variations of temperature at the melt layer interface can influence the splashing of Be melt. The Computational Fluid Dynamics (CFD) model based on the OpenFOAM toolbox, a free open source CFD software package, was developed to treat the coupled flow of liquid Be metal and its vapor. The vapor-melt interface is modeled using the volume of fluid (VOF) approach implemented in the interCondensatingEvaporatingFoam solver that solves the continuity, momentum, heat conduction, and VOF equations. This CFD model is capable to predict the hydrodynamic effects of Be vapor on the melt layer motion, splashing, non-linear growth of melt waves, and ejection of molten droplets. The modeling accounts for the effects of thermal, viscous, gravitational, and surface tension forces at the vapor-melt interface. In this research, we used the interCondensatingEvaporatingFoam solver to simulate the effects of Be phase change and the development of melt motion with the droplets ejected from the surface. The CFD model accounts for inter-phase change between Be liquid and Be vapor. The evaporation model was validated against the Stefan phase-change problems. The influence of heat and mass transfer across the vapor-melt interface on melt layer stability is also investigated. The results provide an understanding of how the rate of phase change affects the development of melt structures and waves at the vapor-melt interface.


Introduction
The examination and understanding the material erosion and migration is one of the important factors for a sustainable and safe operation of the ITER tokamak. Carbon-based first-wall materials have been studied for a while in tokamak reactors [1]. However, the graphite first wall brought several problems. Tritium as one of the fuels for fusion reaction, together with deuterium, can be highly retained in the graphite [2]. As a result, high fuel retention and severe erosion occur on the carbon-based material [3,4]. The flows of these particles are mainly directed towards the inner divertor region where their deposition takes place [5]. In this process, particle transport continued to take place, which led to the degradation of plasma facing component (PFC) due to high fuel retention in the deposited layers, and this has significantly affected the efficiency of carbon-based materials for ITER applications [6].
The other materials as candidates for a first wall have been investigated. The selection of materials was greatly assisted using the plasma-wall interaction computer codes such as ERO [7] and WallDYN [8]. The replacement of PFC, as in ASDEX Upgrade from graphite to tungsten [9], and recently in JET from the carbon fibre material to beryllium (Be) tiles was performed to form a first-wall in the main chamber [1]. In these experiments, Be was examined due to its low atomic number and tolerance of its impurity concentration in the fusion plasma, which

OpenFOAM CFD model
OpenFOAM is an open-source toolbox with extensive multi-physics implementations [57]. It is under active development with the capabilities of commercial CFD software. It allows both 2D or 3D structured and unstructured mesh and parallel running. The solvers in OpenFOAM are implemented in C++ programming language with the source code freely available. These solvers can be used as templates and modified for further development. Many ordinary and partial differential equations can be easily converted into C++ code in OpenFOAM. For our research related to the motion and splashing of Be melt layer, interCondensatingEvaporatingFoam solver was used. It treats the flow of two incompressible and immiscible fluids with the addition of energy equation and the model for phase change (evaporation-condensation) between the melt and its vapor.

Governing equations
The interface between the melt and its vapor phases is resolved using the VOF approach. The location of interface is defined by a scalar function, melt phase volume fraction, represented by where V is the volume of cell and V m is the volume occupied by melt in this cell. Subscript m denotes melt. Melt phase volume fraction is a = 1 m in cells filled with melt, a = 0 m in cells filled with vapor, and between these two values (0 and 1) at the interface. Vapor phase volume fraction can be easy calculated as a a = -1 .
v m Subscript v denotes vapor. The equations of the conservation for mass, momentum, energy, and volume fraction solved by interCondensatingEvaporatingFoam solver are presented below.
Assuming the change in density is negligible, the continuity equation for incompressible fluids can be written as is the velocity vector with u m and u v are the velocity vectors of melt and vapor, respectively.
The momentum equation balances the pressure, viscosity, surface tension, and gravitational forces acting on a fluid element and can be represented by where s is the surface tension and the interface curvatures of melt and vapor in terms of the gradient of a m are given by For numerical convenience to account for the hydrostatic pressure contribution that is important in the multiphase cases, the modified pressure r =r ( · ) g h p p gh rather than the total pressure p is implemented in the solver. Here g is the gravity vector and h is the position vector. The gradient of r p gh can be written as gh Therefore, the term r - + g p in equation (2) can be replaced by - r p .

gh
The energy equation is formulated in terms of temperature T as r r v v is the effective thermal conductivity with melt and vapor conductivities k m and k , v + a a = C C C p m p m v p v is the specific heat capacity with melt and vapor specific heat capacities C pm and C , pv and  Q is the heat source term accounting for the changes in temperature during the phase transformation process. It will be detailed later.
In the VOF interface capturing method, the transport equation for the volume fraction of melt is given by a a a a ¶ ¶ where an artificial compression term a a ·( ) u m v c is introduced in equation (4) to limit the smearing of the interface [59]. It contributes only in the region of the interface due to a non-zero factor a a .
m v The compressive velocity u c is calculated in the normal direction to the interface. The source term  M accounts for the rate of mass transfer during evaporation and condensation and can be expressed as [60] where  m e and  m c are the mass transfer rates per unit volume due to evaporation and condensation between phases. Subscripts e and c stand for evaporation and condensation, respectively. Dot above  m means the rate. The rates of mass transfer at constant pressure are described by the Lee's model as [61] where T sat is the saturation temperature, C e and C c are the coefficients characterizing the intensity of the mass transfer in the evaporation and condensation processes. In the Lee's model, the mass transfer at constant pressure depends on the saturation temperature. Also, the values of coefficients C e and C c depend on material properties and should be selected from the available experimental data. The theoretical expressions are also available for their evaluations [62,63]. The source term  Q in equation (3) is expressed as the product of latent heat (latent energy or heat of phase transformation) H mv and the rate of mass transfer The interCondensatingEvaporatingFoam solver employs the finite volume discretization scheme to solveequations (1) to (4). The pressure-velocity coupling PIMPLE scheme which is a combination of PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) is used for solving iteratively continuity and momentumequations (1) and (2). A first-order implicit Euler scheme is used to discretize the time derivatives. The time step is controlled by setting the maximum Courant number below 0.5. A second-order Gaussian linear integration scheme with van Leer limiter is used for discretizing the divergence and gradient terms in equations (1) to (4). The solution converges when the residuals reach 10 −9 for the pressure and 10 −7 for the velocity and temperature corrections.

Temperature-dependent coefficients
The following equations were used to calculate the thermodynamic and transport properties of liquid and vapor Be at different temperature [63][64][65][66][67]. The vapor pressure of Be at the melt surface for the given temperature is evaluated using the parametric equation where T is used in units of Kelvin in all formula. It is about ∼ Pa 5 at the melting temperature of 1560 K. The specific heat at constant pressure can be calculated using = +´-- For molten and vapor Be, their mass densities are obtained using The mass density of molten Be at melting point is / kg m 1690 . 3 The perfect gas law is used to calculate the mass density of Be vapor. Dynamic viscosity of Be melt is determined using Kinematic viscosity for Be vapor treated as an ideal gas can be computed directly using Newton's law where m is the mass of Be atom, k B is the Boltzmann constant, s » 0.228 nm Be is the hard-sphere diameter of Be atom, and W m is the collision integral (assumed to be equal to 1 for an ideal gas (non-interacting atoms)). The viscosity of vapor is independent of pressure and depends only on temperature.
With two sets of measured surface tension data [65], one can find the surface tension gradient as they are in linear relation. The surface tension of Be melt is calculated as Derivation of this equation for surface tension can be found in appendix. These parameters are temperature dependent, which allow to determine properties of liquid and vapor Be, which vary with temperature at the surface of the PFC. The calculated parameters in the range between the melting and boiling temperature of Be are shown in table 1. Finally, it is useful to provide some constants for Be material such as the latent heat of evaporation is / 53.55 kcal mol, the heat of sublimation is / 76.56 kcal mol, and the heat of fusion is / 2.8 kcal mol. The other basic properties of Be can be found in [65].

Setup of geometry, initial and boundary conditions
For the geometry setup, the computational domain with dimensions of 10 cm by 3 cm was created (figure 1). Such dimensions are chosen to provide a clear visualization of the splashing motion near the PFC surface as well as the molten Be pool. The domain was discretized into a mesh. At = t 0, the volume fractions of Be melt in a pool with a length of 7 cm and vapor in the rest of domain are illustrated. The thickness of Be melt is approximately 2 mm according to JET reports [68]. In the experiment conducted in the DIII-D tokamak, the plasma velocity was found to be~/ 0.5 km s in the poloidal direction and~-/ 10 20 km s in the toroidal direction during ELMs [69]. Hence, Be vapor velocities of 1, 5 and / 10 km s are used to investigate their effects on the splashing. The melt velocity is set at / 2 m s as the TEXTOR experiments show the motion of melt with the velocity approximately / 1.7 m s [69]. As described earlier, properties of Be melt and its vapor, which provided in table 1, are temperature dependent and calculated in section 2.2. For example, at temperature of 2019 K, vapor pressure and mass density at the surface are calculated to be 10 Pa 3 and -0.000537 kg m , 3 respectively. For each simulation, a

Results
The interCondensatingEvaporatingFoam solver was benchmarked against the wave oscillation problem [70] without considering the phase change effects. For this purpose, the coefficients of evaporation and condensation are set to zero in interCondensatingEvaporatingFoam solver. In order to ensure the interCondensatingEvaporatingFoam solver is valid, Stefan phase-change problem is also used to model Be melt motion and splashing with account for the evaporation and condensation effects. The solver is then used to model Be melt splashing from a pool.

Wave oscillation benchmark test problem for interCondensatingEvaporatingFoam solver
In order to ensure the validity of simulations, the comparison between numerical and analytical data of kinetic energy of a dimming wave was performed. Liquid with a wavelength of 10 −4 m will initiate and leave to dim without any external force applied. As the wave oscillate, its kinetic energy is lost through the process and the numerical results of the simulation can be used to compare with the analytical results. The setup and properties of the wave is shown in figure 2 where the wavelength is set as 10 −4 m.
As the wave oscillates, its kinetic energy at different time moments will be calculated during the simulation. These results are used to compare with the theoretical values for validation. The rate of decay of kinetic energy KE due to viscous effect is calculated using the following equation [70]  r sk To improve the accuracy, the different grid resolution in term of ratio of wavelength to grid step λ/dx was used. The results of the comparison are shown in figure 3. It is seen that the calculated kinetic energy ratio for λ/dx = 10 and 20 converge with the exponential decay of the theoretical decay rate. Energy peaks of the calculated oscillations display close agreement with the theoretical curve. Hence, it provides the validation of interCondensatingEvaporatingFoam solver for our simulations of Be melt splashing without taking into account the phase-change effects.

Be melt motion and splashing without taking into account phase change
Motion of the Be melt was investigated at different time moments. First, we study the case when properties of Be melt are kept at temperature of 2717 K, which is its melting temperature, and its corresponding pressure is calculated to be 10 5 Pa (equation (9)). Vapor velocity and mass density used in the simulation were set to 10 km s −1 and 0.0452 kg m −3 , respectively. The volume fraction of the Be-melt is shown in figure 4    The time evolution of velocity field, which is affected by the motion across the surface, is also studied at temperature of 2717 K. Figure 5 displays the evolution of velocity field with time. At time of 0.20 ms, the disturbance of melt layer and appearance of a ligament at upstream have caused a large vortex to be formed after the edge of the pool (top panel). As the Be-melt splashes out of the pool, velocity region at the surface are indicating more uniform flows. However, the void created by the washed off Be melt still results in a significant velocity difference at a boundary~1 cm above the surface (bottom panel).
As a wide range of vapor velocity was observed in the DIII-D tokamak experiment [71], it suggests the splashing motion of melt at different vapor velocities needs to be examined. Figure 6 shows the motion of melt for vapor velocity of 1, 5 and 10 km s −1 at time of 1.12 ms. Under melt temperature of 2717 K, the evaluated pressure of 10 5 Pa, and mass density of 0.0401 kg m −3 , the surface waves develop much faster at higher vapor  speeds. At 1 km s −1 , single Be-melt ligament is formed, and relatively smooth waves occurred. For 5 km s −1 , Bemelt is splashed from half of the pool, whereas the collision and coalescence of long ligaments are seen in the other half of the pool. When the vapor velocity is 10 km s −1 , Be-melt layer is completely removed from the pool at 1.12 ms. As a result, it implies that high-speed flow of vapor can significantly thinner large portions of melt layer at a shorter period.
It was also observed that Kelvin-Helmholtz instabilities are developed, where liquid ligaments with droplets are carried away by the Be-vapor flow due to Be-vapor pressure on the melt surface. Such phenomenon will enhance the rate of splashing of the melt layer away from the surface and into the reactor chamber.
Other than velocity of the flow, vapor density on the Be-melt layer also plays an important role in development of the splashing. Figure 7 shows the snapshots of the melt structure for vapor velocity of 5 km s −1 at time of 1.62 ms. Different temperatures and their respective pressures were set in the runs. At T = 2019 K and vapor pressure of 10 3 Pa that corresponds to vapor mass density of 0.000537 kg m −3 , a single, large Be-melt vortex is developed and there is a smooth melt layer in the rest of the pool (top panel). However, as the temperature is set to T = 2316.5 K which resulted in vapor pressure of 10 4 Pa and mass density of 0.0047 kg m −3 , the irregular surface with Be-melt protrusions is formed (middle panel). Accumulation of the melt is also observed at the edge of the pool. At T = 2717 K and 10 5 Pa that corresponds to vapor mass density of 0.0401 kg m −3 , Be-melt is completely splashed out from the pool (bottom panel). The high vapor pressure has caused the splashing to be relatively close to the surface with less influence of Kelvin-Helmholtz instability.

Stefan phase-change problem for validating interCondensatingEvaporatingFoam solver
The one-dimensional Stefan phase-change problem [72] is selected for validating the Lee's phase change model implemented in interCondensatingEvaporatingFoam solver. Among the different methods that were suggested for simulation of phase change between liquid and vapor, VOF method is chosen [73][74][75][76][77][78]. As VOF method contains inherent mass-conservation property, it allows the interface of different states to be more distinctive [76,77]. Figure 8 shows the geometry setup of the Stefan phase-change problem. Computational domain is constructed in dimension of 1.6 cm by 0.16 cm. On the left boundary, a hot wall of 383.15 K with vapor insults the liquid from the heated wall. A cold wall of 373.15 K is set as the right boundary condition. At initial state,  In order to validate the numerical results with the analytical values, the following equations were used to calculate the interface position at the given temperatures [79] = - Equation (18) is used to obtain the Stanton number (St), which is the ratio of heat transferred into a liquid to the thermal capacity of liquid, at the given specific heat, latent heat of vapor, wall temperature and saturation temperature. By substituting the β value from equations (19) into (20), it gives the position of the interface at a specific thermal conductivity of vapor, d v at any given time. This implies that the boundary between vapor and liquid region changes as time progresses. Equation (21) calculates the values of temperature at any specific time and position within the setup domain. These will provide sets of analytical values to plot interface position and temperature against time. For computations, the geometry set in figure 8 was used in interCondensatingEvaporatingFoam solver together with density, viscosity, specific heat, latent heat of vapor, liquid heat conductivity, surface tension and C e as 1 kgm −3 , 10 −7 m 2 s −1 , 1000 kJkg −1 K −1 , 1 Wm −1 K −1 , 10 6 kJkg −1 , 0.1 Nm −1 and 1000 s −1 respectively. Different gas heat conductivities of 0.1 and 0.01 Wm −1 K −1 are set as the variables for comparison. To achieve higher accuracy, mesh size is set to be 16 μm. Computational run starts at equilibrium with left and right boundaries set as wall with zero initial velocity within the system. Liquid state evaporates into vapor and results in the rightwards movement of the phase boundary away from the heated wall. The progression of vapor-liquid interface in the simulation of Stefan problem is shown in figure 9 which displays the evolution of volume fraction of vapor region (blue) and liquid region (red) along x-axis with time. From the heated wall at the left boundary, the volume fraction is uniform throughout the vapor. Along the x-axis, there is a sharp decrease in the expansion of vapor region to liquid phase in a relatively distinctive interface. Temperatures of the phases are constant throughout their respective regions. Temperature of the vapor region is fixed at saturation temperature and variation of temperature only occurs at the fine domain bounded by the two phases. Figure 9 also shows the decrease in expansion rate of the vapor.
Numerical values of interface position and temperature distribution at different time steps can be extracted from the interCondensatingEvaporatingFoam simulations. Together with the analytical results, obtained by using equations (18) through (21), the interface position as a function of time for different thermal conductivities is illustrated in figure 10. It shows a comparison between the results of the theoretical calculations and simulations using interCondensatingEvaporatingFoam.
At a lower thermal conductivity of 0.01 Wm −1 K −1 , interface position tends to achieve an equilibrium during a shorter time period. The red solid line and green dots, which represent the analytical and numerical data respectively, show excellent agreement between the results. Despite there is decrease in rate of expansion of the vapor, the curves for higher thermal conductivity still show room for expansion even after 100 s. This implies that vapor with higher thermal conductivity will expand and occupy a larger region between the superheated wall and the liquid. Overall, the analytical results can be used to validate the numerical calculations as they are in a relatively high agreement with each other. Figure 11 displays temperature profiles along the x-axis at different times of 20, 60 and 100 s. The analytical curves are in good agreement with the data points which were generated using the numerical method. The constant temperature of heated wall at 383.15 K ensures consistent heat transfer into the vapor region. As the time progresses, the temperature gradient is lowered, and this implies the rightward movement of the interface position between the two phases. The observation also testifies the shift of interface position with time which was displayed in figure 10.
Apart from the benchmark for evaporation, Stefan phase-change problem can also be used for validating the condensation. For condensation, the geometry setup will be inverted where liquid will separate the vapor from the cold wall. Before the simulation run, the liquid and vapor regions are in quiescent equilibrium. As the time runs, condensation of the vapor phase into liquid phase will cause expansion of the liquid region. The interface is at the saturation temperature and the liquid phase will push the vapor away from the cold wall until it reaches an equilibrium position where the rate of evaporation of liquid matches the condensation rate of vapor. This behavior shares a similarity with the evaporation case which was shown in figure 9. After validation of Stefan phase-change problem, the inter-phase change model implemented in interCondensatingEvaporatingFoam solver can be applied to the Be-melt.

Be melt motion and splashing with taking into account phase change
In this section, the Lee's phase change model, which was validated against the Stefan phase-change problem, is applied to the Be-melt with the initial setup shown in figure 1. Since the temperatures within the computational domain are unchanged, the equations which are used to calculate the properties of Be remain the same. By including the phase change, different C e constants are used to determine its effect on the splashing motion. The saturation temperature is set to 2744 K, which is the boiling point of Be. This allows to closer mimic the realistic environment of melt behavior within fusion chamber.
For the first case, Be melt is kept at temperature of 2717 K and pressure is calculated to be 10 5 Pa. Vapor velocity, temperature, pressure and mass density used in the simulation were set to 10 km s −1 , 2744 K, 1.15 × 10 5 Pa, and 0.0452 kg m −3 , respectively. Different melt topologies for different coefficients of evaporation are

Conclusions
The main objective of this research is to study the splash motions of Be melt under different hydrodynamic conditions and the influence of the phase change. Validations were done against the wave oscillation and Stefan phase-change problems to provide the benchmarks for interCondensatingEvaporatingFoam solver.
Temperature and other respective properties were set for the molten Be. As the simulation initiated, different splashing motions were observed for the cases with and without the phase change. With no phase change involved, discrete ligaments are formed on the surface, while the cases with phase change show that the Be melt is nebulized during the splashing. The significant portion of Be melt is splashed out of the pool with time. With increase in vapor velocity, there is a faster growth of waves and ripples on the surface, and it takes less time for splashing of Be melt from a pool for the cases which include the phase change. With increase in mass density, ligaments are compressed closer to PFC surface and elongated by the vapor flow as they are removed from a pool. By comparing figure 4 and figure 13, it is observed that the molten layer is localized closer to the surface when it splashes out of the pool. Hence, it suggests a higher portion of the molten Be is ejected into plasma when the phase change is considered. Depending on the vapor velocity and mass density, complete removal of the Be melt from a pool can take around 1-3 ms. With the phase change, the rate of evaporation depends on the coefficient of evaporation. As vapor velocity increases, there will be a faster rate of depletion of the Be melt layer as thin film of inter-phase region above the melt is removed from the surface. Lastly, according to the simulations and depending on the vapor velocity, complete removal of the Be melt from a pool will occur at a faster rate which only takes up to  0.5 ms. Through these simulations, they provided the estimated duration for complete removal of the 2 mm depth molten layer. The influence of phase change on Be melt splashing has significantly decreased its removal time from 1-3 ms to less than 0.5 ms. The topological structures of the splashed melt also show less discrete waves when the phase change is taken into account.