Analysis of molten metal spreading and solidification behaviors utilizing moving particle full-implicit method

Abstract

To retrieve the fuel debris in Fukushima Daiichi Nuclear Power Plants (1F), it is essential to infer the fuel debris distribution. In particular, the molten metal spreading behavior is one of the vital phenomena in nuclear severe accidents because it determines the initial condition for further accident scenarios such as molten core concrete interaction (MCCI). In this study, the fundamental molten metal spreading experiments were performed with different outlet diameters and sample amounts to investigate the effect of the outlet for spreading-solidification behavior. In the numerical analysis, the moving particle full-implicit method (MPFI), which is one of the particle methods, was applied to simulate the spreading experiments. In the MPFI framework, the melting-solidification model including heat transfer, radiation heat loss, phase change, and solid fraction-dependent viscosity was developed and implemented. In addition, the difference in the spreading and solidification behavior due to the outlet diameters was reproduced in the calculation. The simulation results reveal the detailed solidification procedure during the molten metal spreading. It is found that the viscosity change and the solid fraction change during the spreading are key factors for the free surface condition and solidified materials. Overall, it is suggested that the MPFI method has the potential to simulate the actual nuclear melt-down phenomena in the future.

Key words： molten metal spreading; solidification; particle method; severe accident; fuel debris; decommissioning

1 Introduction

The accident of Fukushima Daiichi Nuclear Power Plants (1F) has highlighted the importance of further understanding of nuclear severe accidents. At present, there are a lot of reports from utilities on fuel debris in 1F. According to the reports with respect to the observation inside primary containment vessel (PCV) [1,2], the part of the fuel was found in Unit 2 while a part of the control rod was found in Unit 3. Furthermore, the accident progression analysis using SAMPSON [3] indicated that the shell attack, which is the molten corium invasion against the PCV wall, might happen in Unit 1. These reports suggest that different types of fuel debris might have generated in each unit of 1F. Therefore, the transient phenomena in the generation process must be studied to understand the existing fuel debris in 1F. The molten corium spreading on the PCV bottom is one of the important phenomena because it determines the initial distribution and composition of the fuel debris in the PCV. The initial condition may affect the progression of the molten core concrete interaction (MCCI), which is a long-term interaction between the corium and containment floor. Therefore, the spreading behavior is one of the keys for 1F decommissioning. However, the corium spreading is quite complicated because corium would consist of the solid-state, the liquid state, and the solid-liquid mixture state, and the spreading and solidification behavior affect each other. To deeply understand such complicated phenomena, both experiment and simulation approaches are essential to accumulate the knowledge of spreading phenomena.

Several spreading experiments using simulant materials and prototypical corium have been conducted by various organizations [4–14]. The majority of these experiments were mainly focused on the molten metal spreading behavior for core catcher design [4,15]. The main concept of the core catcher system is to prevent further accident progression by receiving the corium and removing heat from the corium. Thus, in these experiments, simulant materials were injected into the surface of the receiving plate from the bottom or side without considering downward jet flow and collision between molten metal and the receiving plate. The experiments for considering the collision effects were conducted after Fukushima accident [16–19]. Ogura et al. [16] conducted the experiments using a high-speed thermo-camera to investigate the spreading behavior for different scale effects represented by outlet diameters and falling height. In their experiments, molten copper was dropped on the stainless-steel surface, the spreading behavior was observed and recorded via high-speed thermal camera, and a dimensionless analysis was performed to determine the main parameters for molten metal spreading behavior such as solidified area and solidified thickness. They found that the outlet diameters and the falling height were crucial parameters for the melt spreading behavior.

The authors of the present paper conducted melt spreading experiments with various outlet diameters and molten metal amounts [20]. In 1F accident, it is indicated that the degradation condition at the bottom of the reactor pressure vessels (RPV) in different units might be different [21]. The understanding of the dependence of melt spreads on the outlet diameters will help understand the fuel debris spreading in 1F.

The experimental works [16–19] also provided some empirical equations to predict the edge portion and the thickness of molten metal. Although they could empirically approximate the spreading and solidification, the mechanism of the phenomena has not been clearly understood. Numerical simulation is useful in predicting the phenomena based on mechanical understanding.

Some spreading and solidification models [22,23] and spreading codes, such as MELTSPREAD [23,24], CORFLOW [25], THEMA [26], and LAVA [27], have been developed and validated against the experiment results. However, these codes include the empirical parameters to be determined in obtaining reasonable results. In addition, N-S equations in these codes are integrated into the thickness direction in order to use the height function to represent the free surface, which cannot express vertical distribution of the temperature and composition. Therefore, the method to calculate multi-dimensional free surface flow is desired for capturing the behavior of the spreading in more detail.

The moving particle semi-implicit method (MPS) is one of the particle methods to calculate the motion of the incompressible fluid developed by Koshizuka and Oka [28]. The melting and solidification model is also developed in the MPS framework [29–34]. However, the MPS method cannot calculate the fluid which does not have a free surface because the boundary condition is set by determining the surface of the fluid. In addition, it cannot handle the behavior of the solid properly because the angular momentum is not conserved. Moreover, it cannot treat the problem where energy conservation is important because the kinematic energy in MPS is unstable.

The moving particle full-implicit method (MPFI) [35,36] is also one of the particle methods for incompressible fluid, which can satisfy the above requirements. In comparison with the MPS method, MPFI has several advantages in the analysis of fluids with melting and solidification.

In this study, fundamental experiments on molten metal spreading and solidification with various outlets [20] were conducted using the MPFI method. It is expected that the understanding of the effect of the outlet for the melt spreading and solidification behavior facilitate the understanding of the spreading and solidification of fuel debris in 1F. In the experiment, a melted U-alloy was used as a molten corium simulant, and its spreading on a stainless-steel plate was recorded with a video camera at 60 fps at varying nozzle sizes and sample amounts. The area and thickness of the solidified materials were measured by image analysis. In the 3D simulation, the comparison of experiment and simulation results on the solidified areas history and the temperature profile were discussed in both quantitative and qualitative to further understand the solidification behavior during the molten metal spreading. Furthermore, to investigate the free surface condition and the thickness with the different outlet diameters and the sample amounts, several 2D simulations were performed.

2 Experimental setup

In this experiment, U-alloy 70 which has a melting point of 70°C was introduced as the simulant material of molten corium in terms of similarity of its density with UO₂. The density is one of the important parameters in spreading because it affects the dimensionless parameters such as Reynolds number, Weber number, and Froud number [20]. The test facility and test conditions are described in Fig. 1 and Table 1, respectively. The test section consists of the vessel-heating part and spreading part. The vessel was heated by the ribbon heater up to 200°C, which was controlled by the temperature controller. To measure the temperature of the U-alloy 70, thermocouples were welded to the center of the materials, and the temperature history was recorded via a data logger. When the temperature of the sample reached 170°C, the plug was pulled out and the liquid metal flew out from the vessel and spread on the stainless-steel plate. In the spreading section, the temperature of the stainless-steel-plate was controlled using another temperature controller and heater. To adjust the falling height, a height adjuster was installed below the heater. The spreading behavior was recorded using video cameras at 60 fps (HANDYCAM HDR-CX485, SONY). To measure the temperature history of the molten metal in the spreading section, 2 thermocouples were set on the flat plate, one configurated just under the outlet center (Fig. 1, TC1) while the other set 2 cm apart from the center (Fig. 2, TC2). After complete solidification of deposited metal, the spreading area of the solidified material was estimated via image analysis (Image J), and the thickness of the solidified material was measured using the caliper. The spreading experiments were performed 3 times at least to confirm the experimental reproducibility. The detailed experiments results were described in Ref. [35].

Fig.1 Schematic image of experiment facility.

Full size|PPT slide

Tab.1 Experimental values for spreading experiment

Parameters	Value
Molten metal mass/g	100, 200
Initial temperature of molten metal/°C	170
Outlet nozzle diameter/mm	3, 5, 7, 9
Fall height/mm	50
Stainless-steel temperature/°C	25
Density of molten metal/(\mathrm{kg}\cdot {\mathrm{m}}^{-\mathrm{3}})	9580

3 Numerical method

3.1 MPFI method

For the numerical simulation, the MPFI method [35,36], which is thermodynamically consistent, was applied. In the MPFI method, the following Navier-Stokes (N-S) equation was adopted [35,36] to simulate the incompressible free-surface fluid flow, which is expressed as

(1)\rho {\displaystyle \frac{\mathrm{d}{u}_i}{\mathrm{d}t}}={\displaystyle \frac{\partial }{\partial x_i}}\mu {\dot{\epsilon }}_{ij}+{\displaystyle \frac{\partial }{\partial x_i}}{P}_i+\rho g_i,

where g, P, u, ρ, μ, and ε are gravity acceleration, pressure, velocity, fluid density, fluid viscosity, and strain rate, respectively. Besides, the subscripts i and j denote the particle i and particle j, respectively. The first term on the right-hand side is the viscosity term and the second one is the pressure term. Here, the pressure term is expressed as

(2)P=\left(-\lambda \dot{\epsilon }+\kappa {\displaystyle \frac{\rho -{\rho }_{\mathrm{0}}}{{\rho }_{\mathrm{0}}}}\right),

where λ and κ are bulk viscosity and bulk modules, respectively. The governing equations (Eqs. (1) and (2)) will approach the incompressible N-S equation in the case where the parameter λ and κ are set large enough. In the particle methods [28,37], the governing equations are replaced by particle interaction forces. The particle interaction is limited in a finite range using an effective radius and a weight function. In the MPFI method, the particle interaction models for gradient, divergence, and Laplacian operators are formulated as

(3)\nabla \varphi ={\displaystyle {\sum }_j({\varphi }^j+{\varphi }^i)}{r}_{\square }^{ij}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}},

(4)\nabla ·A={\displaystyle {\sum }_j}(A^j-A^i)r^{ij}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}},

(5){\nabla }^{\mathrm{2}}\varphi ={\displaystyle {\sum }_j}({\varphi }^j-{\varphi }^i){\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}},

where φ, A, r^ij, and d^ij are an arbitrary scalar, an arbitrary vector, the relative position between particles i and j, and the particle distance, respectively; and {w^{\prime }}^{ij} is a differential of the weight function defined as

(6)

where r_e is the effective radius, and the subscript e denotes the initial of the effective. The particle interaction models are formulated similarly in comparison with the smoothed particle hydrodynamics (SPH) [37]. However, the weight function of Eq. (6) instead of the SPH kernel function is used in the MPFI method.

When Eq. (1) is discretized with the particle interaction models, the force acting on the particle is formulated as

(7)\rho {\displaystyle \frac{\mathrm{d}{u}_i}{\mathrm{d}t}}=\mathrm{2}\mu {\displaystyle {\sum }_j}(u^j-u^i){\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}-\lambda {\displaystyle {\sum }_j}\left({\dot{\epsilon }}^j+{\dot{\epsilon }}^i\right)r^{ij}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}+\rho g_i,

where

(8)P^i=-\lambda {\displaystyle \sum _j(u^j-u^i)}{r}^{ij}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}+\kappa {\displaystyle \frac{n_i-n_{\mathrm{0}}}{n_{\mathrm{0}}}}

(9)n_i={\displaystyle {\sum }_j}{w}^{ij},

where w^ij and n_i are the weight function and particle number density, which is a summarization of the weight function. Since the potential energy and dispersive function for the discretized equation can be described [36], the thermodynamic consistency of the particle system after discretization is assured. This implies that the mechanical energy of the system monotonically decreases, and this property is crucial for stable calculation. However, simply applying the Laplacian model to the viscosity term does not allow the high-viscosity material to be rotated because the pairwise force against the relative velocity of two particles yields torque against rotational motion. To allow the rotation, the Laplacian model was modified by introducing angular velocity and subtracting rotational element from relative velocity. Using the modified Laplacian model, the discretized equation was re-formulated as

(10)\rho {\displaystyle \frac{\mathrm{d}{u}_i}{\mathrm{d}t}}=\mathrm{2}\mu {\displaystyle \sum _j\left(u^j-u^i-{\displaystyle \frac{{\omega }^j+{\omega }^i}{\mathrm{2}}}{r}^{ij}\epsilon \right)}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}-\lambda {\displaystyle \sum _j\left({\dot{\epsilon }}^j+{\dot{\epsilon }}^i\right)}{r}^{ij}{\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}+\rho g,

where the angular velocity of the particle is calculated from the angular momentum conservation equation as

(11){\displaystyle {\sum }_j}(u^j-u^i-{\omega }^i{r}^{ij}{\epsilon }_i)r^{ij}{\epsilon }_j{\displaystyle \frac{w^{ij}}{d^{ij}}}\mathrm{\Delta }t=\mathrm{0.}

Equations (8), (10), and (11) will be a linear matrix equation having unknown values of velocity u, angular velocity ω, and pressure P. Since its coefficient matrix is symmetric, it can be solved by utilizing the conjugated residual method.

3.2 Heat transfer model

To simulate the melting-solidification, heat transfer and temperature calculations are essential. The temperature field is calculated with the energy conservation equation

(12){\displaystyle \frac{Dh}{Dt}}=k{\nabla }^{\mathrm{2}}T+Q,

where h, k, T, and Q are enthalpy, thermal conductivity, temperature, and heat input, respectively. The equation is discretized by using the particle interaction model for the Laplacian operator in MPFI (Eq. (5)) as

(13){\displaystyle \frac{D{h}_i}{Dt}}={\displaystyle \sum _j{k}_{ij}}(T^j-T^i){\displaystyle \frac{{w^{\prime }}^{ij}}{d^{ij}}}+Q.

The thermal conductivity k between two particles i and j are evaluated as the harmonic average of the two particles [29,32]. The modified energy conservation equation is expressed as

(14){\displaystyle \frac{D{h}_i}{Dt}}={\displaystyle \sum _j\frac{\mathrm{2}{k}_i{k}_j}{k_i+k_j}}(T^j-T^i){\displaystyle \frac{{w^{\prime }}^{ij}}{d_{ij}}}+Q.

Regarding the heat transfer problem between the solid phase and liquid phase with flow, it is considered that the heat transfer will be accelerated due to the flow. In this simulation, the Nusselt number under the flat surface condition is introduced to estimate the heat transfer as

(15)Nu=\mathrm{0.664}R{e}^{\mathrm{1}/\mathrm{2}}P{r}^{\mathrm{1}/\mathrm{3}}\left(Re\hspace{0.17em}\le \hspace{0.17em}\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{5}}\right),

(16)Nu=\mathrm{0.037}R{e}^{\mathrm{0.8}}P{r}^{\mathrm{1}/\mathrm{3}}\left(Re\hspace{0.17em}>\hspace{0.17em}\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{5}}\right).

In addition, the heat loss due to the radiation from the free surface is considered. In the MPFI method, the radiation heat transfer is simply modeled as the heat removal from free surface particles based on Stephan-Boltzmann’s law,

(17)q_{\mathrm{radiation}}=\sigma \eta A\left(T^{\mathrm{4}}-T_{\mathrm{e}}^{\mathrm{4}}\right),

where σ, ε, A, and T_e denote the Stephan-Boltzmann constant, emissivity, surface area of the particle, and environment temperature, respectively. The radiation heat transfer was applied only on the free surface particles. In this study, the free surface detection was conducted based on the number of neighbor particles as Farmer et al. [24] adopted in the MPS method [28].

3.3 Phase change and viscosity model

Modeling the phase change is a key issue of simulating the spreading-solidification behaviors of molten metal. In the MPFI method, phase changes are expressed by changing the particle viscosity according to their enthalpy and solid fraction. The temperature which is calculated as a function of enthalpy and a solid fraction (γ) of the particle is calculated as

(18)

(19)

where C_ps, C_pm, C_pl, h₀, h₁, and T_m are specific heat for solid state, specific heat at the melting point, specific heat for liquid, solidifying enthalpy, liquefying enthalpy, and melting point, respectively. When γ = 0 , the particle is in a complete fluid state; when γ = 1, the particle is in a complete solid-state; when 0<γ <1 , the particle is a mixture of sold-liquid state.

Evaluation of viscosity depending upon the solid fraction is important. Previous investigations indicate that if the solid fraction is over 50%, the molten metal stops its spreading process due to the drastic increase of the viscosity [38,39]. On the nuclear field, analyses of melt spreading experiments for stainless steel SPREAD suggest that 55% of the solid fraction is needed to stop the spreading process. Moreover, it is suggested that the solid fraction depends on the change of the mass fraction, shear stress, and cooling ratio [38,39]. If the inertia force derived from increasing the nozzle diameter increases, the solid fraction to stop the further spreading also changes. Thus, in the MPFI, an empirical coefficient of Ramacciottii’s correlation [40], which depends on the solid fraction is used depending on the solid fraction.

(20)

(21)

where μ and C are viscosity coefficient and a constant value which is experimentally determined. In general, the C value is set from 4 to 8. In this study, the C value is determined as 6.4 based on Ref. [32].

3.4 Artificial cohesive force model

One of the significant issues in particle methods is that some of the particles are scattered when they move around. In this study, to prevent particle scattering, an artificial cohesive force was introduced based on the Many-body potential [41]. According to the report from Kondo [41], this artificial cohesive force has an effect of the surface tension. In the method, the parameter “a” having the dimension of bulk viscosity controls the magnitude of the force. The detail is shown in Ref. [41].

Table 2 is the summary of the boundary conditions and physical parameters in both 2D and 3D simulation. The simulation geometries in both 2D and 3D are demonstrated in Fig. 2.

Fig.2 Simulation geometry in 2D and 3D.

Full size|PPT slide

Tab.2 Calculation condition of MPFI simulation

Condition	2D	3D
Particle spacing/mm	0.004	0.009
Impact radius/mm	0.009	0.017
A	2.5e–7	2.5e–7
Sample amount/g	100, 200	100, 200
Outlet diameters/mm	3, 5, 9	3, 5, 9
Cps, Cpl/(J∙(kg∙K)–1)	200	200
Melting point/°C	70	70
Thermal conductivity of U-alloy70/(W∙(m∙K)–1)	36	36
Thermal conductivity of stainless steel/(W∙(m∙K)–1)	30	30
Latent heat/(J∙kg–1)	25000	25000
Density/(kg∙m–3)	9580	9580
Molten metal viscosity at 170°C/(Pa∙s)	3.3e–4	3.3e–4
Solid viscosity/(Pa∙s)	1e4	1e4
Stephan-Boltzmann constant	5.67e–8	5.67e–8
Emissivity	0.15	0.15
Surrounding temperature/°C	30	30

4 Results and discussion

4.1 Visualization of molten metal spreading behavior at different nozzle diameters

The snapshots of the molten metal spreading behavior at different nozzle outlets are described in Fig. 3. In the case of an outlet of 9 mm in Fig. 3(a), the molten metal rapidly spread and solidified at the beginning, hindered the further spreading of the melt, and formed a pool in the center region. Consequently, the solidified metal was rough and thick in the edge region and was smooth and thin in the center region. On the other hand, in the 3 mm outlet case (Fig. 3(b)), the molten metal gradually spread, keeping the thickness of the melt, and only the thick smooth area appeared after solidification.

Fig.3 Snapshots of spreading at different outlets.

Full size|PPT slide

4.2 MPFI simulation results

4.2.1 3D calculation at an outlet of 9 mm and a sample amount of100 g

In this section, the 9 mm outlet case is presented. Figure 4 is the snapshots of simulation results compared with the experiment [20]. Based on Figs. 4(b) and 4(c), the solidification procedure is summarized as follows: ① Up to around 1.0 s, the molten metal rapidly spread on the flat surface. During spreading and solidification, no noticeable change in solid fraction and viscosity is observed. ② Around 1.5 s to 2.0 s, the circle shape is formed. The solid fraction distribution (Fig. 6(b)) and viscosity distribution (Fig. 4c)) indicate that solidification begins at the edge of the spread melt, stops the further spreading, and forms a molten metal pool at the center. Furthermore, the drastic viscosity change in the vicinity of the interface is confirmed from Fig. 4(c). Due to this effect, further spreading is terminated. ③ Up to complete solidification, the solidification expands from the edge to the central region rapidly.

Besides, the quantitative comparison between the experiment and the simulation at a 9 mm outlet is described on the growth of solidified areas and temperature history in Figs. 5 and 6, respectively. As exhibited in Fig. 5, the growth of the solidified areas drastically increases up to around 0.5 s. The simulation result is in good agreement with the experiment. Meanwhile, the temperature changes between the molten metal and the receiving plate is measured at 2 measurement points. Figure 6 is a comparison of the experiment and the calculation. The simulation results basically follow the same tendency as the experiment. However, an overestimation of the temperature decrease is observed compared to the experiment. After drastic temperature elevation, the temperature almost remains constant in the experiment while a temperature decrease is calculated in the simulation. One possible reason is the setting of the heat transfer coefficient. In this model, the heat transfer coefficient is empirically determined via the Nusselt number, and it may include deviation. The other reason is the overestimation of radiation heat loss. Since the sensitivity with respect to the emissivity and the environmental temperature is large, their deviation may largely affect the simulation results. Therefore, the error in the heat transfer coefficient, emissivity, and environmental temperature could be a reason for the overestimation of the heat loss.

Fig.4 Snapshots of spreading at an outlet of 9 mm and a sample amount of 100 g.

Full size|PPT slide

Fig.5 Comparison of averaged solidified material of experiment and MPFI at an outlet of 9 mm and a sample amount of 100 g.

Full size|PPT slide

Fig.6 Comparison of temperature history of experiment and simulation at an outlet of 9 mm and a sample amount of100 g.

Full size|PPT slide

4.2.2 3D calculation at an outlet of 3 mm and a sample amount of 100 g

The qualitative comparison of the experiment and numerical simulation at an outlet of 3 mm is depicted in Fig. 7. As in the experiment, the molten metal gradually spreads on the flat surface. The spreading-solidification behavior is summarized as follows: ① The molten metal spreads gradually up to 3.0 s with the areas increasing. ② Around 3.0 s, a molten metal pool, which is clearly recognized by the red particle in Fig. 7 (b), is formed in the center region. At the same time, the melt begins to solidify, and the spreading almost stops. ③ The thickness of the molten metal gradually increases while the solidified area remains almost the same size. ④ Even when the spreading of the molten melt completely stops, the melt pool still exists. Eventually, it gradually cools down.

Comparing the cases at an outlet of 9 mm (Fig. 4) with that at an outlet of 3 mm (Fig. 7), the solidification behavior is found to be obviously different. A larger outlet results in a rapid solidification, while a gradual solidification is observed at the smaller outlet. The possible reason for this difference is the heat removal rate due to the interfacial area in contact with the flat plate surface.

Fig.7 Snapshots of spreading snapshots at an outlet of 3 mm and a sample amount of 100 g.

Full size|PPT slide

A quantitative comparison of the experiment and the simulation at an outlet of 3 mm is presented with respect to the solidified areas and temperature history in Figs. 8 and 9, respectively. The calculation results basically showed a similar trend to the experimental results. However, it is observed that the solidified areas in the simulation are slightly smaller than those in the experiment (Fig. 8). The possible reason is the setting of the heat transfer coefficient. Since the heat transfer coefficient used in this model is calculated based on the Nusselt number on the flat plate, some error may be contained in the heat transfer coefficient when the pool is formed. Figure 9 is a comparison of the experiment and the calculation on the temperature profile. The simulation showed a similar trend to the experiment. As in the case at an outlet of 9 mm (Fig. 6), an overestimation of heat removal is observed.

Fig.8 Comparison of averaged solidified material between experiment and MPFI at an outlet of 3 mm and a sample amount of 100 g.

Full size|PPT slide

Fig.9 Comparison of temperature history of experiment and simulation at an outlet of 3 mm and a sample amount of 100 g.

Full size|PPT slide

4.2.3 3D calculation at an outlet of 5 mm and a sample amount of 200 g

In Sections 4.2.1 and 4.2.2, the outlet effect on the spreading and solidification behavior was analyzed. In this section, the effect of sample amount on the melt spreading and solidification behavior is analyzed. The snapshots of melt spreading of both the experiment and numerical simulation are depicted in Fig. 10. Compared to the spreading behavior in Fig. 4 (an outlet of 9 mm and a sample amount of 100 g). It is observed that the molten metal does not spread rapidly even though the amount of the sample in Fig. 10 is twice as much as that in Fig. 4. The spreading behavior is the difference from both the case of the outlet of 9 mm and the case of the outlet of 3 mm. The outlet scales rather than the amount of the sample play a significant role in controlling the melt spreading behavior. The following molten metal spreading behavior is observed: The molten metal spreads up to 1.5 s while keeping its liquid state. In addition, from 1.5 s to 2.0 s, the further melt spreading is suspended due to the gradual solidification of the melt from the edge to the center. Moreover, from the 2.5 s to 4.0 s, although the melt discharged downward was terminated, the melt pool still existed. Furthermore, the molten pool areas gradually shrink.

The quantitative comparisons in both the history of the solidified areas and the temperature profile are depicted in Figs. 11 and 12, respectively. The simulation results were basically in good agreement with the experiments.

Fig.10 Snapshots of spreading at an outlet of 5mm outlet and a sample amount of 200 g.

Full size|PPT slide

Fig.11 Comparison of averaged solidified material of experiment and MPFI at an outlet of 5 mm and a sample amount of 200 g.

Full size|PPT slide

Fig.12 Comparison of temperature history of experiment and simulation at an outlet of 5 mm outlet and a sample amount of 200 g.

Full size|PPT slide

4.2.4 Free surface formation during spreading and solidification in 2D calculation

In the 3D calculation, the simulation results make it possible to understand the influence of the outlet scale for the melt spreading and solidification behavior. However, it is difficult to analyze the free surface condition and the thickness depending on the outlet scales from the 3D calculation because of the limitation of particle sizes. Therefore, in this study, to investigate the influence of the outlet and the amount of the sample on the free surface condition and the thickness of the solidified metal, the 2D simulation results are addressed.

The experimental results in Fig. 3 demonstrate that the size of the outlet drastically affects the free surface condition after solidification. In the case of the larger outlet case illustrated in Fig. 3(a), the rough free surface is formed close to the edge of the solidified material. In contrast, only the smooth free surface remains at the smaller outlet (Fig. 3(b)).

The calculation results at outlet of 9 mm are shown in Fig. 13, where the solid fraction and viscosity are presented in Figs. 13(a) and 13(b), respectively. The spreading material at the edge keeps the liquid state up to around 0.1 s. However, once the solid fraction value reaches 0.3, which is a critical solid fraction, the further spreading is terminated forcibly due to the existence of the solid particle and a rough surface is observed after solidification near the edge. These behaviors are equivalent to those observed in the experiment (Fig. 3). From the viscosity distribution in Fig. 13(b), it is found that the viscosity at the edge position drastically increases from 10^–4 Pa·s to 10⁴ Pa·s which is large enough to terminate the spreading. With this manner of solidification, the rough surface condition is formed close to the edge.

Fig.13 Snapshots of spreading behavior at an outlet of 9 mm.

Full size|PPT slide

The solidification behavior at an outlet of 3 mm is depicted in Fig. 14. In contrast to the case at an outlet of 9 mm (Fig. 4), the smooth free surface is formed. Up to around 1.4 s, the thickness of the spread material increases as the melt flows down, and a solidified material with the smooth free surface is produced. This tendency is in good agreement with the experiment in Fig. 3(b).

These results via the 2D simulation indicate that the outlet size may affect the formation of the free surface.

Fig.14 Snapshots of spreading behavior at an outlet of 3 mm.

Full size|PPT slide

4.2.5 Thickness comparison of different outlet diameters and samples

The solidified metal at different outlet diameters and the sample amounts is shown in Fig. 15. The thickness tends to increase as the nozzle diameter decreases regardless of the sample amounts. This tendency is observed from the experiment results as well. Figure 16 shows the thickness distribution in the horizontal direction obtained via the calculations. In the case at the outlets of 3 mm and 5 mm, the distribution is almost flat. However, in the case at an outlet of 9 mm, the thickness in the central region is thin while that in the edge region is thick. This tendency with respect to the thickness of the solidified metal is also observed in the experiment.

Fig.15 Thickness of solidified metals at different outlet scales.

Full size|PPT slide

Fig.16 Thickness distribution in the horizontal axis.

Full size|PPT slide

5 Conclusions

In this study, fundamental molten metal spreading and solidification behavior was experimentally and numerically investigated for further understanding of the fuel debris distribution which is essential to retrieve it from 1F. In the experiment, the spreading experiments considering downward jet flow and collision with a receiving plate were performed at different outlet diameters and sample amounts. During the experiments, the spreading behavior was recorded vis 60 fps video camera to analyze the spreading progression of the molten metal using image analysis. In the numerical investigation, the MPFI method which is one of the particle methods, was applied. The melting and solidification model including energy conservation, phase change, and temperature-dependent viscosity model was developed in the MPFI framework, and the spreading experiment was simulated with the model. In the case at an outlet of 9 mm, the molten metal spread rapidly and thinly solidified in a large area, while in the case at an outlet of 3 mm, the metal gradually spread, and solidified material remained thicker in a smaller area. From the simulation results, the order of the solidification during the spreading could be understood. In the case at an outlet of 9 mm, the molten metal drastically spread in the initial stage of the entire spreading behavior, and the molten metal pool was formed around the central region because the solidification started from the edge of the spread metal and blocked the further spreading, which later resulted in a rough solidification surface at the edge. In the case at an outlet of 3 mm, the molten metal gradually spread on the receiving plate with keeping the liquid state on the free surface. Consequently, the molten metal was solidified with a smooth free surface in a relatively small area. Besides, from the simulation results in the case at an outlet of 5 mm and a sample amount of 200 g, the outlet scale played a significant role rather than the sample amount. The 2D simulation was performed to reveal the influence of the outlet and the sample amount on the free surface and the thickness. The free surface formation was different depending on the outlet diameters. It was found that the thickness was almost flat in the case at an outlet of 3 mm and 5 mm regardless of the sample amounts. In contrast, it was found that the thickness was thin around the center while it was thick at the edge position in the case at an outlet of 9 mm. These behaviors were also observed from the experiments. Both the qualitative and the quantitative analysis indicated that the melting and solidification model had the potential to simulate the practical meltdown and spreading phenomena in Fukushima Daiichi Nuclear Power Plant (1F).

Acknowledgment

This study was partially supported by the Mitsubishi Heavy Industry (MHI).

References

1 TEPCO. Unit 3 Primary Containment Vessel Internal Investigation. 2017, available at the website of tepco.co.jp

2 TEPCO. Fukushima Daiichi Nuclear Power Station Unit 2 primary containment vessel internal investigation results. 2018, available at the website of tepco.co.jp

3 Pellegrini M, Dolganov K, Herranz E, Benchmark study of the accident at the Fukushima Daiichi NPS: best-estimate case comparison. Nuclear Technology, 2016, 196(2): 198–210

4 Sehgal B R. Nuclear Safety in Light Water Reactors-severe Accident Phenomenology. USA: Academic Press, 2012

5 Veteau J M, Wittmaack R. CORINE experiments and theoretical modelling. In: Van Goetem G, Balz W, Della Loggia E, eds. FISA 95 EU Research on Severe Accidents, Office Official Publ. Luxembourg: Europ. Communities, 1996, 271–285

6 Green G A, Finrock C, Klages J, Experimental studies on melt spreading, bubbling heat transfer, and coolant layer boiling. In: Proceedings of 16th Water Reactor Safety Meeting, Gaithersburg, MD, USA, 1988

7 Dinh T N, Konovalikhin M J, Sehgal B R. Core melt spreading on a reactor containment floor. Progress in Nuclear Energy, 2000, 36(4): 405–468

8 Suzuki H, Matsumoto T, Sakaki I, Fundamental experiment and analysis for melt spreading on concrete floor. In: Proceedings 2nd ASME/JSME International Conference 1, San Francisco, USA, 1993, 403–407

9 Eppinger B, Fieg G, Schuetz W, Stegmaier U. KATS experiments to simulate corium spreading in the EPR core catcher concept. In: Proceedings of the 9th International Conference on Nuclear Engineering (ICONE-9), Nice Acropolis, France, 2001, 32068804

10 Alsmeyer H, Cron T, Messemer G, ECOKATS-2: a large scale experiment on melt spreading and subsequent cooling by top flooding. In: Proceedings ICAPP'04 (International Congress Advances in Nuclear Power Plants), Pittsburg, USA, 2004

11 Alsmeyer H, Crin T, Foit J J, Test report of the melt spreading tests ECOKATS-V1 and ECOKATS-1. Constract FIKS – CT 1999 – 003. Ex-vessel core melt stabilization research. Karlsruhe: Wissenschaftliche Berichte, 2004, SAM-ECOSTAR-D15/FZKA 7064

12 Steinwarz W, Alemberti A, Häfner W, Investigations on the phenomenology of ex-vessel core melt behavior (COMAS). Nuclear Engineering and Design, 2001, 209(1–3): 139–146

13 Tromm W, Foit J J, Magallon D. Dry and wet spreading experiments with prototypic materials at the FARO facility and the theoretical analysis, 2000, Germany. FZKA, 2000, 6475: 178–188

14 Journeau C, Boccaccio E, Brayer C, Ex-vessel corium spreading: results from the VULCANO spreading tests. Nuclear Engineering and Design, 2003, 223(1): 75–102

15 Alsmeyer H, Albrecht G, Fieg G, Controlling and cooling core melts outside the pressure vessel. Nuclear Engineering and Design, 2000, 202(2–3): 269–278

16 Ogura T, Matsumoto T, Miwa S, Experimental study on molten metal spreading and deposition behaviors. Annals of Nuclear Energy, 2018, 118: 353–362

17 Matsumoto T, Sakurada K, Miwa S, Scaling analysis of the spreading and deposition behaviors of molten-core-simulated metals. Annals of Nuclear Energy, 2017, 108: 79–88

18 Ogura T, Matsumoto T, Miwa S, Experimental study on molten metal spreading and deposition behaviors on wet surface. Progress in Nuclear Energy, 2018, 106: 72–78

19 Sahboun N, Miwa S, Sawa K, A molten metal jet impingement on a flat spreading surface. Journal of Nuclear Science and Technology, 2020, 57(9): 1111–1120

20 Yokoyama R, Suzuki S, Okamoto K, Scale effect of amount of molten corium and outlet diameters on corium spreading. Progress in Nuclear Energy, 2020, 130: 103535

21 TEPCO. Locating fuel debris inside the unit 3 reactor using a muon measurement technology at Fukushima Daiichi nuclear power station. 2017–9–28, available at the website of tepco.co.jp

22 Huppert H. The propagation of two-dimensional and axisymmetric viscous gravity currents over a rigid horizontal surface. Journal of Fluid Mechanics, 1982, 121: 43–58

23 Foit J J. Spreading under variable viscosity and time-dependent boundary conditions: estimate of viscosity from spreading experiments. Nuclear Engineering and Design, 2004, 227(2): 239–253

24 Farmer T, Sienicki J, Chu C, The MELTSPREAD-1 computer code for the analysis of transient spreading and cooling of high temperature melts. Report EPRITR-103413,1993

25 Wittmaack R. Coreflow: a code for the numerical simulation of free-surface flow. Nuclear Technology, 1997, 119(2): 158–180

26 Spindler B, Veteau J M. The simulation of melt spreading with THEMA code. Nuclear Engineering and Design, 2006, 236(4): 415–424

27 Allelein H, Breest A, Spengler C. Simulation of core melt spreading with LAVA: theoretical background and status of validation. In: Proceedings of the OECD Workshop on Ex-vessel Debris Coolability, Karlsruhe, Germany, 1999, FZKA6747

28 Koshizuka S, Oka Y. Moving particle semi-implicit method for fragmentation of incompressible fluid. Nuclear Science and Engineering, 1996, 123(3): 421–434

29 Kawahara T, Oka Y. Ex-vessel molten core solidification behavior by moving particle semi-implicit method. Journal of Nuclear Science and Technology, 2012, 49(12): 1156–1164

30 Li G, Oka Y, Furuya M. Experimental and numerical study of stratification and solidification/melting behaviors. Nuclear Engineering and Design, 2014, 272: 109–117

31 Chai P, Erkan N, Kondo M, Experimental research and numerical simulation of moving molten metal pool. Mechanical Engineering Letters, 2015, 1: Paper No.15–00367

32 Yamaji A, Li X. Development of MPS method for analyzing melt spreading behavior and MCCI in severe accident. Journal of Physics: Conference Series, 2016, 739: 012002

33 Chai P, Kondo M, Erkan N, Numerical simulation of MCCI based on MPS method with different types of concrete. Annals of Nuclear Energy, 2017, 103: 227–237

34 Jubaidah D G, Duan G, Yamaji A, Investigation on corium spreading over ceramic and concrete substrates in VULCANO VE-U7 experiment with moving particle semi-implicit method. Annals of Nuclear Energy, 2020, 141: 107266

35 Kondo M. A physically consistent particle method for incompressible fluid flow calculation. Computational Particle Mechanics, 2021, 8(1): 69–86

36 Kondo M, Ueda S, Okamoto K. Melting simulation using a particle method with angular momentum conservation. In: Proceedings of the 2017 25th International Conference on Nuclear Engineering, Shanghai, China, 2017

37 Monaghan J J. Simulating free surface flows with SPH. Journal of Computational Physics, 1994, 110(2): 399–406

38 Spencer D B, Mehrabian R, Flemings M C. Rheological behavior of Sn-15 Pct Pb in the crystallization range. Metallurgical and Materials Transactions B, Process Metallurgy and Materials Processing Science, 1972, 3(7): 1925–1932

39 Joly P A, Mehrabian R. The rheology of a partially solid alloy. Journal of Materials Science, 1976, 11(8): 1393–1418

40 Ramacciotti M, Journeau C, Sudreau F, Viscosity models for corium melts. Nuclear Engineering and Design, 2001, 204(1–3): 377–389

41 Kondo M. Development of surface tension model with many-body potential. In: International Conference on Particle-based Methods – Fundamentals and Applications (Particles 2017), Hannover, Germany, 2017, 461–470