- Email: [email protected]
Numerical analysis of the melt behavior in a fuel support piece of the BWR by MPS
Annals of Nuclear Energy 102 (2017) 422–439
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Numerical analysis of the melt behavior in a fuel support piece of the BWR by MPS Ronghua Chen a, Lie Chen a, Kailun Guo a, Akifumi Yamaji b, Masahiro Furuya c, Wenxi Tian a, G.H. Su a,⇑, Suizheng Qiu a a b c
School of Nuclear Science and Technology, Xi’an Jiaotong University, Xi’an 710049, China Cooperative Major in Nuclear Energy, Graduate School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan Central Research Institute of Electric Power Industry (CRIEPI), 2-11-1 Iwado-kita, Komae, Tokyo 201-8511, Japan
a r t i c l e
i n f o
Article history: Received 23 June 2016 Received in revised form 5 January 2017 Accepted 8 January 2017 Available online 17 January 2017 Keywords: MPS method Fuel support piece Severe accident Freezing behavior
a b s t r a c t The fuel support piece in a boiling water reactor (BWR) is used to brace fuel assemblies. The channel within the fuel support piece is determined to be a potential corium relocation path from the core region to the lower head during the severe accident of BWR. In the present study, the improved ⁄⁄⁄Moving Particle Semi-implicit (MPS) method was adopted to simulate the flow and solidification behavior of the melt in a fuel support piece. The MPS method was first validated against the Pb-Bi plate ablation test that was performed by CRIEPI. The predicted ablation mass of the plate agreed well with the experimental results. Then the flowing and freezing behaviors of molten stainless steel (SS) and zircaloy in the fuel support piece were simulated by MPS method with a three dimensional particle configuration, respectively. In this study, the flow and solidification behavior of SS was simulated first. After all the SS passed through the channel, the flowing behavior of Zr in the fuel support piece was simulated. The simulation results indicated that the crust layer formed on the inner surface of the fuel support piece during the melt discharging process. The fuel support piece was plugged by the solidified zircaloy particles in the lower initial temperature case. The fuel support piece kept intact in all the calculation that were performed under the assumed order of melt injection. The present results could help to reveal the progression of a BWR severe accident. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Due to the battery power depletion and station blackout, few plant data such as pressure, temperature, and water level changes can be obtained during the accident at the Fukushima Dai-ichi Nuclear Power Station (Sakai et al., 2014). After the accident, local measurement is impossible due to high level radiation in the reactor building. Because of the above reasons, it is still unclear how the corium flows into the lower plenum and finally causes the vessel failure. During the late phase of sever accident in a boiling water reactor, a melt pool may form in the reactor core. There are five potential corium relocation paths from the core region to the lower head, as shown in Fig. 1, including the fuel support piece, control rod guide tube, instrumentation tube, the shroud wall and core support plate failure (Sakai et al., 2014). Because of the impossibility of local measurement, a numerical calculation and analysis is needed ⇑ Corresponding author. E-mail address: [email protected] (G.H. Su). http://dx.doi.org/10.1016/j.anucene.2017.01.007 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
and the result will help to find out the progress of the severe accident occurred in the Fukushima Dai-ichi nuclear power station. During the flowing of high temperature corium in fuel support pieces, great thermal attack will be brought into these channels. The corium may ablate the wall of channel which leads the fuel support piece to be damaged, or it may freeze and plug the channel if the heat dissipation condition of corium is quite well. The calculation of the flowing and heat transfer of corium can assess the integrity of fuel support piece under the severe accident. On the other hand, the evaluation of the flowing and solidification behavior of corium in fuel support piece will help to reveal the relocation path of the corium in the core, which could add knowledge on the research of severe accident in BWR. Whether the molten core debris will lead to the fuel support pieces failure or freeze and plug these channels is a significant issue that need to be investigated. In respect of this issue, several relevant experiments have been performed. Berthoud and Duret (1989) designed GEYSER which is an upward melt injection experiment to research the flowing and solidification behavior of molten UO2 in a small tube. Sandia National Laboratories performed the XR2-1 test to study the
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
423
Nomenclature A C T W a d f ~ f ! g h k l0 m n p ! r r t ! u
surface area, m2 specific heat, J/(kgK) temperature, K weight function thermal diffusivity, m2/s dimension number solid fraction of particle surface tension per unit mass gravity acceleration, m/s2 enthalpy, J/kg curvature, m1/coefficient of thermal conductivity, W/ (mK) particle size, mm mass, kg particle number density pressure, Pa particle location vector particle distance, m time, s particle velocity vector
corium relocation pathways and processes in BWR core under the condition of severe accident. XR2-1 test results indicated that the molten corium pool had not formed in the core under the ‘‘dry core” condition (Gauntt and Humphries, 1997). Rahman et al. (2007) performed a series experiments to figure out the basic physical characteristic in freezing mode of molten metal and characteristic parameters of the freezing during penetration into the structure. Experimental results such as penetration length and width of frozen metal on stainless steel and brass structures were compared under the air and water coolant conditions. Apart from experiments of the melt flowing and freezing behavior, an inspection in the TMI-2 nozzle and instrument tube specimens have been performed (Wolf et al., 1994) and the results indicated that the molten core debris had traveled a long distance outside the reactor pressure vessel through the instrument tube, and the coolant in the annulus prevented the core melt from melting through the in-core housing tube and leaking into the containment (Wolf et al., 1994). For numerical simulation of the corium flowing and solidification behavior, several models and methods have been put forward.
Fig. 1. Corium flow paths from the core in BWR.
Greek letters heat source, W/kg e emissivity / radiation heat transfer / scalar variable k diffusion coefficient l viscosity, Pas q density, kg/m3 r surface tension/stefan-boltzmann constant m kinematic viscosity, mm2/s
U
Superscript/subscript a ambient e effective i; j particle No. l liquid phase s solid phase min minimum w wall 0 initial condition ⁄ temporal
Cheung et al. (1992) developed an effective diffusivity approach to simulate the heat transfer in a horizontal heat-generating layer. Bui and Dinh (1996) put forward the Effective Convectivity Conductivity Model (ECCM) to describe the heat transfer. Ostensen (1973) developed the bulk freezing model and the conduction layer model. These two models were used to predict the penetration length of the molten core debris. In addition to these mechanism models, some analyses based on severe accident integral code (e.g. ASTEC Groudev et al., 2013, MAAP (Users Manual, 1999) and MELCOR (Gauntt et al., 2005) have been done as well as some CFD analyses. Morita et al. (2011) used COMPASS code which is a kind of particle methods to simulate the GEYSER and THEFIS melt penetration experiments. Li et al. (2016) analyzed the melt behavior in a BWR lower head during LOCA and SBO accident with MELCOR 2.1. By introducing a potential lower head failure mechanism, the melt behavior in BWR lower head and leakage paths under postulated LOCA with complete loss of ECCS and SBO accident were sensitively analyzed. In the 1990s, Koshizuka and Oka proposed the Moving Particle Semi-implicit (MPS) method Koshizuka and Oka, 1996. In this kind of particle method, the kernel function and the Lagrange approach are used to trace the particle movements in the dispersed phase. It is suitable for the calculation of incompressible fluid. With the continuous improvement of this method, it has been successfully applied to various fields to solve the problems with lager deformation (Koshizuka et al., 1998). In recent years, the MPS method has a lot of development. The weak compressibility is adopted by Ikeda et al. (2001) to enhance the stability, and a higher order source term for the Pressure Poisson Equation (PPE) is deduced by Khayyer and Gotoh (2009). A higher order Laplacian model for discretizing the Poisson Pressure Equation and a corrective matrix for pressure gradient model were proposed to improve the stability and accuracy of the moving particle semi-implicit method by Khayyer and Gotoh (2011). Jeong et al. (2013) applied an appropriate gradient model and 4th-order Runge–Kutta time integration scheme to carry out numerical simulations for the Rayleigh–Taylor instability and found that the simulation results became closer to the experimental ones. The improved MPS method is also utilized in multiphase-flow, an Euler Lagrange method MPS-MAFL is introduced by Yoon et al., 2001, and Tian et al. (2009), Tian et al. (2010), Chen et al. (2011), Li et al. (2013) used the MPS-MAFL to do a lot of
424
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
research on bubble dynamics. An enhanced stabilized MPS method was developed by Khayyer and Gotoh (2013) for simulation of multiphase flows characterized by high density ratios. Shirakawa et al. (2002) analyzed flows around a BWR spacer with a twofluid particle interaction method (TF-MPS) and the result shows an agreement with the experiment which suggests a certain applicability of this method. Chen and Oka (2014) analyzed the melt behavior in a cold tube with MPS method in which the change of melt viscosity with temperature was taken into account. Beside these, the turbulence model is also adopted in MPS method, many researchers simulated the complex flow by MPS-LES (Gotoh et al., 2001; Song and Hitoshi, 2005; Shao et al., 2006; Duan and Chen, 2015). In this study, the improved MPS method was used to analyze the corium flowing and solidification behavior in a fuel support piece during a sever accident of BWR. At first, the models used in the improved MPS were validated against the analytical solutions and then the thin plate ablation test performed by CRIEPI was simulated and the feasibility of analyzing corium flowing and solidification behavior for this method was validated by comparing the calculation results with experimental results. Then the fuel support piece was simulated by a three-dimensional particle configuration and the behavior of corium in this structure was analyzed under different initial conditions. In this calculation, the SS was assumed to be melted first, and then the Zr was melted and flowed into the fuel support piece when all the SS passed through the channel. 2. The MPS method The basic theory of Moving Particle Semi-implicit (MPS) method is that the macroscopical fluid is represented by a series of discrete particles and the differential operator in governing equation is expressed by forms of particle interaction through corresponding computational model. 2.1. Governing equations and computational models As is often the case in problems of incompressible viscous fluid, the continuity, Navier–Stokes and energy equations are adopted to describe the flow characteristic or the flow field. The continuity equation is presented as follow,
Dq ¼0 Dt
ð1Þ
In MPS method, the Lagrange reference coordinates is used, therefore there is no convection term in momentum equation which is advantageous to the numerical calculation.
q
D~ l l þ q~f þ q~ g ¼ rP þ lr2 ~ Dt
ð2Þ
In the above momentum equation, from left to right, the terms are transient term, pressure gradient term, viscous term, surface tension term, and gravity term. In the calculation, the surface and the gravity terms were solved explicitly. In order to guarantee good convergence and accuracy, the pressure and the viscous terms were solved implicitly.
weight function (Kondo and Koshizuka, 2011) is expressed as follow,
(
WðrÞ ¼
ni ¼
X ! ! Wðj rj r i jÞ
ð5Þ
j–i !
!
r i and r j are the position vectors of particle i and j, respectively. For the incompressible fluid, the density is a constant. In the MPS method, fluid density is characterized by the particle number density, therefore the particle number density for incompressible fluid which is represented by n0 is constant as well.In 1996, Koshizuka and Oka developed the gradient model and Laplacian model (Koshizuka and Oka, 1996) for discretizing the differential operator in the conservation equations.
hruii ¼
! ! d X uj ui ! ! ðr j r i ÞWðrj r i Þ n0 j–i j r!j r!i j2
hr2 uii ¼
! ! 2d X ðu ui ÞWðr j r i Þ kn0 j–i j
ð6Þ
ð7Þ
To ensure the astringency of this method for high viscosity calculation, an implicit scheme was introduced to solve the viscous term. By using the Laplacian model, the viscous term was implicitly discretized as follow, ! u i
!
¼ u ni þ tDt
! ! 2d X ! ! ðu u i ÞWðj r nj r ni jÞ kn0 j–i j
ð8Þ
!
The temporal particle velocity, u ; can be obtained by solving the Eq. (8). On the other hand, to avoid the pressure oscillation during the calculation, the solving of Poisson equation of pressure in MPS method is improved (Kondo and Koshizuka, 2011). The source terms of Poisson equation is changed as,
1
q0
r2 P ¼
1 b n 2nn þ nn1 b c n nn c n n0 þ þ 2 2 2 0 0 n n n0 Dt Dt Dt ð9Þ
Where b and c are used to adjust the compensating parts and they should obey the relationship 0 6 c 6 b 6 1. A more proper source term of Poisson Pressure Equation was proposed by Khayyer and Gotoh (2011) with dynamic coefficients as functions of instantaneous flow field. The differential operator r2 P was discretized by using the Laplacian model mentioned above,
r2 P ¼
! ! 2d X ðP j Pi ÞWðr j r i Þ kn0 j–i
ð10Þ
and the final discretized Poisson equation is as follows,
X ! ! ðP i Pj ÞWðr j ri Þ
ð3Þ
Terms in the energy equation are transient term of enthalpy, conduction term and internal heat source term. In the theory of MPS method, the scope of action for one particle is limited by an effective radius re : Considering the influence of the distance from the central particle, a weight function is introduced into the calculation of this method. In this study, the adopted
ð4Þ
re P r
0;
A parameter named particle number density which is proportional to the fluid density is defined as,
j–i
@h ¼ kr2 T þ U @t
ð1 r=r e Þ2 ; 0 6 r < r e
¼
kq0 n0 1 b n 2nn þ nn1 b c n nn c n n0 þ þ 2 2 2 0 0 0 2d n n n Dt Dt Dt ð11Þ
The surface tension model in the program was developed by Kondo et al. (2007) which was based on the inter-particle potential force, !
f ¼ Kðr r min Þðr re Þ=m
ð12Þ
425
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
where K; r, r min are fitting coefficient, particle distance, and initial minimum particle distance, respectively. The potential force coefficients between two fluid particles and between a fluid particle and a solid particle are given as follows,
X
Kf ¼
1 3
2rr 2min
r 32 r min þ 12 re ðr re Þ2
ð13Þ
i2A;j2B;jr ij j
K fs ¼
1 K f ð1 þ cos hÞ 2
ð14Þ
Where h is the contact angle. In Kondo’s model, the surface tension coefficient is same as the surface energy. The surface energy is the potential energy which is stored on a unit surface. In other words, the required work on creating a unit surface is equivalent to the surface energy. As shown in Fig. 2, A and B are the two areas that are used to create two new fluid surfaces for the calculation of surface tension in Kondo’s model as presented in Kondo et al. (2007). In order to simulate the melting and solidification process, the heat transfer and phase change model was incorporated into the program. Based on the energy equation, the enthalpy of particles is calculated with the equations given as follows, kþ1
hi
k
¼ hi þ
2d X kij ðT kj T ki Þwðj~ rj ~ r i jÞ Dt þ U Dt n0 k j–i
ð15Þ Fig. 3. Calculation algorithm of MPS method.
where, the kij is the average coefficient of thermal conductivity,
kij ¼
2ki kj ki þ kj
ð16Þ
The temperature of every particle will be calculated based on its enthalpy.
8 T þ hhs0 > > < s qC ps T ¼ Ts > > : T þ hhs1 s
qC pl
ðh < hs0 Þ ðhs0 6 h 6 hs1 Þ
ð17Þ
ðhs1 < hÞ
T s is the melting point and hs0 , hs1 are the enthalpies when the particle begins and stops melting. The phase change model was introduced as follow, 0
f ¼
8 > <1
hs1 h > hs1 hs0
:
0
ðhs1 < hÞ
2.2. Calculation process The flow chart of MPS method is shown in Fig. 3. In one time step, the energy equation is solved firstly to obtain the temperature of particles. Based on the particle temperature, the viscosity of the particle is calculated. Then the viscous term of the Navier– Stokes equation is discretized as Eq. (8) and solved to obtain the !
temporal velocity u . Then gravity term and surface tension term !
are explicitly calculated to update theu . With the temporal veloc!
ðh < hs0 Þ ðhs0 6 h 6 hs1 Þ
0
f is the solid fraction of particle and different types of particle were introduced to characterize different phase states.
ð18Þ
ityu , the temporal particle position is calculated by the following equation. !
!
!
r ¼ r n þ Dtu
ð19Þ
!
After the r is obtained, the temporal particle number density n is solved by Eq. (5). Then, the particle pressure is updated by implicitly solving the Poisson equation of pressure with conjunction gradient method. And lastly, the particle velocity and position are updated using the new pressure gradient as follows, ! nþ1
u
! nþ1
r
!
¼ u !
¼ r
Dt
q0
rP
Dt 2
q0
rP
ð20Þ
ð21Þ
3. Validation of the MPS method
Fig. 2. The areas of A and B.
In this study, some changes had been made to the MPS method, including implicit viscous term calculation model and the stability improvement model. The improved MPS method was employed to simulate the classical dam break problem to validate its correctness. Fig. 4 shows the initial particle configuration of the dam break calculation. The particle total number of this simulation is
426
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 4. Particle configuration for dam break problem.
2164, and particle diameter is 0.01 m, the wall condition is set as no-slip boundary condition and a free surface boundary was applied at the surface of the fluid. When the number density of a particle was less than bn0 ; the particle was considered as a free surface particle. Here b was set to be 0.97 and n0 was initial particle number density for inner particles. For the free-surface particles, the Dirichlet boundary condition was used. As the incompressible fluid is not sensitive to absolute pressure, the pressure of the free surface particle was set to be zero. A more accurate method for detecting free-surface particles was proposed by Gotoh and Khayyer (2011) which includes two steps. The first step is to calculate the number density of every fluid particle and find the particle whose density is smaller than the predetermined threshold. The second step is to check symmetry of the particle which is chosen from the first step and the asymmetrical particle is set to be free surface particle. The details for this method is presented in the Gotoh’s reference. This calculation were performed on the computer with Intel Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory. For the original MPS, the initial time step is 104 s, and the time step was adjusted by the CFL condition during the calculation. For the 5.0 s simulation time, the total calculation time is 2 min and 42.19 s. The improved MPS method almost does not influence the computational efficiency. Fig. 5
shows the time histories of the pressure at the test particle, which was marked by circle in Fig. 4. They are predicted by the original MPS method and the improved MPS method with different values of b and c: As shown in Fig. 5, the pressure oscillation can be suppressed by the improved MPS method compared with the original MPS method. The differences of the pressure histories among the three parameter combinations are slight. Thus, the values of b and c were set as 0.5 and 0.4 respectively, in this study. To access the capability of the improved MPS to analyze the heat transfer and the phase change, the one-dimensional solidification problem was simulated by the improved MPS method. Fig. 6 shows the particle configuration of the one-dimensional solidification problem. The wall particle is on the left side with a constant temperature 173 K, the initial temperature of water is 283 K. The heat transfer boundary condition of the water is thermal insulating. The particle diameter shown in Fig. 6 is 0.002 m, and the sensitivity analysis of the particle size is also performed with different particle diameters including 0.001 m, 0.002 m and 0.0025 m. In this simulation, the pressure calculation is not involved. This simulation is performed on the computer with Intel Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory. When l0 equals 0.001 m, the particle total number is 20,000, initial time step is 103 s, the total calculation time is 29 min and 1.003 s for 4.0 s simulation time. When l0 equals 0.002 m, the particle total number is 5000, initial time step is 103 s, the total calculation time is 5 min and 4.244 s for 4.0 s simulation time. When l0 equals 0.0025 m, the particle total number is 3200, initial time step is 103 s, the total calculation time is 3 min and 31.227 s for 4.0 s simulation time. Due to release heat to the cold wall, the water gradually changes its phase from liquid to solid. The analytical solutions of the position of interface, temperature distribution of solid phase, and liquid phase are as follows, Position of interface
pffiffiffiffiffiffiffi nðtÞ ¼ 2R a1 t
ð22Þ
Temperature of solid phase:
T 1 ðt; xÞ ¼ T w þ ðT m T w Þ
erf ð pxffiffiffiffiffiÞ 2
a1 t
erf ðRÞ
Temperature of liquid phase:
pxffiffiffiffiffi 2 a2 t qffiffiffiffi T 2 ðt; xÞ ¼ T 1 þ ðT m T 1 Þ 1 erf R aa12
ð23Þ
1 erf
ð24Þ
Original MPS γ =0.4 β=0.5 γ =0.3 β=0.5 γ =0.4 β=0.6
Pressure/Pa
16000
12000
8000
4000
0
0
1
2
3
4
5
Time/s Fig. 5. Time history of the pressure at test particle.
Fig. 6. The particle configuration of the one-dimensional solidification problem.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439 300
280 273
Freezing point
Temperature/K
260
240
220
l0=0.001m
200
l0=0.0025m
l0=0.002m Analytical
180 0.00
0.05
0.10
0.15
0.20
x/m Fig.7. Results of solidification validation.
where, a1 , a2 are the thermal diffusivity of the wall and water. The result calculated with the program and the analytical solution at 4 s (y = 0.05 m) is presented in the Fig. 7.
427
As shown in the Fig. 7, the temperature distribution predicted by MPS agrees well with the analytical solution. With time going by, water with high temperature froze due to the cold wall at the left side. The solidification process of water and the temperature field of the case during the calculation is shown in Fig. 8. In Fig. 9, the temperature distribution at x = 0.1 m shows that the temperature near the free-surface is not consistent, this is because the number density of particle is less than n0 , thus the Laplacian model should be modified to compensate the PND (particle number density) near free-surface. It should be noted that the Laplacian model introduced by Eq. (7) cannot obtain accurate and consistent results for incomplete compact supports such as free surface. In 2013, Souto-Iglesias et al. improved the Laplacian model based on the analysis of the consistency of the MPS approximation to the Laplacian operator (Souto-Iglesias et al., 2013) and in 2016, Gotoh and Khayyer proposed a more consistent Laplacian model when they calculated the surface tension by considering BI (boundary integrals) (Gotoh and Khayyer, 2016). In this study, we carried out several calculations to evaluate the influence of this kind of Laplacian model. As shown in Fig. 9, we can conclude from Fig. 9 that the error on free-surface is decreasing when the l0 (particle size) is smaller. When l0 = 0.002 m, the relative error is 0.037%. For this reason, in the engineering simulation, the error on freesurface is ignored in this paper.
Fig. 8. Results calculated by MPS (l0 = 0.002 m), (a) solidification process, (b) temperature field.
428
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
l0=0.001m
274.4
l0=0.002m
274.2
l0=0.0025m Analytical
Temperature(K)
274.0 273.8
Table 1 Properties of materials in the thin plate ablation test. Properties 3
Density [kg/m ] Specific heat [J/kg K] Kinematic viscosity [m2/s] Heat conductivity [W/m K] Initial temperature [K]
273.6
Pb-Bi (44.5%, 55.5%)
Wood
Air
10050.0 147.0 1.29E07 13.8 300.0 (plate) Experimental value (inlet)
300.0 1300.0 – 0.069 300.0
1.205 1005 1.511E5 0.0257 300.0
273.4 273.2 273.0 272.8 0.00
0.02
0.04
0.06
0.08
0.10
y(m) Fig.9. Temperature distribution at x = 0.1 m.
4. Validation against the Pb-Bi plate ablation test 4.1. Introduction of the Pb-Bi plate ablation test Central Research Institute of Electric Power Industry (CRIEPI) in Japan had performed a series of plate ablation tests to add the knowledge on the RPV wall ablation mechanism, as shown in Fig. 10. In the plate ablation test, the Pb-Bi was employed as the substitute of the corium and was heated to 350 °C by the electric heaters. At the beginning of the test, the molten Pb-Bi was poured into a wooden vessel. Subsequently, the molten PbBi was discharged from the wooden vessel through the pipe that was connected to the bottom of the vessel. When the molten PbBi left the pipe, it dropt on the Pb-Bi plate. The thermocouples, located in both the wooden vessel and the pipe outlet, were used to measure the temperature of the molten Pb-Bi. Meanwhile, the ablation histories of the Pb-Bi plate were recorded by the video cameras from the side and the top during the experiment processes. The detailed parameters of material properties in this test is shown in Table 1.
Fig. 11. Computational domain of the Pb-Bi plate ablation test.
4.2. Numerical simulation by MPS method 4.2.1. Computational domain and boundary condition As shown in Fig. 11, a three-dimensional particle configuration was constructed for analyzing the plate ablation test. The inlet boundary particles were placed in the same position as the pipe outlet to inject melt particles. The angle between the melt injection direction and the plate is 5 degree. The angle between the plate and the horizontal plane is 15 degree. The plate is 5 mm in thickness, 17 mm in width and 37 mm in length. In this study, the particle diameter is 0.5 mm, and totally 52,626 particles are employed. At the end of the plate, the edge is chamfered which is the same as the experiment.
Fig. 10. Schematic of the plate ablation test facility.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
As to boundary conditions, nonslip conditions were applied at the surface of the Pb-Bi plate, and a free surface boundary was applied at the surface of the fluid. For the surface particles, nature convection heat transfer between particles and air and radiation heat transfer between the surfaces particle and the ambient were taken into consideration. The nature convection heat transfer (Yang and Zhang, 1994) was calculated as,
( Nu ¼
0:59ðPr GrÞ1=4 ;
for laminar flow
0:11ðPr GrÞ1=3 ;
for tubulent flow
ð25Þ
The radiation heat transfer between the surface and the ambient was calculated as,
/ ¼ eArðT 4w T 4a Þ
ð26Þ
where, e is the emissivity and r is the stefan-boltzmann constant. The inflow temperature was set according to the fitting curve based on the experimental data which was shown below,
8 0:0s < t 6 0:35s > < 403; TðKÞ ¼ 18:519t3 123:94t2 þ 293:54t þ 46:008; 0:35s < t 6 3:0s > : 0:0598t3 1:4404t 2 þ 10:578t þ 279; 3:0s < t 6 10:0s ð27Þ The inflow mass rate was approximately set to be 15.22 g/s.
429
4.2.2. Comparison between the numerical calculation and experimental result Fig. 12 shows the comparison of the ablation histories predicted by MPS and the experimental results. In the test, the plate side wall was melted through by the melt at about 4.0 s. The similar phenomenon of side wall melt through was successfully reproduced by MPS as show in Fig. 12. Because of the overshadowing of boundary particles, the melt-though of the sides of the plate couldn’t be presented in the figure. In the simulation of MPS, the edge of the plate was melted though at 4.0 s, and the melt inventory within the plate began to discharge. At 4.5 s, almost all the inventory melt flowed out. On the other hand, the time for the edge being melt through and all the inventory melt flowing out was little later in the experiment and there are two reasons that cause the difference between the simulation and the experiment. First, in the experiment, the transient mass flow rate at the inlet is hard to be accurately measured, thus it is hard to give the exactly same inlet condition in the simulation; Second, a film forms on the surface of the molten Pb-Bi when it is freezing, this film will influence the shape of the molten material. But in the simulation, limited by the computing capabilities, the influence of the film is not considered unless the size of the particle is small enough. Generally, the predicted plate ablation histories agreed with the experimental observation. Fig. 13 shows the temperature and pressure distribution during the simulation period.
Fig. 12. The plate ablation histories predicted by MPS and measured in experiment.
430
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 13. Temperature and pressure distribution of the ablation experiment.
0.14
Inflow mass MPS Outflow mass MPS Outflow mass Exp.
0.12
Outflow mass(kg)
0.10 0.08 0.06 0.04 0.02 0.00
0
1
2
3
4
5
6
7
Time(s) Fig. 14. The mass profiles predicted by MPS and measured in the test.
Fig. 14 shows the inflow mass profile and the outflow mass profile predicted by MPS, and the experimental measured result was also shown in Fig. 14 for comparison. Within three seconds, small difference was observed from the three profiles. After four seconds, the predicted out flow mass began to increase quickly, and finally the outflow mass is consistent with the experimental value. The comparison results indicated that the present MPS method is able to simulate the thin plate ablation test reasonably. Next, the MPS method was employed to analyze the flow and solidification behavior of the melt in the fuel support piece, which has the similar mechanism as the plate ablation test.
Fig. 15. Schematic of fuel support piece.
5. Analysis of the corium behavior in the fuel support piece Each fuel support piece carries four fuel assemblies as shown in Fig 15. Accordingly, at the bottom of the fuel support piece there are four orifices which allow the coolant to flow into the fuel support piece and then pass through it and entry into the fuel assemblies. In the BWR core melting down accident, passing though the fuel support piece is a possible way for the corium to
431
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 16. Computation domain of the fuel support piece.
Table 2 The main parameters of the size for the fuel support piece. Parameter
Value [cm]
Total height of computational domain Outer diameter of upper barrel Height of transition region Outer diameter of lower barrel Height of inclined section Outlet diameter of bottom outlet
9.2 6.0 0.8 4.6 5.4 1.8
relocate from the core region to the lower head. In this study, the MPS method was used to analyze the melt behavior in the fuel support piece.
and a few other parts are made by stainless steel. And the claddings of fuel are made up of zircaloy (Zr). During the severe accident, the control blade would be inserted into the fuel assemblies and the sheath which is made up of SS is possible to reach the core. As the melting point of SS is lower than that of the zircaloy, it may be melted first and flows into the fuel support pieces during severe accident of BWR and with the temperature of the core increase, other materials will be melted later. In the simulation, a melt pool was assumed to be formed within the fuel assembly. The melt was injected into the fuel support piece from the top driven by the pressure difference which was calculated according to the depth of the melt pool. In the calculation, the flow and solidification behavior of SS was simulated first and after all the SS passed through the channel, the flowing behavior of Zr in the fuel support piece was simulated. Thus during the calculation, there would be only one molten material. As the assumption mentioned above, the reaction between SS and Zr was ignored and from the calculation results it can be seen that the amount of crust formed on the wall of the pieces is small in this study. Thus the reaction between SS and Zr would have a less influence on the whole process and to simplify the calculation, the reaction between SS and Zr was not considered. Other interactions such as interactions between the B4C and stainless steel was not considered as well and the neglect of all these interactions simplified the calculation and the further research will be done in the future. The mass of SS was 18.56 kg which was calculated based on the data of one fuel assembly of Fukushima Daiichi-Unit 1 Nuclear Reactor (Tanabe, 2011), and the mass of the melt Zr was 47.75 kg which is equal to the mass of cladding of each assembly (Tanabe, 2011). In order to perform sensitivity analysis of the flow and freezing behaviors of the melt, several cases had been calculated. For the first case, the melt was set up as SS with different initial temperature or diameter of particles. For the second case, the sensitivity analysis has been done with different material of the melt. The initial condition for this calculation is listed in Table 3. To make the results comparable, the variable-controlling approach was adopted and only one initial condition was different in the cases for comparison.
5.1. Computational domain and boundary condition
5.2. Results and discussion
In this calculation, considering the symmetry of the fuel support piece, only one fourth of the fuel support was simulated. The part of the fuel support piece that above the core plate was neglected to reduce calculation amount. As shown in Fig. 16a, a threedimensional particle configuration was constructed for analyzing the melt flowing behavior within the fuel support piece. The geometric parameters of the fuel support piece are listed in Table 2. Fig. 16b is the cross section view of the computation domain. As shown in Fig. 16b, the inlet boundary particles were added at the top of the channel to inject melt particles. As to boundary conditions, nonslip conditions were applied at the inner surface of the fuel support piece, and a free surface boundary was applied at the surface of the fluid. The temperature of particles of the fuel support piece was set to be 600 K. It is difficult to accurately determine the heat transfer condition of the outer surface of the fuel support piece due to the complex environment and uncertain accident scenario. In this study, the radiation heat transfer boundary condition was applied to the outer surface of the fuel support piece. The core of the BWR is made up of several materials. The sheath of the control blade, the tie plate of the fuel assembly
Fig. 17 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten SS and the fuel support piece equal to 1701 K and 600 K, respectively. In Fig. 17a, the particle size is 2.0 mm and in Fig. 17b, the particle is 2.5 mm. For the case with 2.0 mm particle size, the particle total number of this simulation is 243,465. This calculation was performed on the computer with Intel Xeon E52697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory and during the 3.55 s simulation time, the total calculation time is 32 h 28 min and 0.27 s (initial time step is 5 104 s). For the
Table 3 The initial condition for the calculation. Parameter
Melt Material
Material (melt/wall)
Stainless steel (melt)
Zr (melt)
Stainless steel (wall)
Initial temperature [K] Material melting point [K] Melt mass [kg] Diameter of particle [mm]
1701, 1710, 1750 1700 18.56 2.0, 2.5
2099, 2118 2098 47.75 2.0
600 1700 – 2.0, 2.5
432
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 17. Topology sequences of the fuel support piece (initial temperature of SS is 1701 K).
mass stopped increasing when all the molten SS passed through the orifice of the fuel support piece. The temperature distribution of the wall of the fuel support piece is shown in Fig. 19. Initially, the inner surface of the fuel support piece was heated to a higher temperature. Then the
20
15
Outflow mass(kg)
case with 2.5 mm particle size, the particle total number of this simulation is 149,524. During the same simulation time, the total calculation time is 26 h 55 min and 1.22 s with the same initial time step. The particle size (distance between adjacent particles) can affect the results calculated by MPS method. In present study, the sensitivity analysis of the particle size was done. In the cases, all the conditions were the same except the size of particles. Fig. 17 shows the topology sequences of the fuel support piece calculated by MPS with different particle sizes. In both cases, the result indicated that the molten SS particle near the inner surface of the fuel support piece gradually solidified and formed a crust layer on the inner surface of the fuel support piece. The green one represents the wall of the fuel support pieces, the blue one represents the fluid, the red one represents the particle which is freezing, the grey one represents the formed crust and the yellow one represents the inlet boundary particle. Finally, as the mass of the melt is small, most of the molten SS passed through the fuel support piece within a few seconds and a thin crust formed on the wall of the fuel support piece which would influence the flow and heat transfer of the zircaloy that came later. Fig. 18 shows the profiles of the outflow mass predicted by MPS with different particle sizes. The results show that similar results could be achieved by setting the particle size to be 2.0 mm and 2.5 mm, as shown in Fig. 18. The rate of the outflow mass decreased with the time due to the formation of the crust layer on the inner surface of the fuel support piece. Finally, the outflow
10
5
2mm 2.5mm 0
0
1
2
3
4
5
Time(s) Fig. 18. The profiles of outflow mass predicted by MPS (initial temperature of SS is 1701 K).
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 19. The temperature distribution of the fuel support piece (initial temperature of SS is 1701 K).
Fig. 20. Comparison of pressure for different particle sizes.
433
434
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 21. Comparison of temperature for different particle sizes.
heat gradually transferred from the inner surface to the outer part by conduction. On the outer surface of the whole device, temperature increased quickly at the entrance and exit. At the entrance, the melt was of high temperature which led to the high temperature of neighboring wall. The complex structure of the exit had a great influence on the flow and heat transfer of the melt, which cause the exit of the fuel support piece to become hot. With the time going by, the heat transferred to the surrounding wall and the temperature distribution became uniform. The comparison of pressure and temperature for different particle sizes in MPS calculation is shown in Figs. 20 and 21, and similar results were obtained. In the initial stage, because of the injection of the molten SS, the pressure of the fuel support piece was high and with the time going by, the amount of melt SS at the inlet decreased which led to reduction of the pressure at the upper part. At the end of the calculation, all the molten SS passed through the fuel support piece and the pressure reduced. The sensitivity analysis of initial temperature of the molten SS were performed by setting it as 1710 K and 1750 K to estimate the effect of the initial temperature on the melt behavior within the fuel support piece. Fig. 22a shows the topology sequences of the fuel support piece predicted by MPS with the SS initial temperature equals to 1710 K and Fig. 22b is for the case in which the initial temperature equals to 1750 K. As shown in the figures, a crust layer gradually formed on the
inner surface of the fuel support piece as the molten SS discharged through the fuel support piece. In both of the cases, the fuel support piece had not been plugged by the solidified SS particles when all the SS drained off. Fig. 23 shows the comparison among the outflow mass profiles of 1701 K case, 1710 K case and that of 1750 K case. As the amount of the molten SS was small and almost all the SS passed through the fuel support piece during the process of simulation, the melt initial temperature has a slight effect on the outflow mass of the three cases. As shown in the Fig. 23, the outflow mass increasing rate and the final outflow mass are almost the same. For the case of 1750 K SS, due to the higher temperature, it’s harder to form the curst on the wall and the final quality of outflow mass is a little bit more than the other two cases. With the temperature of the core increasing, the zircaloy in the assembly will be melted and fall into the fuel support piece. To simulate this process, the material of the melt was set up as molten zircaloy in the next cases. Considering the influence of the temperature field and crust that formed during the flow of molten SS, these zircaloy cases were calculated based on the results of previous SS case. Fig. 24 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten Zr and the fuel support piece equal to 2099 K and 600 K, respectively. The particle total number of this simulation is 243,465, and particle diameter is 0.002 m.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
435
Fig. 22. Topology sequences of the fuel support piece ((a) initial temperature of SS is 1710 K (b) initial temperature of SS is 1750 K).
This calculation were performed on the computer with Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) 64 GB memory. During the 22.13 s simulation time, the calculation time is 146 h 45 min and 31.23 s (initial time
Intel and total step
is 5 104 s). Fig. 24 indicated that the molten Zr particle near the inner surface of the fuel support piece gradually solidified and formed a crust layer on the inner surface of the fuel
20
SS 1701K SS 1710K SS 1750K
Outflow mass(kg)
15
10
5
0
0
1
2
3
4
Time(s) Fig. 23. The comparison between the SS outflow mass of different initial temperature.
support piece. Finally, the fuel support piece began to be plugged by the solidified Zr particles at the orifice before all the molten Zr drained off. The pressure distribution calculated by the MPS is shown in Fig. 25. At the beginning of the calculation, a lot of molten Zr accumulated at the top of fuel support piece which led to a high pressure and with the melt passing through, the amount of molten Zr at the inlet decreased and the pressure decreased at the upper part. Finally, due to the continuous heat transfer, the molten Zr changed into crust and the pressure reduced as well. Fig. 26 shows the temperature distribution of the wall of the fuel support piece during the simulation process. At first, the wall of the fuel support piece was heated by the hot melt. Then the heat transferred from the inner surface to the outer part by conduction and the wall gradually cooled down. The initial temperature of the molten Zr is uncertain when it discharges into the fuel support piece. In this study, another case with initial melt temperature of 2118 K was performed to estimate the effect of the initial temperature on the melt behavior within the fuel support piece. Fig. 27 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten Zr and the fuel support piece equal to 2118 K and 600 K, respectively. A crust layer gradually formed on the inner surface of the fuel support piece as the molten Zr discharged through the fuel support piece. However, in this case, the fuel support piece had not been plugged by the solidified Zr particles when all the Zr drained off.
436
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 24. Topology sequences of the fuel support piece (initial temperature of Zr is 2099 K).
Fig. 25. Pressure distribution of the fuel support piece (initial temperature of Zr is 2099 K).
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
437
Fig. 26. The temperature distribution of wall (initial temperature of Zr is 2099 K).
Fig. 28 shows the comparison between the outflow mass profile of 2099 K case and that of 2118 K case. Initially, the melt initial temperature has slight effect on the outflow mass increasing rate. After about four seconds, the outflow mass of
the 2099 K case increased more slowly than that of 2118 case, because for the lower initial temperature case the crust layer grew more quickly, in turn it obstruct the flow more significantly.
438
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig.27. Topology sequences of the fuel support piece (initial temperature of Zr is 2118 K).
50
Zr 2118K Zr 2099K
Outflow mass(kg)
40
30
20
10
0
0
5
10
15
2
Time(s) Fig. 28. The comparison between the Zr outflow mass of different initial temperature.
performed with stainless steel and zircaloy. In the calculation, the SS in the core of reactor was assumed to be melted first, and then the Zr was melted and flowed into the fuel support piece when SS drained out of it. Thus, the whole calculation was divided into two calculation cases. For the SS cases, the crust layer formed on the inner surface of the fuel support piece during the melt discharging process. However, due to the small amount of molten SS, the crust formed in the fuel support piece was thin and the fuel support piece had not been plugged before all the melt drained off. For the zircaloy cases, the crust layer also formed on the inner surface of the fuel support piece during the melt discharging process. For the lower initial temperature (2099 K) case, the fuel support piece was plugged by the solidified Zr particles, whereas for the 2118 K case, the fuel support was almost plugged when all the melt drained off. In both cases, the fuel support piece kept intact till the fuel support piece was plugged or all the melt drained off. In this study, all the results indicated that the molten SS and Zr was unlikely to destroy the fuel support piece during severe accident of BWR under the assumed order of melt injection. It should be noted that as a new method for numerical calculation, MPS method can easily go for most engineering calculation and is worth a further develop.
6. Conclusion Acknowledgments In the present study, the MPS method was adopted to analyze the flow and freezing behavior of the melt in the fuel support piece during the severe accident of BWR. Considering that the reactor core contains several types of material, the simulation was
A part of this study is the result of ‘‘Mechanistic study of melt behavior in lower RPV head” carried out under the Initiatives for Atomic Energy Basic and Generic Strategic Research by the
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Ministry of Education, Culture, Sports, Science and Technology of Japan. The professional guidance and help from Dr. Yoshiaki Oka, Emeritus Professor of the University of Tokyo, is highly appreciated. The MPS code of the present study was developed based on the original code, MPS-SW-Main-Ver2.0 that was kindly provided by S. Koshizuka and K. Shibata. References Berthoud G., Duret B. 1989. The freezing of molten fuel: reflections and new results. In: Proc. 4th Topical Meeting on Nuclear Reactor Thermal-hydraulics, pp. 675– 681. Bui, V.A., Dinh, T.N., 1996. Modeling of heat transfer in heated-generating liquid pools by an effective diffusivity-convectivity approach. In: Proceedings of 2nd European Thermal-Sciences Conference, Rome, Italy, pp. 1365–1372. Chen, R., Oka, Y., 2014. Numerical analysis of freezing controlled penetration behavior of the molten core debris in an instrument tube with MPS. Ann. Nucl. Energy 71, 322–332. Chen, R., Tian, W., Su, G.H., et al., 2011. Numerical investigation on coalescence of bubble pairs rising in a stagnant liquid. Chem. Eng. Sci. 66 (21), 5055–5063. Cheung F.B., Shiah S.W., Cho D.H., et al., 1992. Modeling of heat transfer in a horizontal heat-generating layer by an effective diffusivity approach. ASMEHTD192, pp. 55–62. Duan, G., Chen, B., 2015. Large Eddy Simulation by particle method coupled with Sub-Particle-Scale model and application to mixing layer flow. Appl. Math. Model. 39 (10–11), 3135–3149. Gauntt R.O., Humphries L.L., 1997. Final Result of the XR2-1 BWR Metallic Melt Relocation Experiment, NUREG/CR-6527 SAND97-1039. Gauntt R.O., Cash J.E., Cole R.K., et al., 2005. MELCOR computer code manual. Core (COR) package reference manuals, NUREG/CR-6119, 2(2). Gotoh H., Khayyer A., 2011. Method and device for determining interface particle used in particle method, and program for determining interface particle. United States Patent Application Publication. US 2011/0172948 A1. Gotoh, H., Khayyer, A., 2016. Current achievements and future perspectives for projection-based particle methods with applications in ocean engineering. J. Ocean Eng. Mar. Energy 2 (3), 251–278. Gotoh, H., Shibahara, T., Sakai, T., 2001. Sub-particle-scale turbulence model for the MPS method Lagrangian flow model for hydraulic engineering. Adv. Methods Comput. Fluid Dyn. 9–4, 339–347. Groudev, P., Atanasova, P., Chatterjee, B., et al., 2013. ASTEC investigations of severe core damage behaviour of VVER-1000incase of loss of coolant accident along with Station-Black-Out. Nucl. Eng. Des. 272, 237–244. Ikeda, H., Koshizuka, S., Oka, Y., et al., 2001. Numerical analysis of jet injection behavior for fuel-coolant interaction using particle method. J. Nucl. Sci. Technol. 38 (3), 174–182. Jeong, S.M., Nam, J.W., Hwang, S.C., et al., 2013. Numerical prediction of oil amount leaked from a damaged tank using two-dimensional moving particle simulation method. Ocean Eng. 69 (5), 70–78. Khayyer, A., Gotoh, H., 2009. Modified Moving Particle Semi-implicit methods for the prediction of 2D wave impact pressure. Coast. Eng. 56 (4), 419–440. Khayyer, A., Gotoh, H., 2011. Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230 (8), 3093–3118.
439
Khayyer, A., Gotoh, H., 2013. Enhancement of performance and stability of MPS mesh-free particle method for multiphase flows characterized by high density ratios. J. Comput. Phys. 242 (1), 211–233. Kondo, M., Koshizuka, S., 2011. Improvement of stability in moving particle semiimplicit method. Int. J. Numer. Meth. Fluids 65, 638–654. Kondo M., Koshizuka S., Suzuki K., et al., 2007. Surface tension model using interparticle force in particle method. In: Proceedings of FEDSM 2007 5thJoint ASME/JSME Fluids. Engineering Conference July 30 – August 2, SanDiego, California USA. Koshizuka, S., Oka, Y., 1996. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123, 421–434. Koshizuka, S., Nobe, A., Oka, Y., 1998. Numerical analysis of breaking waves using the moving particle semi-implicit method. Int. J. Numer. Meth. Fluids 26, 751– 769. Li, X., Tian, W., Chen, R., et al., 2013. Numerical simulation on single Taylor bubble rising in LBE using moving particle method. Nucl. Eng. Des. 256 (17), 227–234. Li, G., Liu, M., Wang, J., et al., 2016. MELCOR 2.1 analysis of melt behavior in a BWR lower head during LOCA and SBO accident. Ann. Nucl. Energy 90, 195–204. Morita, K., Zhang, S., Koshizuka, S., et al., 2011. Detailed analyses of keyphenomena in core disruptive accidents of sodium-cooled fast reactors by the COMPASS code. Nucl. Eng. Des. 241, 4672–4681. Ostensen R.W., 1973. Extended fuel motion study, AHI-RDP-18. Rahman, M.M., Hino, T., Morita, K., et al., 2007. Experimental investigation of molten metal freezing onto a structure. Exp. Therm. Fluid Sci. 32, 198–213. Sakai, N., Horie, H., Yanagisawa, H., et al., 2014. Validation of MAAP model enhancement for Fukushima Dai-chi accident analysis with Phenomena Identification and Ranking Table (PIRT). J. Nucl. Sci. Technol. 51 (7–8), 951–963. Shao, S., Ji, C., Graham, D.I., et al., 2006. Simulation of wave overtopping by an incompressible SPH model. Coast. Eng. 53 (9), 723–735. Shirakawa, N., Yamamoto, Y., Horie, H., et al., 2002. Analysis of flows around a BWR spacer by the two-fluid particle interaction method. J. Nucl. Sci. Technol. 39, 572–581. Song, D., Hitoshi, G., 2005. Turbulence particle models for tracking free surfaces. J. Hydraul. Res. 43 (3), 276–289. Souto-Iglesias, A., Macià, F., González, L., et al., 2013. On the consistency of MPS. Comput. Phys. Commun. 184 (3), 732–745. Tanabe, F., 2011. Analysis of core melt accident in Fukushima Daiichi-Unit 1 nuclear reactor. J. Nucl. Sci. Technol. 48, 1135–1139. Tian, W., Ishiwatari, Y., Ikejiri, S., et al., 2009. Numerical simulation on void bubble dynamics using moving particle semi-implicit method. Nucl. Eng. Des. 239 (11), 2382–2390. Tian, W., Ishiwatari, Y., Ikejiri, S., et al., 2010. Numerical computation of thermally controlled steam bubble condensation using Moving Particle Semi-implicit (MPS) method. Ann. Nucl. Energy 37 (1), 5–15. MAAP4 Users Manual, 1999. Fauske Associated Inc., 2. Wolf J.R., Rempe J.L. Stickler A., et al., 1994.TMI-2 vessel investigation project integration report (NUREG/CR-6197). Yang S.M., Zhang Z.Z., 1994. An experimental study of nature convection of natural convection heat transfer from a horizontal cylinder in high Rayleigh number laminar and turbulent regions. In: Proc. of the 10th International Heat Transfer Conference, Brighton, 7, pp. 185–189. Yoon, H.Y., Koshizuka, S., Oka, Y., 2001. Direct calculation of bubble growth, departure, and rise in nucleate pool boiling. Int. J. Multiph. Flow 27 (2), 277– 298.
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Numerical analysis of the melt behavior in a fuel support piece of the BWR by MPS Ronghua Chen a, Lie Chen a, Kailun Guo a, Akifumi Yamaji b, Masahiro Furuya c, Wenxi Tian a, G.H. Su a,⇑, Suizheng Qiu a a b c
School of Nuclear Science and Technology, Xi’an Jiaotong University, Xi’an 710049, China Cooperative Major in Nuclear Energy, Graduate School of Advanced Science and Engineering, Waseda University, Tokyo 169-8555, Japan Central Research Institute of Electric Power Industry (CRIEPI), 2-11-1 Iwado-kita, Komae, Tokyo 201-8511, Japan
a r t i c l e
i n f o
Article history: Received 23 June 2016 Received in revised form 5 January 2017 Accepted 8 January 2017 Available online 17 January 2017 Keywords: MPS method Fuel support piece Severe accident Freezing behavior
a b s t r a c t The fuel support piece in a boiling water reactor (BWR) is used to brace fuel assemblies. The channel within the fuel support piece is determined to be a potential corium relocation path from the core region to the lower head during the severe accident of BWR. In the present study, the improved ⁄⁄⁄Moving Particle Semi-implicit (MPS) method was adopted to simulate the flow and solidification behavior of the melt in a fuel support piece. The MPS method was first validated against the Pb-Bi plate ablation test that was performed by CRIEPI. The predicted ablation mass of the plate agreed well with the experimental results. Then the flowing and freezing behaviors of molten stainless steel (SS) and zircaloy in the fuel support piece were simulated by MPS method with a three dimensional particle configuration, respectively. In this study, the flow and solidification behavior of SS was simulated first. After all the SS passed through the channel, the flowing behavior of Zr in the fuel support piece was simulated. The simulation results indicated that the crust layer formed on the inner surface of the fuel support piece during the melt discharging process. The fuel support piece was plugged by the solidified zircaloy particles in the lower initial temperature case. The fuel support piece kept intact in all the calculation that were performed under the assumed order of melt injection. The present results could help to reveal the progression of a BWR severe accident. Ó 2017 Elsevier Ltd. All rights reserved.
1. Introduction Due to the battery power depletion and station blackout, few plant data such as pressure, temperature, and water level changes can be obtained during the accident at the Fukushima Dai-ichi Nuclear Power Station (Sakai et al., 2014). After the accident, local measurement is impossible due to high level radiation in the reactor building. Because of the above reasons, it is still unclear how the corium flows into the lower plenum and finally causes the vessel failure. During the late phase of sever accident in a boiling water reactor, a melt pool may form in the reactor core. There are five potential corium relocation paths from the core region to the lower head, as shown in Fig. 1, including the fuel support piece, control rod guide tube, instrumentation tube, the shroud wall and core support plate failure (Sakai et al., 2014). Because of the impossibility of local measurement, a numerical calculation and analysis is needed ⇑ Corresponding author. E-mail address: [email protected] (G.H. Su). http://dx.doi.org/10.1016/j.anucene.2017.01.007 0306-4549/Ó 2017 Elsevier Ltd. All rights reserved.
and the result will help to find out the progress of the severe accident occurred in the Fukushima Dai-ichi nuclear power station. During the flowing of high temperature corium in fuel support pieces, great thermal attack will be brought into these channels. The corium may ablate the wall of channel which leads the fuel support piece to be damaged, or it may freeze and plug the channel if the heat dissipation condition of corium is quite well. The calculation of the flowing and heat transfer of corium can assess the integrity of fuel support piece under the severe accident. On the other hand, the evaluation of the flowing and solidification behavior of corium in fuel support piece will help to reveal the relocation path of the corium in the core, which could add knowledge on the research of severe accident in BWR. Whether the molten core debris will lead to the fuel support pieces failure or freeze and plug these channels is a significant issue that need to be investigated. In respect of this issue, several relevant experiments have been performed. Berthoud and Duret (1989) designed GEYSER which is an upward melt injection experiment to research the flowing and solidification behavior of molten UO2 in a small tube. Sandia National Laboratories performed the XR2-1 test to study the
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
423
Nomenclature A C T W a d f ~ f ! g h k l0 m n p ! r r t ! u
surface area, m2 specific heat, J/(kgK) temperature, K weight function thermal diffusivity, m2/s dimension number solid fraction of particle surface tension per unit mass gravity acceleration, m/s2 enthalpy, J/kg curvature, m1/coefficient of thermal conductivity, W/ (mK) particle size, mm mass, kg particle number density pressure, Pa particle location vector particle distance, m time, s particle velocity vector
corium relocation pathways and processes in BWR core under the condition of severe accident. XR2-1 test results indicated that the molten corium pool had not formed in the core under the ‘‘dry core” condition (Gauntt and Humphries, 1997). Rahman et al. (2007) performed a series experiments to figure out the basic physical characteristic in freezing mode of molten metal and characteristic parameters of the freezing during penetration into the structure. Experimental results such as penetration length and width of frozen metal on stainless steel and brass structures were compared under the air and water coolant conditions. Apart from experiments of the melt flowing and freezing behavior, an inspection in the TMI-2 nozzle and instrument tube specimens have been performed (Wolf et al., 1994) and the results indicated that the molten core debris had traveled a long distance outside the reactor pressure vessel through the instrument tube, and the coolant in the annulus prevented the core melt from melting through the in-core housing tube and leaking into the containment (Wolf et al., 1994). For numerical simulation of the corium flowing and solidification behavior, several models and methods have been put forward.
Fig. 1. Corium flow paths from the core in BWR.
Greek letters heat source, W/kg e emissivity / radiation heat transfer / scalar variable k diffusion coefficient l viscosity, Pas q density, kg/m3 r surface tension/stefan-boltzmann constant m kinematic viscosity, mm2/s
U
Superscript/subscript a ambient e effective i; j particle No. l liquid phase s solid phase min minimum w wall 0 initial condition ⁄ temporal
Cheung et al. (1992) developed an effective diffusivity approach to simulate the heat transfer in a horizontal heat-generating layer. Bui and Dinh (1996) put forward the Effective Convectivity Conductivity Model (ECCM) to describe the heat transfer. Ostensen (1973) developed the bulk freezing model and the conduction layer model. These two models were used to predict the penetration length of the molten core debris. In addition to these mechanism models, some analyses based on severe accident integral code (e.g. ASTEC Groudev et al., 2013, MAAP (Users Manual, 1999) and MELCOR (Gauntt et al., 2005) have been done as well as some CFD analyses. Morita et al. (2011) used COMPASS code which is a kind of particle methods to simulate the GEYSER and THEFIS melt penetration experiments. Li et al. (2016) analyzed the melt behavior in a BWR lower head during LOCA and SBO accident with MELCOR 2.1. By introducing a potential lower head failure mechanism, the melt behavior in BWR lower head and leakage paths under postulated LOCA with complete loss of ECCS and SBO accident were sensitively analyzed. In the 1990s, Koshizuka and Oka proposed the Moving Particle Semi-implicit (MPS) method Koshizuka and Oka, 1996. In this kind of particle method, the kernel function and the Lagrange approach are used to trace the particle movements in the dispersed phase. It is suitable for the calculation of incompressible fluid. With the continuous improvement of this method, it has been successfully applied to various fields to solve the problems with lager deformation (Koshizuka et al., 1998). In recent years, the MPS method has a lot of development. The weak compressibility is adopted by Ikeda et al. (2001) to enhance the stability, and a higher order source term for the Pressure Poisson Equation (PPE) is deduced by Khayyer and Gotoh (2009). A higher order Laplacian model for discretizing the Poisson Pressure Equation and a corrective matrix for pressure gradient model were proposed to improve the stability and accuracy of the moving particle semi-implicit method by Khayyer and Gotoh (2011). Jeong et al. (2013) applied an appropriate gradient model and 4th-order Runge–Kutta time integration scheme to carry out numerical simulations for the Rayleigh–Taylor instability and found that the simulation results became closer to the experimental ones. The improved MPS method is also utilized in multiphase-flow, an Euler Lagrange method MPS-MAFL is introduced by Yoon et al., 2001, and Tian et al. (2009), Tian et al. (2010), Chen et al. (2011), Li et al. (2013) used the MPS-MAFL to do a lot of
424
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
research on bubble dynamics. An enhanced stabilized MPS method was developed by Khayyer and Gotoh (2013) for simulation of multiphase flows characterized by high density ratios. Shirakawa et al. (2002) analyzed flows around a BWR spacer with a twofluid particle interaction method (TF-MPS) and the result shows an agreement with the experiment which suggests a certain applicability of this method. Chen and Oka (2014) analyzed the melt behavior in a cold tube with MPS method in which the change of melt viscosity with temperature was taken into account. Beside these, the turbulence model is also adopted in MPS method, many researchers simulated the complex flow by MPS-LES (Gotoh et al., 2001; Song and Hitoshi, 2005; Shao et al., 2006; Duan and Chen, 2015). In this study, the improved MPS method was used to analyze the corium flowing and solidification behavior in a fuel support piece during a sever accident of BWR. At first, the models used in the improved MPS were validated against the analytical solutions and then the thin plate ablation test performed by CRIEPI was simulated and the feasibility of analyzing corium flowing and solidification behavior for this method was validated by comparing the calculation results with experimental results. Then the fuel support piece was simulated by a three-dimensional particle configuration and the behavior of corium in this structure was analyzed under different initial conditions. In this calculation, the SS was assumed to be melted first, and then the Zr was melted and flowed into the fuel support piece when all the SS passed through the channel. 2. The MPS method The basic theory of Moving Particle Semi-implicit (MPS) method is that the macroscopical fluid is represented by a series of discrete particles and the differential operator in governing equation is expressed by forms of particle interaction through corresponding computational model. 2.1. Governing equations and computational models As is often the case in problems of incompressible viscous fluid, the continuity, Navier–Stokes and energy equations are adopted to describe the flow characteristic or the flow field. The continuity equation is presented as follow,
Dq ¼0 Dt
ð1Þ
In MPS method, the Lagrange reference coordinates is used, therefore there is no convection term in momentum equation which is advantageous to the numerical calculation.
q
D~ l l þ q~f þ q~ g ¼ rP þ lr2 ~ Dt
ð2Þ
In the above momentum equation, from left to right, the terms are transient term, pressure gradient term, viscous term, surface tension term, and gravity term. In the calculation, the surface and the gravity terms were solved explicitly. In order to guarantee good convergence and accuracy, the pressure and the viscous terms were solved implicitly.
weight function (Kondo and Koshizuka, 2011) is expressed as follow,
(
WðrÞ ¼
ni ¼
X ! ! Wðj rj r i jÞ
ð5Þ
j–i !
!
r i and r j are the position vectors of particle i and j, respectively. For the incompressible fluid, the density is a constant. In the MPS method, fluid density is characterized by the particle number density, therefore the particle number density for incompressible fluid which is represented by n0 is constant as well.In 1996, Koshizuka and Oka developed the gradient model and Laplacian model (Koshizuka and Oka, 1996) for discretizing the differential operator in the conservation equations.
hruii ¼
! ! d X uj ui ! ! ðr j r i ÞWðrj r i Þ n0 j–i j r!j r!i j2
hr2 uii ¼
! ! 2d X ðu ui ÞWðr j r i Þ kn0 j–i j
ð6Þ
ð7Þ
To ensure the astringency of this method for high viscosity calculation, an implicit scheme was introduced to solve the viscous term. By using the Laplacian model, the viscous term was implicitly discretized as follow, ! u i
!
¼ u ni þ tDt
! ! 2d X ! ! ðu u i ÞWðj r nj r ni jÞ kn0 j–i j
ð8Þ
!
The temporal particle velocity, u ; can be obtained by solving the Eq. (8). On the other hand, to avoid the pressure oscillation during the calculation, the solving of Poisson equation of pressure in MPS method is improved (Kondo and Koshizuka, 2011). The source terms of Poisson equation is changed as,
1
q0
r2 P ¼
1 b n 2nn þ nn1 b c n nn c n n0 þ þ 2 2 2 0 0 n n n0 Dt Dt Dt ð9Þ
Where b and c are used to adjust the compensating parts and they should obey the relationship 0 6 c 6 b 6 1. A more proper source term of Poisson Pressure Equation was proposed by Khayyer and Gotoh (2011) with dynamic coefficients as functions of instantaneous flow field. The differential operator r2 P was discretized by using the Laplacian model mentioned above,
r2 P ¼
! ! 2d X ðP j Pi ÞWðr j r i Þ kn0 j–i
ð10Þ
and the final discretized Poisson equation is as follows,
X ! ! ðP i Pj ÞWðr j ri Þ
ð3Þ
Terms in the energy equation are transient term of enthalpy, conduction term and internal heat source term. In the theory of MPS method, the scope of action for one particle is limited by an effective radius re : Considering the influence of the distance from the central particle, a weight function is introduced into the calculation of this method. In this study, the adopted
ð4Þ
re P r
0;
A parameter named particle number density which is proportional to the fluid density is defined as,
j–i
@h ¼ kr2 T þ U @t
ð1 r=r e Þ2 ; 0 6 r < r e
¼
kq0 n0 1 b n 2nn þ nn1 b c n nn c n n0 þ þ 2 2 2 0 0 0 2d n n n Dt Dt Dt ð11Þ
The surface tension model in the program was developed by Kondo et al. (2007) which was based on the inter-particle potential force, !
f ¼ Kðr r min Þðr re Þ=m
ð12Þ
425
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
where K; r, r min are fitting coefficient, particle distance, and initial minimum particle distance, respectively. The potential force coefficients between two fluid particles and between a fluid particle and a solid particle are given as follows,
X
Kf ¼
1 3
2rr 2min
r 32 r min þ 12 re ðr re Þ2
ð13Þ
i2A;j2B;jr ij j
K fs ¼
1 K f ð1 þ cos hÞ 2
ð14Þ
Where h is the contact angle. In Kondo’s model, the surface tension coefficient is same as the surface energy. The surface energy is the potential energy which is stored on a unit surface. In other words, the required work on creating a unit surface is equivalent to the surface energy. As shown in Fig. 2, A and B are the two areas that are used to create two new fluid surfaces for the calculation of surface tension in Kondo’s model as presented in Kondo et al. (2007). In order to simulate the melting and solidification process, the heat transfer and phase change model was incorporated into the program. Based on the energy equation, the enthalpy of particles is calculated with the equations given as follows, kþ1
hi
k
¼ hi þ
2d X kij ðT kj T ki Þwðj~ rj ~ r i jÞ Dt þ U Dt n0 k j–i
ð15Þ Fig. 3. Calculation algorithm of MPS method.
where, the kij is the average coefficient of thermal conductivity,
kij ¼
2ki kj ki þ kj
ð16Þ
The temperature of every particle will be calculated based on its enthalpy.
8 T þ hhs0 > > < s qC ps T ¼ Ts > > : T þ hhs1 s
qC pl
ðh < hs0 Þ ðhs0 6 h 6 hs1 Þ
ð17Þ
ðhs1 < hÞ
T s is the melting point and hs0 , hs1 are the enthalpies when the particle begins and stops melting. The phase change model was introduced as follow, 0
f ¼
8 > <1
hs1 h > hs1 hs0
:
0
ðhs1 < hÞ
2.2. Calculation process The flow chart of MPS method is shown in Fig. 3. In one time step, the energy equation is solved firstly to obtain the temperature of particles. Based on the particle temperature, the viscosity of the particle is calculated. Then the viscous term of the Navier– Stokes equation is discretized as Eq. (8) and solved to obtain the !
temporal velocity u . Then gravity term and surface tension term !
are explicitly calculated to update theu . With the temporal veloc!
ðh < hs0 Þ ðhs0 6 h 6 hs1 Þ
0
f is the solid fraction of particle and different types of particle were introduced to characterize different phase states.
ð18Þ
ityu , the temporal particle position is calculated by the following equation. !
!
!
r ¼ r n þ Dtu
ð19Þ
!
After the r is obtained, the temporal particle number density n is solved by Eq. (5). Then, the particle pressure is updated by implicitly solving the Poisson equation of pressure with conjunction gradient method. And lastly, the particle velocity and position are updated using the new pressure gradient as follows, ! nþ1
u
! nþ1
r
!
¼ u !
¼ r
Dt
q0
rP
Dt 2
q0
rP
ð20Þ
ð21Þ
3. Validation of the MPS method
Fig. 2. The areas of A and B.
In this study, some changes had been made to the MPS method, including implicit viscous term calculation model and the stability improvement model. The improved MPS method was employed to simulate the classical dam break problem to validate its correctness. Fig. 4 shows the initial particle configuration of the dam break calculation. The particle total number of this simulation is
426
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 4. Particle configuration for dam break problem.
2164, and particle diameter is 0.01 m, the wall condition is set as no-slip boundary condition and a free surface boundary was applied at the surface of the fluid. When the number density of a particle was less than bn0 ; the particle was considered as a free surface particle. Here b was set to be 0.97 and n0 was initial particle number density for inner particles. For the free-surface particles, the Dirichlet boundary condition was used. As the incompressible fluid is not sensitive to absolute pressure, the pressure of the free surface particle was set to be zero. A more accurate method for detecting free-surface particles was proposed by Gotoh and Khayyer (2011) which includes two steps. The first step is to calculate the number density of every fluid particle and find the particle whose density is smaller than the predetermined threshold. The second step is to check symmetry of the particle which is chosen from the first step and the asymmetrical particle is set to be free surface particle. The details for this method is presented in the Gotoh’s reference. This calculation were performed on the computer with Intel Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory. For the original MPS, the initial time step is 104 s, and the time step was adjusted by the CFL condition during the calculation. For the 5.0 s simulation time, the total calculation time is 2 min and 42.19 s. The improved MPS method almost does not influence the computational efficiency. Fig. 5
shows the time histories of the pressure at the test particle, which was marked by circle in Fig. 4. They are predicted by the original MPS method and the improved MPS method with different values of b and c: As shown in Fig. 5, the pressure oscillation can be suppressed by the improved MPS method compared with the original MPS method. The differences of the pressure histories among the three parameter combinations are slight. Thus, the values of b and c were set as 0.5 and 0.4 respectively, in this study. To access the capability of the improved MPS to analyze the heat transfer and the phase change, the one-dimensional solidification problem was simulated by the improved MPS method. Fig. 6 shows the particle configuration of the one-dimensional solidification problem. The wall particle is on the left side with a constant temperature 173 K, the initial temperature of water is 283 K. The heat transfer boundary condition of the water is thermal insulating. The particle diameter shown in Fig. 6 is 0.002 m, and the sensitivity analysis of the particle size is also performed with different particle diameters including 0.001 m, 0.002 m and 0.0025 m. In this simulation, the pressure calculation is not involved. This simulation is performed on the computer with Intel Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory. When l0 equals 0.001 m, the particle total number is 20,000, initial time step is 103 s, the total calculation time is 29 min and 1.003 s for 4.0 s simulation time. When l0 equals 0.002 m, the particle total number is 5000, initial time step is 103 s, the total calculation time is 5 min and 4.244 s for 4.0 s simulation time. When l0 equals 0.0025 m, the particle total number is 3200, initial time step is 103 s, the total calculation time is 3 min and 31.227 s for 4.0 s simulation time. Due to release heat to the cold wall, the water gradually changes its phase from liquid to solid. The analytical solutions of the position of interface, temperature distribution of solid phase, and liquid phase are as follows, Position of interface
pffiffiffiffiffiffiffi nðtÞ ¼ 2R a1 t
ð22Þ
Temperature of solid phase:
T 1 ðt; xÞ ¼ T w þ ðT m T w Þ
erf ð pxffiffiffiffiffiÞ 2
a1 t
erf ðRÞ
Temperature of liquid phase:
pxffiffiffiffiffi 2 a2 t qffiffiffiffi T 2 ðt; xÞ ¼ T 1 þ ðT m T 1 Þ 1 erf R aa12
ð23Þ
1 erf
ð24Þ
Original MPS γ =0.4 β=0.5 γ =0.3 β=0.5 γ =0.4 β=0.6
Pressure/Pa
16000
12000
8000
4000
0
0
1
2
3
4
5
Time/s Fig. 5. Time history of the pressure at test particle.
Fig. 6. The particle configuration of the one-dimensional solidification problem.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439 300
280 273
Freezing point
Temperature/K
260
240
220
l0=0.001m
200
l0=0.0025m
l0=0.002m Analytical
180 0.00
0.05
0.10
0.15
0.20
x/m Fig.7. Results of solidification validation.
where, a1 , a2 are the thermal diffusivity of the wall and water. The result calculated with the program and the analytical solution at 4 s (y = 0.05 m) is presented in the Fig. 7.
427
As shown in the Fig. 7, the temperature distribution predicted by MPS agrees well with the analytical solution. With time going by, water with high temperature froze due to the cold wall at the left side. The solidification process of water and the temperature field of the case during the calculation is shown in Fig. 8. In Fig. 9, the temperature distribution at x = 0.1 m shows that the temperature near the free-surface is not consistent, this is because the number density of particle is less than n0 , thus the Laplacian model should be modified to compensate the PND (particle number density) near free-surface. It should be noted that the Laplacian model introduced by Eq. (7) cannot obtain accurate and consistent results for incomplete compact supports such as free surface. In 2013, Souto-Iglesias et al. improved the Laplacian model based on the analysis of the consistency of the MPS approximation to the Laplacian operator (Souto-Iglesias et al., 2013) and in 2016, Gotoh and Khayyer proposed a more consistent Laplacian model when they calculated the surface tension by considering BI (boundary integrals) (Gotoh and Khayyer, 2016). In this study, we carried out several calculations to evaluate the influence of this kind of Laplacian model. As shown in Fig. 9, we can conclude from Fig. 9 that the error on free-surface is decreasing when the l0 (particle size) is smaller. When l0 = 0.002 m, the relative error is 0.037%. For this reason, in the engineering simulation, the error on freesurface is ignored in this paper.
Fig. 8. Results calculated by MPS (l0 = 0.002 m), (a) solidification process, (b) temperature field.
428
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
l0=0.001m
274.4
l0=0.002m
274.2
l0=0.0025m Analytical
Temperature(K)
274.0 273.8
Table 1 Properties of materials in the thin plate ablation test. Properties 3
Density [kg/m ] Specific heat [J/kg K] Kinematic viscosity [m2/s] Heat conductivity [W/m K] Initial temperature [K]
273.6
Pb-Bi (44.5%, 55.5%)
Wood
Air
10050.0 147.0 1.29E07 13.8 300.0 (plate) Experimental value (inlet)
300.0 1300.0 – 0.069 300.0
1.205 1005 1.511E5 0.0257 300.0
273.4 273.2 273.0 272.8 0.00
0.02
0.04
0.06
0.08
0.10
y(m) Fig.9. Temperature distribution at x = 0.1 m.
4. Validation against the Pb-Bi plate ablation test 4.1. Introduction of the Pb-Bi plate ablation test Central Research Institute of Electric Power Industry (CRIEPI) in Japan had performed a series of plate ablation tests to add the knowledge on the RPV wall ablation mechanism, as shown in Fig. 10. In the plate ablation test, the Pb-Bi was employed as the substitute of the corium and was heated to 350 °C by the electric heaters. At the beginning of the test, the molten Pb-Bi was poured into a wooden vessel. Subsequently, the molten PbBi was discharged from the wooden vessel through the pipe that was connected to the bottom of the vessel. When the molten PbBi left the pipe, it dropt on the Pb-Bi plate. The thermocouples, located in both the wooden vessel and the pipe outlet, were used to measure the temperature of the molten Pb-Bi. Meanwhile, the ablation histories of the Pb-Bi plate were recorded by the video cameras from the side and the top during the experiment processes. The detailed parameters of material properties in this test is shown in Table 1.
Fig. 11. Computational domain of the Pb-Bi plate ablation test.
4.2. Numerical simulation by MPS method 4.2.1. Computational domain and boundary condition As shown in Fig. 11, a three-dimensional particle configuration was constructed for analyzing the plate ablation test. The inlet boundary particles were placed in the same position as the pipe outlet to inject melt particles. The angle between the melt injection direction and the plate is 5 degree. The angle between the plate and the horizontal plane is 15 degree. The plate is 5 mm in thickness, 17 mm in width and 37 mm in length. In this study, the particle diameter is 0.5 mm, and totally 52,626 particles are employed. At the end of the plate, the edge is chamfered which is the same as the experiment.
Fig. 10. Schematic of the plate ablation test facility.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
As to boundary conditions, nonslip conditions were applied at the surface of the Pb-Bi plate, and a free surface boundary was applied at the surface of the fluid. For the surface particles, nature convection heat transfer between particles and air and radiation heat transfer between the surfaces particle and the ambient were taken into consideration. The nature convection heat transfer (Yang and Zhang, 1994) was calculated as,
( Nu ¼
0:59ðPr GrÞ1=4 ;
for laminar flow
0:11ðPr GrÞ1=3 ;
for tubulent flow
ð25Þ
The radiation heat transfer between the surface and the ambient was calculated as,
/ ¼ eArðT 4w T 4a Þ
ð26Þ
where, e is the emissivity and r is the stefan-boltzmann constant. The inflow temperature was set according to the fitting curve based on the experimental data which was shown below,
8 0:0s < t 6 0:35s > < 403; TðKÞ ¼ 18:519t3 123:94t2 þ 293:54t þ 46:008; 0:35s < t 6 3:0s > : 0:0598t3 1:4404t 2 þ 10:578t þ 279; 3:0s < t 6 10:0s ð27Þ The inflow mass rate was approximately set to be 15.22 g/s.
429
4.2.2. Comparison between the numerical calculation and experimental result Fig. 12 shows the comparison of the ablation histories predicted by MPS and the experimental results. In the test, the plate side wall was melted through by the melt at about 4.0 s. The similar phenomenon of side wall melt through was successfully reproduced by MPS as show in Fig. 12. Because of the overshadowing of boundary particles, the melt-though of the sides of the plate couldn’t be presented in the figure. In the simulation of MPS, the edge of the plate was melted though at 4.0 s, and the melt inventory within the plate began to discharge. At 4.5 s, almost all the inventory melt flowed out. On the other hand, the time for the edge being melt through and all the inventory melt flowing out was little later in the experiment and there are two reasons that cause the difference between the simulation and the experiment. First, in the experiment, the transient mass flow rate at the inlet is hard to be accurately measured, thus it is hard to give the exactly same inlet condition in the simulation; Second, a film forms on the surface of the molten Pb-Bi when it is freezing, this film will influence the shape of the molten material. But in the simulation, limited by the computing capabilities, the influence of the film is not considered unless the size of the particle is small enough. Generally, the predicted plate ablation histories agreed with the experimental observation. Fig. 13 shows the temperature and pressure distribution during the simulation period.
Fig. 12. The plate ablation histories predicted by MPS and measured in experiment.
430
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 13. Temperature and pressure distribution of the ablation experiment.
0.14
Inflow mass MPS Outflow mass MPS Outflow mass Exp.
0.12
Outflow mass(kg)
0.10 0.08 0.06 0.04 0.02 0.00
0
1
2
3
4
5
6
7
Time(s) Fig. 14. The mass profiles predicted by MPS and measured in the test.
Fig. 14 shows the inflow mass profile and the outflow mass profile predicted by MPS, and the experimental measured result was also shown in Fig. 14 for comparison. Within three seconds, small difference was observed from the three profiles. After four seconds, the predicted out flow mass began to increase quickly, and finally the outflow mass is consistent with the experimental value. The comparison results indicated that the present MPS method is able to simulate the thin plate ablation test reasonably. Next, the MPS method was employed to analyze the flow and solidification behavior of the melt in the fuel support piece, which has the similar mechanism as the plate ablation test.
Fig. 15. Schematic of fuel support piece.
5. Analysis of the corium behavior in the fuel support piece Each fuel support piece carries four fuel assemblies as shown in Fig 15. Accordingly, at the bottom of the fuel support piece there are four orifices which allow the coolant to flow into the fuel support piece and then pass through it and entry into the fuel assemblies. In the BWR core melting down accident, passing though the fuel support piece is a possible way for the corium to
431
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 16. Computation domain of the fuel support piece.
Table 2 The main parameters of the size for the fuel support piece. Parameter
Value [cm]
Total height of computational domain Outer diameter of upper barrel Height of transition region Outer diameter of lower barrel Height of inclined section Outlet diameter of bottom outlet
9.2 6.0 0.8 4.6 5.4 1.8
relocate from the core region to the lower head. In this study, the MPS method was used to analyze the melt behavior in the fuel support piece.
and a few other parts are made by stainless steel. And the claddings of fuel are made up of zircaloy (Zr). During the severe accident, the control blade would be inserted into the fuel assemblies and the sheath which is made up of SS is possible to reach the core. As the melting point of SS is lower than that of the zircaloy, it may be melted first and flows into the fuel support pieces during severe accident of BWR and with the temperature of the core increase, other materials will be melted later. In the simulation, a melt pool was assumed to be formed within the fuel assembly. The melt was injected into the fuel support piece from the top driven by the pressure difference which was calculated according to the depth of the melt pool. In the calculation, the flow and solidification behavior of SS was simulated first and after all the SS passed through the channel, the flowing behavior of Zr in the fuel support piece was simulated. Thus during the calculation, there would be only one molten material. As the assumption mentioned above, the reaction between SS and Zr was ignored and from the calculation results it can be seen that the amount of crust formed on the wall of the pieces is small in this study. Thus the reaction between SS and Zr would have a less influence on the whole process and to simplify the calculation, the reaction between SS and Zr was not considered. Other interactions such as interactions between the B4C and stainless steel was not considered as well and the neglect of all these interactions simplified the calculation and the further research will be done in the future. The mass of SS was 18.56 kg which was calculated based on the data of one fuel assembly of Fukushima Daiichi-Unit 1 Nuclear Reactor (Tanabe, 2011), and the mass of the melt Zr was 47.75 kg which is equal to the mass of cladding of each assembly (Tanabe, 2011). In order to perform sensitivity analysis of the flow and freezing behaviors of the melt, several cases had been calculated. For the first case, the melt was set up as SS with different initial temperature or diameter of particles. For the second case, the sensitivity analysis has been done with different material of the melt. The initial condition for this calculation is listed in Table 3. To make the results comparable, the variable-controlling approach was adopted and only one initial condition was different in the cases for comparison.
5.1. Computational domain and boundary condition
5.2. Results and discussion
In this calculation, considering the symmetry of the fuel support piece, only one fourth of the fuel support was simulated. The part of the fuel support piece that above the core plate was neglected to reduce calculation amount. As shown in Fig. 16a, a threedimensional particle configuration was constructed for analyzing the melt flowing behavior within the fuel support piece. The geometric parameters of the fuel support piece are listed in Table 2. Fig. 16b is the cross section view of the computation domain. As shown in Fig. 16b, the inlet boundary particles were added at the top of the channel to inject melt particles. As to boundary conditions, nonslip conditions were applied at the inner surface of the fuel support piece, and a free surface boundary was applied at the surface of the fluid. The temperature of particles of the fuel support piece was set to be 600 K. It is difficult to accurately determine the heat transfer condition of the outer surface of the fuel support piece due to the complex environment and uncertain accident scenario. In this study, the radiation heat transfer boundary condition was applied to the outer surface of the fuel support piece. The core of the BWR is made up of several materials. The sheath of the control blade, the tie plate of the fuel assembly
Fig. 17 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten SS and the fuel support piece equal to 1701 K and 600 K, respectively. In Fig. 17a, the particle size is 2.0 mm and in Fig. 17b, the particle is 2.5 mm. For the case with 2.0 mm particle size, the particle total number of this simulation is 243,465. This calculation was performed on the computer with Intel Xeon E52697 v3 CPU (14 cores/28 threads, 2.60 GHz) and 64 GB memory and during the 3.55 s simulation time, the total calculation time is 32 h 28 min and 0.27 s (initial time step is 5 104 s). For the
Table 3 The initial condition for the calculation. Parameter
Melt Material
Material (melt/wall)
Stainless steel (melt)
Zr (melt)
Stainless steel (wall)
Initial temperature [K] Material melting point [K] Melt mass [kg] Diameter of particle [mm]
1701, 1710, 1750 1700 18.56 2.0, 2.5
2099, 2118 2098 47.75 2.0
600 1700 – 2.0, 2.5
432
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 17. Topology sequences of the fuel support piece (initial temperature of SS is 1701 K).
mass stopped increasing when all the molten SS passed through the orifice of the fuel support piece. The temperature distribution of the wall of the fuel support piece is shown in Fig. 19. Initially, the inner surface of the fuel support piece was heated to a higher temperature. Then the
20
15
Outflow mass(kg)
case with 2.5 mm particle size, the particle total number of this simulation is 149,524. During the same simulation time, the total calculation time is 26 h 55 min and 1.22 s with the same initial time step. The particle size (distance between adjacent particles) can affect the results calculated by MPS method. In present study, the sensitivity analysis of the particle size was done. In the cases, all the conditions were the same except the size of particles. Fig. 17 shows the topology sequences of the fuel support piece calculated by MPS with different particle sizes. In both cases, the result indicated that the molten SS particle near the inner surface of the fuel support piece gradually solidified and formed a crust layer on the inner surface of the fuel support piece. The green one represents the wall of the fuel support pieces, the blue one represents the fluid, the red one represents the particle which is freezing, the grey one represents the formed crust and the yellow one represents the inlet boundary particle. Finally, as the mass of the melt is small, most of the molten SS passed through the fuel support piece within a few seconds and a thin crust formed on the wall of the fuel support piece which would influence the flow and heat transfer of the zircaloy that came later. Fig. 18 shows the profiles of the outflow mass predicted by MPS with different particle sizes. The results show that similar results could be achieved by setting the particle size to be 2.0 mm and 2.5 mm, as shown in Fig. 18. The rate of the outflow mass decreased with the time due to the formation of the crust layer on the inner surface of the fuel support piece. Finally, the outflow
10
5
2mm 2.5mm 0
0
1
2
3
4
5
Time(s) Fig. 18. The profiles of outflow mass predicted by MPS (initial temperature of SS is 1701 K).
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 19. The temperature distribution of the fuel support piece (initial temperature of SS is 1701 K).
Fig. 20. Comparison of pressure for different particle sizes.
433
434
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 21. Comparison of temperature for different particle sizes.
heat gradually transferred from the inner surface to the outer part by conduction. On the outer surface of the whole device, temperature increased quickly at the entrance and exit. At the entrance, the melt was of high temperature which led to the high temperature of neighboring wall. The complex structure of the exit had a great influence on the flow and heat transfer of the melt, which cause the exit of the fuel support piece to become hot. With the time going by, the heat transferred to the surrounding wall and the temperature distribution became uniform. The comparison of pressure and temperature for different particle sizes in MPS calculation is shown in Figs. 20 and 21, and similar results were obtained. In the initial stage, because of the injection of the molten SS, the pressure of the fuel support piece was high and with the time going by, the amount of melt SS at the inlet decreased which led to reduction of the pressure at the upper part. At the end of the calculation, all the molten SS passed through the fuel support piece and the pressure reduced. The sensitivity analysis of initial temperature of the molten SS were performed by setting it as 1710 K and 1750 K to estimate the effect of the initial temperature on the melt behavior within the fuel support piece. Fig. 22a shows the topology sequences of the fuel support piece predicted by MPS with the SS initial temperature equals to 1710 K and Fig. 22b is for the case in which the initial temperature equals to 1750 K. As shown in the figures, a crust layer gradually formed on the
inner surface of the fuel support piece as the molten SS discharged through the fuel support piece. In both of the cases, the fuel support piece had not been plugged by the solidified SS particles when all the SS drained off. Fig. 23 shows the comparison among the outflow mass profiles of 1701 K case, 1710 K case and that of 1750 K case. As the amount of the molten SS was small and almost all the SS passed through the fuel support piece during the process of simulation, the melt initial temperature has a slight effect on the outflow mass of the three cases. As shown in the Fig. 23, the outflow mass increasing rate and the final outflow mass are almost the same. For the case of 1750 K SS, due to the higher temperature, it’s harder to form the curst on the wall and the final quality of outflow mass is a little bit more than the other two cases. With the temperature of the core increasing, the zircaloy in the assembly will be melted and fall into the fuel support piece. To simulate this process, the material of the melt was set up as molten zircaloy in the next cases. Considering the influence of the temperature field and crust that formed during the flow of molten SS, these zircaloy cases were calculated based on the results of previous SS case. Fig. 24 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten Zr and the fuel support piece equal to 2099 K and 600 K, respectively. The particle total number of this simulation is 243,465, and particle diameter is 0.002 m.
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
435
Fig. 22. Topology sequences of the fuel support piece ((a) initial temperature of SS is 1710 K (b) initial temperature of SS is 1750 K).
This calculation were performed on the computer with Xeon E5-2697 v3 CPU (14 cores/28 threads, 2.60 GHz) 64 GB memory. During the 22.13 s simulation time, the calculation time is 146 h 45 min and 31.23 s (initial time
Intel and total step
is 5 104 s). Fig. 24 indicated that the molten Zr particle near the inner surface of the fuel support piece gradually solidified and formed a crust layer on the inner surface of the fuel
20
SS 1701K SS 1710K SS 1750K
Outflow mass(kg)
15
10
5
0
0
1
2
3
4
Time(s) Fig. 23. The comparison between the SS outflow mass of different initial temperature.
support piece. Finally, the fuel support piece began to be plugged by the solidified Zr particles at the orifice before all the molten Zr drained off. The pressure distribution calculated by the MPS is shown in Fig. 25. At the beginning of the calculation, a lot of molten Zr accumulated at the top of fuel support piece which led to a high pressure and with the melt passing through, the amount of molten Zr at the inlet decreased and the pressure decreased at the upper part. Finally, due to the continuous heat transfer, the molten Zr changed into crust and the pressure reduced as well. Fig. 26 shows the temperature distribution of the wall of the fuel support piece during the simulation process. At first, the wall of the fuel support piece was heated by the hot melt. Then the heat transferred from the inner surface to the outer part by conduction and the wall gradually cooled down. The initial temperature of the molten Zr is uncertain when it discharges into the fuel support piece. In this study, another case with initial melt temperature of 2118 K was performed to estimate the effect of the initial temperature on the melt behavior within the fuel support piece. Fig. 27 shows the topology sequences of the fuel support piece predicted by the MPS method with the initial temperature of the molten Zr and the fuel support piece equal to 2118 K and 600 K, respectively. A crust layer gradually formed on the inner surface of the fuel support piece as the molten Zr discharged through the fuel support piece. However, in this case, the fuel support piece had not been plugged by the solidified Zr particles when all the Zr drained off.
436
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig. 24. Topology sequences of the fuel support piece (initial temperature of Zr is 2099 K).
Fig. 25. Pressure distribution of the fuel support piece (initial temperature of Zr is 2099 K).
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
437
Fig. 26. The temperature distribution of wall (initial temperature of Zr is 2099 K).
Fig. 28 shows the comparison between the outflow mass profile of 2099 K case and that of 2118 K case. Initially, the melt initial temperature has slight effect on the outflow mass increasing rate. After about four seconds, the outflow mass of
the 2099 K case increased more slowly than that of 2118 case, because for the lower initial temperature case the crust layer grew more quickly, in turn it obstruct the flow more significantly.
438
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Fig.27. Topology sequences of the fuel support piece (initial temperature of Zr is 2118 K).
50
Zr 2118K Zr 2099K
Outflow mass(kg)
40
30
20
10
0
0
5
10
15
2
Time(s) Fig. 28. The comparison between the Zr outflow mass of different initial temperature.
performed with stainless steel and zircaloy. In the calculation, the SS in the core of reactor was assumed to be melted first, and then the Zr was melted and flowed into the fuel support piece when SS drained out of it. Thus, the whole calculation was divided into two calculation cases. For the SS cases, the crust layer formed on the inner surface of the fuel support piece during the melt discharging process. However, due to the small amount of molten SS, the crust formed in the fuel support piece was thin and the fuel support piece had not been plugged before all the melt drained off. For the zircaloy cases, the crust layer also formed on the inner surface of the fuel support piece during the melt discharging process. For the lower initial temperature (2099 K) case, the fuel support piece was plugged by the solidified Zr particles, whereas for the 2118 K case, the fuel support was almost plugged when all the melt drained off. In both cases, the fuel support piece kept intact till the fuel support piece was plugged or all the melt drained off. In this study, all the results indicated that the molten SS and Zr was unlikely to destroy the fuel support piece during severe accident of BWR under the assumed order of melt injection. It should be noted that as a new method for numerical calculation, MPS method can easily go for most engineering calculation and is worth a further develop.
6. Conclusion Acknowledgments In the present study, the MPS method was adopted to analyze the flow and freezing behavior of the melt in the fuel support piece during the severe accident of BWR. Considering that the reactor core contains several types of material, the simulation was
A part of this study is the result of ‘‘Mechanistic study of melt behavior in lower RPV head” carried out under the Initiatives for Atomic Energy Basic and Generic Strategic Research by the
R. Chen et al. / Annals of Nuclear Energy 102 (2017) 422–439
Ministry of Education, Culture, Sports, Science and Technology of Japan. The professional guidance and help from Dr. Yoshiaki Oka, Emeritus Professor of the University of Tokyo, is highly appreciated. The MPS code of the present study was developed based on the original code, MPS-SW-Main-Ver2.0 that was kindly provided by S. Koshizuka and K. Shibata. References Berthoud G., Duret B. 1989. The freezing of molten fuel: reflections and new results. In: Proc. 4th Topical Meeting on Nuclear Reactor Thermal-hydraulics, pp. 675– 681. Bui, V.A., Dinh, T.N., 1996. Modeling of heat transfer in heated-generating liquid pools by an effective diffusivity-convectivity approach. In: Proceedings of 2nd European Thermal-Sciences Conference, Rome, Italy, pp. 1365–1372. Chen, R., Oka, Y., 2014. Numerical analysis of freezing controlled penetration behavior of the molten core debris in an instrument tube with MPS. Ann. Nucl. Energy 71, 322–332. Chen, R., Tian, W., Su, G.H., et al., 2011. Numerical investigation on coalescence of bubble pairs rising in a stagnant liquid. Chem. Eng. Sci. 66 (21), 5055–5063. Cheung F.B., Shiah S.W., Cho D.H., et al., 1992. Modeling of heat transfer in a horizontal heat-generating layer by an effective diffusivity approach. ASMEHTD192, pp. 55–62. Duan, G., Chen, B., 2015. Large Eddy Simulation by particle method coupled with Sub-Particle-Scale model and application to mixing layer flow. Appl. Math. Model. 39 (10–11), 3135–3149. Gauntt R.O., Humphries L.L., 1997. Final Result of the XR2-1 BWR Metallic Melt Relocation Experiment, NUREG/CR-6527 SAND97-1039. Gauntt R.O., Cash J.E., Cole R.K., et al., 2005. MELCOR computer code manual. Core (COR) package reference manuals, NUREG/CR-6119, 2(2). Gotoh H., Khayyer A., 2011. Method and device for determining interface particle used in particle method, and program for determining interface particle. United States Patent Application Publication. US 2011/0172948 A1. Gotoh, H., Khayyer, A., 2016. Current achievements and future perspectives for projection-based particle methods with applications in ocean engineering. J. Ocean Eng. Mar. Energy 2 (3), 251–278. Gotoh, H., Shibahara, T., Sakai, T., 2001. Sub-particle-scale turbulence model for the MPS method Lagrangian flow model for hydraulic engineering. Adv. Methods Comput. Fluid Dyn. 9–4, 339–347. Groudev, P., Atanasova, P., Chatterjee, B., et al., 2013. ASTEC investigations of severe core damage behaviour of VVER-1000incase of loss of coolant accident along with Station-Black-Out. Nucl. Eng. Des. 272, 237–244. Ikeda, H., Koshizuka, S., Oka, Y., et al., 2001. Numerical analysis of jet injection behavior for fuel-coolant interaction using particle method. J. Nucl. Sci. Technol. 38 (3), 174–182. Jeong, S.M., Nam, J.W., Hwang, S.C., et al., 2013. Numerical prediction of oil amount leaked from a damaged tank using two-dimensional moving particle simulation method. Ocean Eng. 69 (5), 70–78. Khayyer, A., Gotoh, H., 2009. Modified Moving Particle Semi-implicit methods for the prediction of 2D wave impact pressure. Coast. Eng. 56 (4), 419–440. Khayyer, A., Gotoh, H., 2011. Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230 (8), 3093–3118.
439
Khayyer, A., Gotoh, H., 2013. Enhancement of performance and stability of MPS mesh-free particle method for multiphase flows characterized by high density ratios. J. Comput. Phys. 242 (1), 211–233. Kondo, M., Koshizuka, S., 2011. Improvement of stability in moving particle semiimplicit method. Int. J. Numer. Meth. Fluids 65, 638–654. Kondo M., Koshizuka S., Suzuki K., et al., 2007. Surface tension model using interparticle force in particle method. In: Proceedings of FEDSM 2007 5thJoint ASME/JSME Fluids. Engineering Conference July 30 – August 2, SanDiego, California USA. Koshizuka, S., Oka, Y., 1996. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123, 421–434. Koshizuka, S., Nobe, A., Oka, Y., 1998. Numerical analysis of breaking waves using the moving particle semi-implicit method. Int. J. Numer. Meth. Fluids 26, 751– 769. Li, X., Tian, W., Chen, R., et al., 2013. Numerical simulation on single Taylor bubble rising in LBE using moving particle method. Nucl. Eng. Des. 256 (17), 227–234. Li, G., Liu, M., Wang, J., et al., 2016. MELCOR 2.1 analysis of melt behavior in a BWR lower head during LOCA and SBO accident. Ann. Nucl. Energy 90, 195–204. Morita, K., Zhang, S., Koshizuka, S., et al., 2011. Detailed analyses of keyphenomena in core disruptive accidents of sodium-cooled fast reactors by the COMPASS code. Nucl. Eng. Des. 241, 4672–4681. Ostensen R.W., 1973. Extended fuel motion study, AHI-RDP-18. Rahman, M.M., Hino, T., Morita, K., et al., 2007. Experimental investigation of molten metal freezing onto a structure. Exp. Therm. Fluid Sci. 32, 198–213. Sakai, N., Horie, H., Yanagisawa, H., et al., 2014. Validation of MAAP model enhancement for Fukushima Dai-chi accident analysis with Phenomena Identification and Ranking Table (PIRT). J. Nucl. Sci. Technol. 51 (7–8), 951–963. Shao, S., Ji, C., Graham, D.I., et al., 2006. Simulation of wave overtopping by an incompressible SPH model. Coast. Eng. 53 (9), 723–735. Shirakawa, N., Yamamoto, Y., Horie, H., et al., 2002. Analysis of flows around a BWR spacer by the two-fluid particle interaction method. J. Nucl. Sci. Technol. 39, 572–581. Song, D., Hitoshi, G., 2005. Turbulence particle models for tracking free surfaces. J. Hydraul. Res. 43 (3), 276–289. Souto-Iglesias, A., Macià, F., González, L., et al., 2013. On the consistency of MPS. Comput. Phys. Commun. 184 (3), 732–745. Tanabe, F., 2011. Analysis of core melt accident in Fukushima Daiichi-Unit 1 nuclear reactor. J. Nucl. Sci. Technol. 48, 1135–1139. Tian, W., Ishiwatari, Y., Ikejiri, S., et al., 2009. Numerical simulation on void bubble dynamics using moving particle semi-implicit method. Nucl. Eng. Des. 239 (11), 2382–2390. Tian, W., Ishiwatari, Y., Ikejiri, S., et al., 2010. Numerical computation of thermally controlled steam bubble condensation using Moving Particle Semi-implicit (MPS) method. Ann. Nucl. Energy 37 (1), 5–15. MAAP4 Users Manual, 1999. Fauske Associated Inc., 2. Wolf J.R., Rempe J.L. Stickler A., et al., 1994.TMI-2 vessel investigation project integration report (NUREG/CR-6197). Yang S.M., Zhang Z.Z., 1994. An experimental study of nature convection of natural convection heat transfer from a horizontal cylinder in high Rayleigh number laminar and turbulent regions. In: Proc. of the 10th International Heat Transfer Conference, Brighton, 7, pp. 185–189. Yoon, H.Y., Koshizuka, S., Oka, Y., 2001. Direct calculation of bubble growth, departure, and rise in nucleate pool boiling. Int. J. Multiph. Flow 27 (2), 277– 298.