Abstract
The moving particle semi-implicit (MPS) method as a Lagrangian method is attracting increasing attention in severe accident analysis. In this paper, we developed an MPS code for the corium behavior analysis with several additional models added: an improved heat transfer model to improve the calculation between different materials, an enthalpy-based viscosity model to realize a smooth transition of viscosity at the solid-liquid interface, and a surface tension model for better simulation of surface shape. Validation of the developed simulation approach is carried out on a classical water column collapse example. The development of the heat transfer model is validated by the example of a one-dimensional semi-infinite plate. A comprehensive example of the melting of “Wood’s alloy” is carried out to verify the capacity of MPS method in the simulation of melting and expansion. The simulation results are in good agreement with the experimental results, which indicates that MPS method promises well in the field of severe accidents.
1. Introduction
During severe accidents, it is important to accurately predict the total thermal hydraulic response of the power plant. And the evolution of the core corium, including the formation, migration, and coolability, is particularly critical. When the molten corium falls into the lower head of the pressure vessel, it reacts violently with the coolant inside the lower head, which is called the fuel-coolant interaction (FCI). After the FCI, the corium breaks down into many small and irregularly shaped particle fragments, which accumulate in the lower head to form the debris bed [1, 2]. Failure of the pressure vessel if the debris bed is not cooled effectively can result in the discharge of the core melt from the pressure vessel, which will cause thermal erosion of the containment structure and eventually lead to loads of late containment failure. Therefore, the research of debris bed is important in the phenomenology of nuclear severe accidents.
The physical phenomena of the debris bed are very complex, and it is difficult to make analytical or semianalytical calculations in theory. The study on the mechanism should be carried out mainly through experiments and numerical simulations. There are many experiments on debris beds. FCI experiments in FARO experimental facility [3], FCI experiments in KROTOS experimental facility [4], and CCM experimental facility in Argonne National Laboratory [5] all concluded that the median diameter of the debris particles formed in FCI varied in a wide range, and the porosity of the debris bed was not uniformly distributed in space. To prevent the debris bed from remelting in the lower head due to decay heat, water is injected into the dry debris bed. Royal Institute of Technology completed the process of reflooding of debris bed on the POMECO experimental facility [6]. Booker Haven National Laboratory [7] analyzed the impact of initial debris temperature and coolant injection flow rate on the reflooding process through experiments. In all, the research studies on the debris bed mainly focus on the formation and cooling of debris bed. The experiments on the cooling of the debris bed also support the development of the simulation program. In the simulation of debris bed formation, the results obtained by IKEJET/IKEMIX program were in good agreement with FARO experiment [8]. IKE developed a program module WABE [9] that can calculate the cooling process of debris bed in the lower head based on the porous media theory. Some well-known severe accident analysis software (such as RELAP5/SCDAP, MAAP, and MELCOR) can also be used for the calculation of debris bed cooling analysis.
In computational fluid dynamics (CFD), the grid method is relatively mature, and almost all mainstream CFD commercial software is developed based on the finite volume method (FVM). The application of the grid method based on Eulerian coordinates is very convenient, but the grid also forms a natural constraint on solving problems such as moving phase boundary and free surface deformation. For the simulation of melting and deformation phenomena which involves free surface and solid-liquid two-phase, the traditional grid method can mainly employ the following methods: finite element method with moving grid [10], fluid volume method (VOF) [11], level set method [12], and variational domain variational method [13], etc. It is necessary to redefine the computational boundary in each time step and rearrange the grids by means of grid reconstruction, dynamic grid, etc., to achieve effective interface capture and boundary layer calculation, which will bring large calculation costs. If the melting involves a variety of immiscible materials, the identification of interfaces between materials will also become a difficulty in grid simulation. Because of numerical dissipation, the clarity of the phase interface and free surface is difficult to guarantee.
Lagrangian methods show advantages in theory for solving these problems. In the Lagrangian coordinates, there is no convection term in the governing equations and therefore no numerical dissipation problem. Particle methods based on Lagrangian coordinates [14, 15] were proposed in the 1970s. The main difference between different particle methods lies in the selection of different approximate function construction methods and the discrete form of the governing equations. The moving particle semi-implicit method (MPS) in this paper is a particle method that uses a kernel function and point collocation to approximate the physical problem. The MPS method was proposed in 1996 by Koshizuka et al. [16], and then first used in the field of nuclear engineering for the simulation of the crushing of molten metal in the process of steam explosion [17]. Tian et al. [18] tracked the topological shape changes when bubbles burst in the liquid phase. In addition, MPS has also been involved in eutectic interactions related to severe accidents in nuclear power plants and the behavior of molten pool in the lower head in recent years [19, 20]. In the study of accidents outside pressure vessels, MPS method is mainly used to simulate the molten core-concrete interaction (MCCI) and the phenomenon of corium spreading [17, 21, 22]. Based on this algorithm, heat transfer, phase transition, melt viscosity, and thermal radiation have been gradually considered. In a word, the current application of MPS in the field of severe accidents in nuclear power plants focuses on its advantages as a particle method, such as dealing with problems with moving two-phase interfaces and with large liquid deformation.
Interested in the potential capability of the MPS in modeling free surface heat transfer problems, we independently developed an MPS code. It solves the Navier–Stokes equations with additional models for the heat transfer, phase transition, viscosity, rigid body motion, and surface tension. The heat transfer model is further improved and validated. The MPS program is verified by a classical example of water column collapse. Finally, the simulation of Wood’s alloy melting is carried out to verify that the multimodel-coupled MPS method can effectively calculate the melting and spreading behavior of materials.
2. Numerical Methods
2.1. The Basis MPS Method
The governing equation of MPS is composed of continuity equation, Navier–Stokes equation, and energy equation.where the vectors u, f, and F, respectively, are velocity, surface tension, and external force (gravity); t, ρ, p, and υ, respectively, are time, density, pressure, and kinematic viscosity; h, K, T, and Q, respectively, are enthalpy, thermal conductivity, temperature, and internal heat source.
Particle interactions are defined by the kernel functionwhere is the distance between particle and particle .
The gradient and Laplacian terms could be discretized using the following schemes.where and represent scalar, is the constant particle density calculated in the initial state, and is the spatial dimension.
2.2. Development of the Heat Transfer Model
The energy equation is discretized aswhere , is the enthalpy of particle at time step , , respectively.
Equation (6) only applies to the case of constant thermal conductivity k, which means it only applies to homogeneous material. For the heat transfer between different materials, Kawahara et al. [23] put forward the numerical calculation equation of MPS
However, the temperature distribution with equation (7) near the contact surface has a little deviation (seenin Figure 1).
Considering this deviation, we use other way to deal with heat transfer between different materials. We can directly modify the thermal conductivity k in (4). In the perspective of the thermal balance, it can be deduced that when the heat flow is across two different materials, it is appropriate to express the thermal conductivity by the harmonic-mean of the thermal conductivity of the two materials.where , is the thermal conductivity of material 1, material 2, respectively. With the modified heat conductivity, the heat conduction of the heterogeneous material can be calculated.
2.3. Enthalpy-Based Phase Change Model
The phase change model reflects the linear relationship between temperature and enthalpywhere , is, respectively, solidus temperature, liquidus temperature; , is, respectively, the enthalpy at the beginning of the melt, at the end of the melt. , is the heat capacity in the solid phase, in the liquid phase, respectively.
The solid fraction is also calculated based on the enthalpywhere means the particle is completely solid, and means the particle is completely liquid.
2.4. Viscosity Model
The viscous term “” on the right of the momentum equation (2) can be calculated explicitly as follows:
In the process of phase change, the viscosity coefficient υ often changes with temperature dramatically. To describe the process, an empirical formula for calculating viscosity is usedwhere is the kinematic viscosity in liquid phase, is the empirical constant, and is the solid fraction of the melt.
2.5. Surface Tension Model
The potential function [24] is as follows:where , , and , respectively, are the distance between particles, the distance between initial particles, and the effective radius, which is generally set as .
is a coefficient related to surface tension, which can be obtained by the surface energy between fluids.where is the coefficient of surface tension.
By the derivative of the potential function, we can get the force between the particles related to the surface tension.
2.6. Calculation Process
The calculation process calculation is shown in Figure 2. The pressure term is implicitly calculated with the pressure Poisson equation:
According to the pressure distribution, the velocity and coordinate of particle in this time step are obtained:
In this paper, conjugate gradient algorithm is adopted as the solver of pressure Poisson equation, so that OpenMP can be used in parallel to improve the computational efficiency.
3. Results and Discussion
In this section, a classic example of water column collapse is calculated with large deformation and low viscosity. Then, to validate the new heat transfer model, the calculation of a one-dimensional semi-infinite plate is carried out. At last, Wood’s alloy melting experiment is simulated to verify the MPS method coupled with heat transfer, phase change, viscosity, and surface tension calculation. The effective simulation of the melting experiment will lay the groundwork for more complex corium calculations with the MPS method.
3.1. Water Column Collapse
The case of water column collapse has been widely used in CFD validation. The height of the container is 29.2 cm, and the length of the bottom is 58.4 cm. The water column, which stands still on the left side of the container, is equal to the height of the container, and the width is 1/4 of the bottom of the container. Water begins to collapse under gravity from . (seen in Figure3)
With the initial particle distance of and time step of , the simulation has a total time of . In Martin’s experiment [25], when and the size of water column (equivalent to parameter L in Figure 2) was 1.125 in and 2.25 in, respectively, were taken as the experimental reference value. The dimensionless number of fluid front positions was compared with the experiment as Figure 4, and the results of volume of fluid (VOF) [11] are also presented as a reference of simulation.
As shown in Figure 4, the trend of fluid front changing with time in MPS calculation is consistent with the experimental results, but the front motion changes slightly faster in simulation. Such a difference is reasonable considering that it takes time for the water column barrier to be extracted in the experiment, and there is friction on the surface of the solid wall. All of these will hinder the process of water column collapse. The simulation results of water column collapse by VOF also reflect this difference. The results of MPS are consistent with those of VOF, and the differences between MPS and experiments are also consistent with the actual situation.
3.2. One-Dimensional Plate
The heat transfer calculation is validated by the case of semi-infinite one-dimensional plate. The geometry of the plate is shown in Figure 5. The temperature at the left boundary is maintained at , and the initial temperature of other position is .
The length of the plate is long enough to be considered semi-infinite in the computational region. The one-dimensional temperature distribution with time can be obtained from an analytical solution. Given the thermal properties, the temperature distribution at the middle position at is shown in Figure 6. It can be seen from the figure that the temperature calculation results of MPS are in good agreement with the analytical solution.
For heat transfer of heterogeneous material, the midst of X direction of the semi-infinite one-dimensional plate is the boundary, and the left and right sides are set with different materials. The initial conditions () of this plate are shown in Figure 7, taking and .
Given the thermal properties, the temperature distribution at the middle position at 10 s is shown as in Figure 1. The temperature distribution calculated by MPS with equation (7) near the contact surface is quite different from the reference value, as shown in the black dotted line. Using equation (6) with harmonic mean of thermal conductivity equation (9), the unsteady thermal conductivity between different materials is calculated again and the results are in the red dotted line. It shows that the temperature distribution near the contact surface of different materials is clearly improved, and the numerical solution agrees well with the analytical solution. Therefore, the improvement of thermal conductivity between different materials is proved to be effective.
Heat transfer coupled with phase transition in MPS is validated by Stefan’s problem [26]. At the initial time, the temperature of the liquid is and the surface of the liquid at suddenly drops to ( is below the melting point) and stays at . Therefore, the liquid-solid interface moves in the positive X direction from . The temperature distribution at the middle position of the one-dimensional semi-infinite plate at and is shown in Figure 8.
From Figure 8, the temperature distribution obtained by MPS is in good agreement with the analytical solution. The positions of the solid-liquid interface can be calculated as and from the analytical solution. It shows that the positions of the solid-liquid interface calculated by MPS perfectly match the calculation results of the analytical solution.
3.3. Three-Dimensional Alloy Melting
The process of alloy melting heated by the bottom plate is simulated as an integrated case. The initial situation of the alloy is shown in Figure 9, and the solid alloy is placed on the heating plate at the initial moment. With a fixed heating temperature, this melting process lasts about 2 s and stops when the melt interface is about to reach the solid center.
A three-dimensional model is used and the convection heat transfer of air is ignored. The solid part of the alloy is modified by the rigid body model. The simulation results are compared with the experiment results conducted by Guo et al. [27]. The front view and top view of the simulation results compared with the experimental images are shown in Figures 10 and 11, respectively. The simulation images are in good agreement with the experimental images, and the simulated particles clearly show the free interface shape of the molten alloy.
(a)
(b)
Figures 12–14 compare the simulation results of MPS method of solid phase alloy height, solid phase alloy central temperature, and liquid phase alloy width with the experimental results, and the simulation results obtained by Guo et al. [27] using the finite volume particle method (FVP) are shown in Figures.
In Figure 12, the simulation results of MPS method for solid phase alloy height are generally close to the experimental results, and the relative error is within an acceptable range. The result of MPS method is close to that of FVP, and the curve of MPS method is less volatile in the first half period. The height of solid alloy is greatly influenced by the rate of phase transition and the thickness of the molten alloy under the solid alloy. However, the experimental solid alloy height did not change much in the initial 0.3 s, while the simulated alloy height decreased from the beginning of the calculation. This is mainly because the experiment timing begins with the alloy being put on the heating plate, but at the initial moment, due to the roughness of the heating plate surface and metal surface, the contact did not achieve full joint; thus the initial thermal conductivity between alloy and the heating plate is actually lower than that of the numerical calculation. This makes the height decline of the alloy slow at the initial stage, and the numerical simulation cannot reflect this detailed process.
Figure 13 shows the change of the central temperature of the solid phase alloy over time. The heat transfer in the solid is relatively slow. The temperature at the center of the solid alloy does not increase significantly until the melt interface is close enough to the center of the solid alloy. In the temperature curve, the first half of the curve is almost parallel to the time axis, while the slope of the second half of the curve increases gradually. The central temperature calculated by MPS method is close to the experimental results, which is clearly better than that calculated by the FVP method.
Figure 14 shows the variation of liquid phase alloy width over time. On the whole, the variation of extension width calculated by the MPS method is more similar to the experimental results than the calculated results of FVP. Especially in the second half of the time, the calculation results of the FVP method show a trend that the slope is close to 0, which is not quite consistent with the experiment. And the calculation results of MPS method show that the elongation width is still increasing. For the difference between the results of MPS method and the experimental results, it is speculated that the exponential increase viscosity model does not fully comply with the change of alloy viscosity in the experiment, which is the main reason. To obtain a more consistent viscosity model, more experiments may be required to fit the viscosity change function.
In the calculation of alloy melting, the coupling of different physical models shows good simulation results. The results of MPS simulation show good similarity with the experimental photos, and the characteristic parameters of solid alloy and molten alloy show good consistency between the simulated and experimental results. This indicates that the MPS method can well reproduce the flow behavior of Wood’s alloy with the phase transition after heating.
4. Conclusions
To simulate the problems related to corium behavior, MPS numerical algorithm is established in this paper, and some models are coupled into the MPS. The heat transfer model of homogeneous material is introduced and with some improvements; the heat transfer model can also be used in heterogeneous materials. To simulate the solid-liquid coupling problem effectively, a phase transition model and viscosity model based on enthalpy are added. To simulate the shape of the free surface effectively, the surface tension model based on potential function is added.
By simulating the water column collapse with low viscosity and large deformation, the effectiveness of MPS method in simulating the free surface problem is validated. The heat transfer model and phase transition model are validated by the simulation of one-dimensional semi-infinite plate. The temperature distribution near the contact surface can be calculated accurately. Finally, the thermal melting process of Wood’s alloy is simulated, and the simulation results show that MPS method can reproduce the process of phase transition and flow behavior of Wood’s alloy that is heated on board.
Overall, the developed MPS code with the added features enables its capability in modeling the phenomena of melt behavior. Further studies on specific analysis of melt behavior in nuclear engineering applications would be considered. Although the development of MPS method is still very immature, it has the advantage of overcoming the grid limitation of traditional Euler methods in the simulation of complex scenarios related to severe accidents. This new algorithm may bring a new perspective to the problems like the formation and remelting of debris beds, the formation and rupture of crust in corium spreading, and the structure and characteristics of porous media in the process of corium cooling.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors are grateful for the financial support of the Fundamental Research Funds for the National Key R&D Program of China (Project no. 2018YFB1900100).