Anni Nuril Hidayati? Abdul Waris ? Asril Pramutadi Andi Mustari ?Dwi Irwanto ? Nur Asiah Aprianti

Abstract In the case of a severe accident involving nuclear reactors,an important aspect that should be considered is the leakage of molten material from the inside of the reactor into the environment. These molten materials damage other reactor components,such as electrical tubes,grid plates and core catchers.In this study,the moving particlesemi-implicit(MPS)method is adopted and improved to analyze the twodimensional downward relocation process of moltenWood’s metal as a representation of molten material in a nuclear reactor. The molten material impinges the Wood’s metal plate(WMP),which is mounted on a rigid dummy stainless steel in a cylindrical test vessel. The breaching process occurs because of heat transfer between the molten material and WMP.The formed breach areas were in good agreement with the experimental results, and they showed that the molten Wood’s metal spread above the WMP. The solid WMP fraction decreased with time until it reached the termination time of the simulation. The present results show that the MPS method can be applied to simulate and analyze the downward relocation process of molten material in the grid plate of a nuclear reactor.

Keywords Heat transfer ?Moving particle semi-implicit ?

1 Introduction

In recent years, melt leakage behaviors have been increasingly recognized as an important, worldwide nuclear technology concern globally, especially when a severe accident occurs. The molten materials that leak during an accident may relocate to other locations inside the nuclear reactor and damage other components [1]. The damage can be avoided by quenching the molten materials using coolants, such as sodium, to decrease their temperature and fragment them.Consequently,some fragments or freeze materials plug or clog and obstruct the other components [2]. Furutani et al. conducted experiments to determine the erosion rate of a solid plate impinged by a liquid jet using paraffin and Wood’s metal material [3].These materials were selected because of their low melting points. The conjugated nature of melt breakup and steam generation has been experimentally investigated using the non-boiling conditions of Wood’s metal [4], which explains why the Kelvin–Helmholtz instability caused jet breakup.Similarly,Wood’s metal has been used to analyze the molten behavior when it is poured into a stainless steel tube and injected into the coolant through a nozzle[5,6].It was found that the frozen layer of metal adhered to the structure and its shape was thin and wide in the air coolant experiment.However,in the water coolant experiment,the frozen metal was broken and fell down as debris. Additional research into molten behavior has been conducted by the Korea Atomic Energy Research Institute by losing some feed water and coolant without safety injection.These scenarios were referred to as level 1 probabilistic safety assessments[7].Recently,molten Wood’s metal was used by Sudha et al. (2018) to investigate the downward relocation to a plate with the same material as the molten Wood’s metal [8].

The MPS method has been used to simulate the free surface and multiphase flow of incompressible media [9],and it was developed by Koshizuka and Oka in 1996.This is a numerical method that does not use meshes, but uses particles as a part of the simulation. The Navier–Stokes equation is mainly used for dynamic particle calculations.The equation is transformed into particle interactions within an effective radius,and it has been used to simulate the penetration rates of molten uranium into steel test sections [10]. It was found that there was a small discrepancy between the result and the experiment because of the assumption of the diffusion coefficient, which was applied in the calculation.Li performed a two-dimensional(2D)MPS calculation to predict the molten behavior of the Fukushima Daiichi Unit 1 [11]. It focused on the lower head tube rupture, resulting in debris being released into the other part of the reactor. Similarly, Chen used it to investigate the melt-freezing process of uranium dioxide in a tube, based on several experiments [12]. In previous studies, the MPS method was used to simulate the process of vertical fluid descent on static fluids such as water and cooking oil using pipes. This has been validated using experimental data [13, 14]. Moreover, the stratification of these fluids was simulated using a 2D MPS code [15].These previous studies were simulated without heat transfer calculations between particles in the geometry because they were conducted at 300 K.

Recently, a molten metal relocation experiment was conducted to investigate the performance in the absence of coolant [8]. The results showed how molten material damages the supporting system when a severe nuclear accident occurs. Although it was not an exact scale model of a nuclear accident, some physical phenomena may be experimentally simulated. However, some behaviors observed when the molten material interacts with WMP were not explained in detail.The physical behavior can be approximated using a numerical simulation.Therefore,this study aims to use the 2D MPS method to analyze the WMP downward relocation process. The focus of the numerical simulation is on the relocation process of molten Wood’s metal and the impact on hole formation of the WMP. The heat transfer process was included to simulate the melting process of the WMP owing to the molten impingement of Wood’s metal. The simulation results were compared and validated by the experimental results. Moreover, some predictions of the behavior of molten Wood’s metal can also be explained numerically.

2 Methods

2.1 MPS method

The MPS method is a fully mesh-less Lagrangian method for simulating incompressible media,such as solids and liquids [9]. The Lagrangian fluid dynamics are represented by particles that possess their own physical properties. For example, the velocity field is not a function of position, as is the case with Eulerian fluids. MPS has two specific definitions,namely the radius of interaction reand the weight function w(r). The radius relimits the interaction of the particles,as shown in Fig.1.Moreover,it can be implemented as the physical influence of a particle on others within a certain distance. The weight function is calculated between two particles at a distance r, as shown in Eq.(1),which satisfies the condition that there is no selfinteraction in one particle.

Fig. 1 Interaction model of MPS method

The summation of all particle weight functions is defined as the particle number density,as shown in Eq. (2).The particle number density was used to evaluate the incompressible fluid density. The integration of the weight function with the particle number density can represent the flow of fluid in the interaction zone and the number of particles in the interaction zone fluctuates. However, the weight function balances it by the definition that w(r )becomes low as the distance of the particle approaches re.Moreover, w(r ) will have an asymptotic value near a zero distance, which prevents fluid particles from passing through other fluids or wall particles [16]. The MPS method has an expression for gradient, divergent and Laplacian models, as shown in Eqs. (3), (4) and (5),respectively.

All three models are used to discretize the differential operator in the MPS governing equation. These three models depend on the weight function. In Eq. (5) of the Laplacian model, the λ value can be approximated by Eq. (6) to be proportional to the analytical solution [9].

The fluid density is assumed to be constant as the incompressibility of the media is implemented,as shown in Eq. (7) [17]. The Navier–Stokes equation is used to calculate the movement of an incompressible fluid represented by the velocity change. The fluid particle moving mechanism is conducted by an external force, which is represented as the gravitational acceleration. Moreover, the viscosity determines the frictional force and momentum transfer between the particles when a collision occurs.Equation (9) shows that the energy conservation equation is composed of the conduction and heat source terms.Following the Laplacian form of the MPS method, Eq. (9)can be rewritten as

Moreover, the enthalpy is a function of temperature, as described in Eq. (11):

Each particle temperature is explicitly calculated based on Eq.(12).The calculation is performed after the particle positions are updated and it does not include the temporal value of each step. The external force and viscosity term are calculated simultaneously with the particle temperature calculation, as shown in Eq. (13). Meanwhile, in Equation (14), the pressure term is calculated implicitly. The temporal position is updated owing to the velocity change,as shown in Eqs. (15) and (16):

Figure 2 shows the flowchart of MPS calculations; it begins with the preparation of the initial conditions. All particles have their initial conditions, such as the position,velocity, pressure and temperature, and these data are stored in a grid data file. Other data that must be included in the initial condition are the average distance between particles, mass density, compressibility and kinematic viscosity. The value of gravity is 9.81 m/s2for a particle influenced by gravity, such as a fluid. After all preparations, the calculation starts with the initial particle number density.There is a temporary velocity that is influenced by the viscosity and gravity. This temporary velocity triggers the movement of some particles. Eventually, the configuration of the particle changes. In addition, the temperature and pressure are calculated for each particle using Eqs.(12)and (14), respectively. The simulation is terminated after the steady-state condition is reached.

2.2 Verification of heat transfer calculation

To examine the accuracy of the MPS method, the heat transfer was verified using a 2D slab model, as shown in Fig. 3a.The MPS calculation was performed assuming that no gravitational acceleration worked on the particles.Both materials are listed in Table 1; the material specifications are close to those of stainless steel.The analytical solutions of the transient temperature profile were obtained using Eqs. (17) and (18) [18]:

Table 1 Material specification for heat transfer model calculation

The thermal diffusion coefficient is shown in Eq. (19).The values T1,iand T2,iare the temperatures of the material at the end of the slab,as shown in Fig. 3.At the beginning of the simulation, the red and blue particles were 600 and 300 K, respectively. The red and blue particles are materials 1 and 2, respectively. Tsis the interface temperature calculated using Eq. (20). The interface temperature depended on the density, heat conductivity, specific heat coefficient and initial temperature of both slabs.

Fig. 2 Improved MPS flowchart

Fig. 3 Heat transfer model of two-dimensional MPS at a 0 s and b 4 s

Fig. 4 Comparison of analytic and MPS heat transfer calculation at 4 s

Fig. 5 Three-dimensional illustration of experiment

Figure 3b depicts the MPS simulation results of the twoslab heat transfer process at 4 s.The temperature tended to be constant at the edges of each slab. However, at the interface, there was a temperature change due to heat transfer from material 1 to material 2. Figure 4 shows a comparison between the MPS and the analytical solution.The discrepancy between the MPS and analytical calculations was miniscule. The discrepancy was 2 K, which was insignificant compared to the initial condition of the slabs with a temperature difference of 300 K. Therefore, this MPS method could be used for the next heat transfer calculation between two materials with a high temperature difference.

2.3 Analysis of experiment

Based on Sudha’s experiment, the initial condition was set in a stainless steel SS304L melt chamber consisting of a cylindrical and conical section [8]. Six ceramic beads and quartz-insulated tape heaters were used to heat the Wood’s metal. Moreover, an auxiliary pressure system could regulate the pressure of the material at 2 atm. The melt chamber system was capable of heating the metal up to 673 K. The chamber was thermally insulated to minimize ambient heat loss during the experiment. Seven K-type thermocouples were placed in various locations to control and monitor the Wood’s metal and vessel test temperatures.The thermocouples were insulated by SS-sheathed MgO.Two thermocouples were placed on the top and bottom surfaces of the first WMP.Then,the molten Wood’s metal was placed inside the melt chamber and preheated until it reached 573 K.It was released through a nozzle 60 mm in diameter from the melt chamber into a cylindrical test

vessel. The stainless steel SS304L cylindrical test vessel was 950 mm in diameter and 1100 mm in height. However, the report did not mention the length of the melt chamber pipe. Therefore, some assumptions were used in the simulation. Two WMPs were mounted on dummy SS inside the cylindrical test vessel. The diameter and thickness of the WMPs were 470 mm and 7 mm, respectively.The distance between the nozzle and the WMPs was 270 mm. The cylindrical test vessel was partially filled with water immediately below the first WMP. This scenario was analogous to a real condition inside a nuclear reactor in which the grid plate had a closed outer boundary.Figure 5 illustrates the experiment.The simulation focused on the red dashed line, which shows the downward relocation of the molten material to the plate.

The experiments were conducted under two scenarios.First,the two WMPs were placed inside the cylindrical test vessel without any barrier restricting the outflow of molten material and this scenario was used for the first to the fourth experiments.This caused the molten material to spill out directly into the water after impinging the WMP.In the second scenario in the fifth to eleventh experiments, there were perspex cylindrical rings to restrict the flow of all molten material on the plate. The sequence of the Wood’s metal impingement process was captured using a 1000 fps high-speed camera. After conducting the first four experiments,the holes that formed on the WMP on both surfaces were measured. The formed hole areas on the top and bottom surfaces were 170 cm2and 20 cm2, respectively.The diameter of each surface was approximately 9.44 cm for the top surface and 5.04 cm for the bottom surface.The formed hole diameter data were used to validate the hole formation in the simulation results.

2.4 Simulation of downward relocation process

In this study,the 2D MPS simulation was mainly used to analyze the heat transfer between the molten Wood’s metal and the first WMP. Then, we approximated the hole formation behavior when the molten Wood’s metal breached the WMP. Therefore, not all geometries were built in the simulation. It focused only on the first to the fourth experiment without the perspex cylindrical rings mounted above the WMP. It is simulated in an x-y Cartesian coordinate system, as shown in Fig. 6. The initial pressure of the molten Wood’s metal inside the melt chamber was set at a zero-bar gauge, as experimentally conducted. Moreover, in the original MPS method, a pressure of 0 Pa was applied at the beginning of the calculation to avoid numerical instability owing to negative pressures[15, 16, 18]. Originally, there were four types of particles in MPS, namely fluid, rigid body, wall and dummy wall particles. The fluid particles were used to represent the molten Wood’s metal.The rigid body particles represented a more viscous fluid than the previous particle.However,in this simulation, the rigid body particle was not used because it was necessary to present only one type of fluid particle. The dummy wall particle was a type of particle that behaved like a wall without the application of pressure and temperature calculations. It was used to give geometrical boundaries of the wall particle to the void condition.The wall particle was a type of solid particle in which the pressure and temperature calculations were performed, as well as fluid particles. Here, it was necessary to represent the WMP and SS304L as wall-type particles.Therefore,to distinguish between the two, a new stainless steel wall particle, namely a second wall, was declared and is represented by the magenta color in Fig. 6. The particles in red, purple and blue represent the molten Wood’s metal,WMP and dummy wall, respectively. The initial temperatures of the molten Wood’s metal,stainless steel and WMP were set to 573 K, 673 K and 300 K, respectively. The diameter of all the particles was 1 mm.The total number of particles in the simulation was 15,875. The ratio of the particle number density was set to be 2.1 times larger than the particle diameter, which means that the revalue was 2.1 mm.

Based on the experimental results,the phase change due to heat transfer occurred quickly.Therefore,the simulation performed the main critical time when the impingement and phase change of the WMP was 2 s. The time step was 0.001 s, so there were 2000 sequence results in one simulation.They can be represented as the sequences provided by the high-speed camera used in the experiment. However, the result sequences were presented with a time step of 0.1 s to minimize the memory utility.The calculation is assumed to be convergent when the temperature gradient in a limiting radius reis sufficiently small and may not conduct the heat transfer between particles, especially for the molten material and the metal plate. Table 2 lists the physical properties of the Wood’s metal and stainless steel SS304L used in the simulation.Subsequently,the obtained experimental temperature data of both main thermocouples placed on the top and bottom surfaces of the first WMP were used to validate the simulation results. Both were obtained by digitizing the reported data particularly for the time evolution of the temperature and the prediction time when the hole on the surface started forming. Some information was also provided, including the remaining solid fraction of the WMP and the formed hole due to the impingement of the molten Wood’s metal. Therefore, it is the validation data for the MPS simulation results.

3 Result and Discussion

Figure 7 depicts the simulation result of the downward relocation process at 0.4 s when the molten material impinged the WMP. It can be seen that there were color changes on some WMP surface particles, indicating heat transfer between the molten Wood’s metal and the plate.However, the hole was not observed at that time because the WMP surface particles remained in the original position. Along with the heat transfer between the two materials, the phase change would occur simultaneously with the velocity change, which causes the particle to move.After flowing down to the WMP,all of the molten Wood’s metal particles spread horizontally on the plate.The spread was counted as the leading-edge distance of the molten Wood’s metal particles.

The detailed sequences of the molten Wood’s metal impingement on the WPM are shown in Fig. 7. Impingement was observed with no temperature change on the plate during 0.2 s of the molten Wood’s metal relocation.The range of velocities of the molten Wood’s metal when it reached the WMP surface was 1.25–2.0 m/s.The breaching process to the bottom surface of the WMP was found at 1.2 s.Figure 8b shows the breaching process of molten Wood’s metal at 1.3 s when the melted WMP started to flow down owing to the impingement and the melting process. It also indicated that the molten Wood’s metal had spread horizontally on the top surface of the WMP. The molten Wood’s metal temperature at the WMP top surface decreased owing to the interaction with the solid particles.Moreover, the hole formations became more prominent at 1.4 s, as shown in Fig. 8c. The melted WMP particles gradually mixed with the molten particles. The particles spread on the top surface tended to move and then flow down through the formed hole,and then reach the stainless steel cylindrical test vessel.Therefore,there was a decrease in the local height of the molten Wood’s metal, which spread qualitatively.Moreover,because of the high melting point of the cylindrical test vessel and the low melting point of the Wood’s metal, it could be assumed that there was no heat transfer between the molten material and the vessel.

Fig. 6 Initial condition of MPS simulation

Figure 9 shows a comparison of the temperature profiles obtained from the simulation and experimental results.There was a significant change in the temperature of the WMP top surface at the beginning of the relocation of the molten Wood’s metal. The experimental data showed that the temperature surpassed the WMP melting point at 0.3 s.Meanwhile,in the simulation results,this occurred at 0.4 s,which indicates a phase change and the start of hole formation on the top surface of the WMP. Both experimental and simulation results showed that the maximum recorded temperature was 553 K with a maximum top side temperature at 0.5 s. Furthermore, there was a sharp temperature decrease in the experimental result, while the simulation results showed a gradual decrease in temperature. When the molten Wood’s metal vertically impinged the WMP,it led to the penetration of the bottom surface of the plate. Moreover, the cooling process of the molten Wood’s metal in air accelerated the temperature decrease on the WMP surface.

In the experiment, the temperature change on the top surface appeared to have a sharply decreasing slope, while the temperature at the bottom surface increased significantly. The same temperature trend was observed in the simulation results. However, at a specific time step, the temperature tended to have a gentle slope until the end of the simulation. This was caused by the cooling process,which depended on the particle interactions and was not simultaneously cooled by the real air condition. The bottom surface temperature also exhibited the same pattern,as shown in Fig.8.The results showed that the WMP melting point was surpassed at 0.6 s and 1.2 s in the experiment andsimulation, respectively. Experimentally, the molten Wood’s metal penetrated the plate in an interval time of 0.3 s.Meanwhile,the penetration time of molten Wood’s metal took 0.8 s in the simulation result. The penetration delay could be caused by the independence of the temperature calculation in the pressure module calculation. When the particles underwent a gradual temperature change over a certain time,the viscosity changed,which led to a pressure change around the particles at a certain effective radius.This change is significant when there is a small temperature gradient between the particles and small fluid thermal conductivity, but a large heat capacity is required to increase the temperature of the solid particles. In the simulation,the temperature gradient was high at the beginning and the penetration occurred after less than 1 s. Therefore,there may be some delay, but the results showed that a similar bottom side maximum temperature was reached at 412 K and 410 K for the experiment and simulation,respectively. The small discrepancy did not significantly impact hole formation during the impingement of molten Wood’s metal on the WMP.

Fig. 7 MPS result of Wood’s metal relocation to the WMP at 0.4 s

Fig. 8 Simulation sequence results for the WMP at a 0.4 s, b 1.3 s, c 1.4 s

Fig. 9 Temperature profile on the WMP surfaces

The phase change of WMP particles due to the heat transfer process first occurred at 0.4 s,as shown in Fig.10.The solid and liquid phases changed linearly and oppositely until the end of the simulation. The results indicated that the molten Wood’s metal phase increased, while that of the solid phase of the WMP decreased. The fraction change became stagnant after 1.6 s,along with the flowing down of molten Wood’s metal to the cylinder chamber.Generally,it is assumed that all molten materials reach the reactor catcher, which is located near the main reactor vessel. However, the experimental results showed that the fraction of molten material that reached the core catcher was only 0.720.Similarly,this simulation showed that only 0.750 of the molten material was formed until the end of the simulation time. Therefore, the simulation result has a 5.50% absolute error compared with the experimental results.

Fig. 10 Variation of WMP phase fraction

The hole size or any impingement trace measurement may be performed after the experiment is completed.However,in the simulation case,it could be predicted step by step. The erosion and heat transfer process simulation results in previous studies were also validated by approaching the results of the experiment [12, 19–21].However, the validation did not clarify the dynamic behavior of the hole-size increment during the impingement process.Therefore,this simulation result revealed the physical trend of diameter size enlargement during the impingement of molten Wood’s metal on the WMP. Figures 11 and 12 depict the diameter of the formed hole due to the melting process in the WMP. The top and bottom holes started to form at 0.4 s and 1.2 s, respectively. The discrete data of hole-size enlargement were approached by the fitting pattern, as shown in Eqs. (21) and (22).The coefficient values were 95 and 50 for the top and hole diameter enlargements, respectively. The calculation may result in a maximum of unity if there is no amplification factor from the coefficient. Therefore, the coefficient values were approximated based on the maximum diameter achieved in the experiment.The approximation was similar to the equation reported by Johnson, Mehl, Avrami and Kolmogorov (JMAK) for the transformation kinetics[22, 23]. In the JMAK equation, k and n characterize the reaction as a function of time, temperature and other variables that have significant influences on the transformation process. Initially, the JMAK equation was used to explain the solid-state transformation kinetics due to isothermal diffusion. However, recent studies have revealed that the transformation also worked under nonisothermal conditions [24–26].

The enlargement of the hole diameter due to the Wood’s metal impingement occurred when there was a high temperature change in a short time, approximately 2 s. The enlargement was also caused by the phase change of the WMP from the solid phase to the liquid phase.Data fitting with Eq.(21)was performed to determine the enlargement behavior on the top surface of the WMP. The calculation result showed a significant increase at 0.6 s, while the simulation result was not the same.However,the pattern of calculation and simulation results showed the same value from 0.7 to 0.9 s. The saturation result started from 1.1 s along with a significant increase in the result from Eq.(22).Physically,it was explained as an indication that the size of the top surface-formed hole significantly increased at the beginning of the simulation until 0.9 s, as shown in Figure 10. After that, it changed slightly owing to the horizontal particle movement so that the vertical penetration decreased. Thus, from 1 s to 1.6 s, the top diameter changed from 80 mm to 95 mm. As a result of vertical penetration, the WMP bottom surface would form a hole,although its size was not as large as that on the top surface,as shown in Fig. 12. The experiment result obtained previously showed that the formed hole diameters on the top and bottom surfaces were 94.43 mm and 50.48 mm,respectively.Likewise,the simulation results indicated that the diameters of the top and bottom sides were 95 mm and 50 mm, respectively. It can be concluded that the absolute errors for the top and bottom formed hole diameters were 1.2% and 1.8 %, respectively.

Fig. 12 Variation in diameter of the WMP’s bottom surface-formed hole

Fig. 13 Molten Wood’s metal leading-edge position

Figure 13 shows the molten Wood’s metal outermost particle distance when it impinged on the WMP plate.This distance was calculated from particles distributed radially to the x-positive and x-negative. The previous research reported that the molten particles spread quickly above the WMP.The MPS simulation showed that at 1.1 s,all molten particles had spread evenly until the edge of the WMP.After that, some molten particles fell into the cylindrical test vessel.At the first time that molten particles impinged the WMP, their maximum velocity was 1.25 m/s. However,when all particles had spread to the edge of the WMP,the maximum velocity decreased to 0.22 m/s.The decrease was caused by a change of momentum and transfer of energy from the molten particles to WMP particles.

4 Conclusion

The improved MPS method has been used to analyze the molten metal behavior during relocation on the plate. The results show that during the spread of molten Wood’s metal, some energy or heat was also transferred to the WMP so there was a phase change from a solid to liquid phase. From 0.4 s to 1.6 s, there was a linear volume fraction change due to the solid–liquid phase change. At the end of the simulation, the total volume fraction of the formed molten WMP was 0.750 and it had a 5.50%absolute error compared to the experiment. Another consequence of heat transfer between the molten Wood’s metal and WMP was hole formation in the center of the WMP. At the top and bottom surfaces of the WMP, holes with respective diameters of 95 mm and 50 mm were formed.The absolute errors for the top and bottom formed hole diameters were 1.2% and 1.8%, respectively. The temperature profile was observed in detail in the simulation results and showed that there was a decrease in the molten Wood’s metal temperature after breaching the plate.All the simulations showed good agreement and provided an optimal prediction of the behavior of the molten material during the relocation and impingement processes, which was not explained by the experiment. Some empirical formulas have been introduced to explain the hole formation mechanism of the metal plate.For some improvement,a 3D simulation of the downward relocation process of Wood’s metal can be applied and compared with the 2D simulation results. Moreover, an improved MPS can be implemented to analyze the molten material behavior inside a nuclear reactor when a severe accident occurs.

AcknowledgementsThe authors wish to thank Prof. S. Koshizuka,Prof.M.Sakai and Dr.K.Shibata of the University of Tokyo for their helpful comments and for providing the basic MPS code for fluids.

Author ContributionsAll authors contributed to the study conception and design. Anni Nuril Hidayati performed the numerical calculation, processed the experimental data, performed the analysis,drafted the manuscript and designed the figures. Abdul Waris, Asril Pramutadi Andi Mustari and Dwi Irwanto were involved in planning and supervised the work. Asril Pramutadi Andi Mustari and Nur Asiah Aprianti aided in interpreting the results and analysis.The first draft of the manuscript was written by Anni Nuril Hidayati and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.