Melting Penetration Simulation of Fe-U System at High Temperature Using MPS_LER

Melting penetration information of Fe-U system is necessary for simulating the molten core behavior during severe accident in nuclear power plants. For Fe-U system, the information is mainly obtained from experiment, i.e. TREAT experiment. However, there is no reported data on SS304 at temperature above 1350°C. The MPS_LER has been developed and validated to simulate melting penetration on Fe-U system. The MPS_LER modelled the eutectic phenomenon by solving the diffusion process and by applying the binary phase diagram criteria. This study simulates the melting penetration of the system at higher temperature using MPS_LER. Simulations were conducted on SS304 at 1400, 1450 and 1500°C. The simulation results show rapid increase of melting penetration rate.


Introduction
This study was made to aid in the safety evaluation of hypothetical meltdown during severe accident. At first phase, high temperature will cause swelling and bending of materials. This condition will provide contact interaction among neighboring materials. At high temperature, surface contact of two materials may decrease their melting point due to eutectic reaction. Melting penetration rate of uranium into iron at high temperature may be important information for computer simulation to simulate melting behavior in nuclear power plant during severe accident.
For SS304, TREAT experiment (Fe-U system) only presented results at and below 1350°C [1]- [3]. Previously developed MPS_LER code has been created to have capability to simulate eutectic reaction. The TREAT experiment has been use to validate the code.
The objective of this study is to predict the melting penetration rate using modified MPS method of SS304 with liquid uranium above 1350°C.

MPS method
The MPS method is the first particle method developed for analyzing incompressible media such as water and metal important for industrial application [4]. The governing equations of MPS method are discretized in terms of particles. A particle interacts with its neighbours using the weight function : where is the distance between two particles and is the effective radius ( Figure 1). To calculate a weighted average, particle number density is to be calculated as, Where is the coordinates of each particle. The governing equations of velocity field are as follows: Equation (3) is the law of mass conservation. This shows that density change is zero in every time. Equation (4), Navier-Stokes equation, is the law of momentum conservation. is the kinematic viscosity. The left side of equation (4) is Lagragian differential of vector velocity. The first term of right hand side is the pressure term then viscosity and gravity, respectively. Implicit calculation of gradient model is applied to pressure. The viscosity is calculated using explicit calculation of Laplacian model. Differential operators, such as gradient and Laplacian, are represented by the following particle interaction models.
Where and ! are the number of spatial dimensions and initial particle number density, respectively.
Differential operators in the governing equations are substituted by the particle interaction models, and then motion equations of the particles are obtained. The particle number density has two meanings. First, it is the normalization factor of the weighted average. Second, it is the value proportional to the fluid density. Since the density should be constant in incompressible fluid dynamics, the particle number density is required to be constant. The constant is denoted by ! . Meaning that the particle number density needs to be kept equals to ! in order to analyze incompressible flow.
To keep the particle number density constant, a semi implicit algorithm is used in the MPS method. This algorithm included explicit and implicit steps. The first is the explicit step in which the forces except for the pressure are calculated, such as gravity, viscosity and surface tension. After the step, the particle number density does not equal to ! .
This issue should be solved in the second step. The Poisson equation of pressure is implicitly solved in the second step. The Poisson equation of pressure is deduced from the implicit mass conservation equation and the implicit pressure gradient term.
, where * is the temporal particle number density after the explicit step. The right hand side of equation (6) represents the source term which is the deviation of the temporal particle number density from ! . The Poisson equation is discretized to a symmetric matrix equation using the Laplacian model (equation (3)). The matrix equation is solved by the conjugate gradient (CG) method. In fluid dynamics, the list of neighboring particle is updated in each time step. Figure 1 Relative coordinates, effective radius and initial minimum particle distance.

Model development
Basically, diffusion coefficient of species and temperature dependent of mixing ratio. In this model, the definition of eutectic melting is the condition where melting occur below its melting point due to mixing of the two materials. The modeling of eutectic phenomenon is consisted of two steps. The first step is the calculation of mass transfer among particles. The solution is provided by solving Fick's second law of mass diffusion. The second step is applying boundary based on binary phase diagram as chemical reaction.
The first step is the utilization of Fick's second law as governing equation for mass diffusion by where m is mass and D is diffusion coefficient of material. Mass diffusion term of right side of the equation (equation 9) is calculated by MPS Laplacian model explicitly, equation (10). And Laplacian radius is written as equation (11).
where ! !!! is mass of particle at the time step + 1, ! ! is mass of particle at the time step , ! ! is mass of j particle at the time step t, ! is position vector of i particle, ! is position vector of j particle, d is the number of space dimensions, and ! is initial particle number density of inner particle which does not include particle detected on free surface in the interaction zone. Diffusion coefficient (D) is obtained from experimental data. Experiments were conducted at several temperatures to record diffusion coefficient. However, for materials that has no experimental data that may due to complex preparation of experiment, molecular dynamic method is one the solution to obtain diffusion coefficient. The second step is correlating mass diffusion result in a particle with eutectic criteria from binary phase diagram. In the phase diagram, the criteria are given in temperature and mass fraction. MPS_LER assumed that each particle is mixture bowl where eutectic reaction takes place. This simplification should be taken because MPS method works in mesoscopic level, thus atomic movement cannot be simulated well. Each particle brings information of its mass density based upon the material and particle size of the particle. Mass diffusion occurs between the particle and its surrounding particles. The MPS_LER has been validated with TREAT experiment and Pb-Sn experiment conducted by CRIEPI [5], [6].

Simulation condition
In this simulation, the minimum particle distance was about 0.0001 m. The test section thickness of the MPS geometry was about 500 µm. The total number of particles for developing the test geometry was 42164. Figure 2 shows the 2D calculation model for TREAT experiment.
Diffusion coefficient was obtain from first principle molecular dynamics (FPMD) code VASP simulation [7]. The calculation was conducted on metal fuel U-Pu-Zr and cladding Fe. Detail information of diffusion coefficient calculation of SS304 and uranium is explained in previous work [5].  Fig. 2 shows the results of MPS_LER calculation on melting of the SS304 test section. The results indicate that after particular time depending on temperature, the outer layer particles of the test section were transforming into molten particle. This suggests that melting is taking place to the outer particles. Fig. 3 shows the melting penetration rate predicted from the geometry of TREAT experiment. Three temperature conditions were simulated, i.e. 1400, 1450 and 1500 ºC. At 1400ºC, the penetration rate increased to about 300 µm/s. The increased of penetration rate is correspond to mass fraction of the system and increase of diffusion coefficient as temperature increases. At this temperature, the fraction of uranium over total mixture is only 0.23. On the other hand, the diffusion coefficient of uranium into iron is about 8.00E-9 m 2 /s while the number for iron into uranium is 2.10E-9 m 2 /s. At  [7]. While diffusion coefficient of iron into uranium is about 2.24E-9 m 2 /s. Penetration rate at 1500ºC increased rapidly to about 1250 µm/s. This is due to decrease of mass fraction and increase of diffusion coefficient. The overall pattern of melting penetration rate increases exponentially.

Conclusion
Moving particle semi-implicit (MPS) code was used to calculate eutectic interaction of solid-liquid system by of the U-Fe, i.e. SS304, system. The prediction of melting penetration may give important information for future investigation on severe accident, since there is no empirical data above 1350°C. It was found that the melting penetration rates were rapidly increased, especially at temperature near melting point of iron.