2D MPS analysis of hydrodynamic fine fragmentation of melt drop with initial steam film during fuel–coolant interaction

Annals of Nuclear Energy 142 (2020) 107378

Contents lists available at ScienceDirect

Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene

2D MPS analysis of hydrodynamic fine fragmentation of melt drop with initial steam film during fuel–coolant interaction Gen Li, Panpan Wen, Haobo Feng, Jun Zhang, Junjie Yan ⇑ State Key Laboratory of Multiphase Flow in Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China

a r t i c l e

i n f o

Article history: Received 24 March 2019 Received in revised form 31 October 2019 Accepted 2 February 2020

Keywords: MPS method Fuel–coolant interaction Melt drop Fragmentation

a b s t r a c t The hydrodynamic fine fragmentation of melt drops, which increases the melt–coolant contact area and affects heat transfer and melt oxidation, is an important phenomenon during fuel–coolant interaction. The presence of a steam film surrounding the melt drops complicates the breakup simulation, including interface instabilities and high-density ratio. In this study, we modeled the breakup behavior of a melt drop with an initial steam film by using the meshless moving particle semi-implicit (MPS) method in the Lagrangian frame. Numerical models were validated step by step by simulating Kelvin–Helmholtz instability, Rayleigh–Taylor instability, and liquid drop breakup. The breakup behavior of a UO2 melt drop with and without initial steam film was simulated in 2D, and the fragment number, size distribution, and melt–coolant contact perimeter were comparatively analyzed. Results indicate that the initial steam film rapidly moved away from the melt drop and separated into two main bubbles, which had insignificant influence on either the melt drop deformation or the fragment number and the perimeter. The total number of fragments increased at the initial moments of the breakup process due to the flow inertia force and then decreased under the effect of surface tension. The number of small fragments (d/d0 = 0.01–0.03) increased and decreased in variations, but other relatively large fragments continuously increased in the entire duration. The perimeter contribution of small fragments was larger than those of other fragment groups in the initial duration of the melt drop breakup at We = 655, but it reversed subsequently. By contrast, the perimeter contribution by fragments with d/d0 = 0.01–0.02 and 0.02–0.03 was constantly larger than those of other groups at We = 2618. Moreover, the summations of perimeter contributions by all fragments in the present statistics were approximately 50% and 60%–80% for cases of We = 655 and 2618, respectively. Ó 2020 Elsevier Ltd. All rights reserved.

1. Introduction Fuel–coolant interaction (FCI) occurs at the in- or ex-vessel molten corium relocation stage during the severe accident of a nuclear reactor. Such contact mode can raise vigorous thermal interaction, intense steam generation, and, potentially, steam explosion. The FCI process can be divided into four phases in general: (i) initial premixing, (ii) triggering, (iii) detonation propagation, and (iv) hydrodynamic expansion (Haraldsson, 2000). Melt fragmentation mainly occurs in the premixing and detonation propagation phases. The premixing phase is characterized by melt jet breakup and melt drop fragmentation when the melt is poured into the coolant pool. The jets and drops break up into millimetric ligaments, droplets, and blobs, and the breakup is mainly driven by the velocity difference between the melt and the coolant surround⇑ Corresponding author. E-mail address: [email protected] (J. Yan). https://doi.org/10.1016/j.anucene.2020.107378 0306-4549/Ó 2020 Elsevier Ltd. All rights reserved.

ing it. During the detonation propagation phase, melt drops are rapidly fragmented into daughter droplets (i.e., secondary breakup) due to the shock wave at the drop surface. Melt drop fragmentation is an important phenomenon of FCI, which can rapidly increase the melt–coolant contact area and further affect heat transfer, melt oxidation, and solidification. If the fragments are cooled down, then their size distributions will influence the porosity and coolability of the debris bed. However, the behavior of melt drop breakup in coolant flow differs from a liquid drop breakup in gas flow that has been widely studied in sprays in medical and agricultural applications and combustion in engines (Hsiang and Faeth, 1992; Han and Tryggvason, 2001; Guildenbecher et al., 2009; Kékesi et al., 2014; Strotos et al., 2016). The analysis of melt drop breakup faces various complications, such as highly nonlinear interactions during FCI that involve steam films, interface hydrodynamic instabilities (Kelvin–Helmholtz [K-H] and Rayleigh–Taylor [R-T] instabilities), high density ratios, and heat and mass transfers.

2

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

The convoluted melt–vapor–coolant configuration during drop fragmentation was visualized experimentally by using a novel instrumentation system that integrates X-ray radiography and high-speed photography. Chen et al. (1999) investigated the fragmentation behavior of molten steel drops forced to explode in sustained pressure fields. When the steel drop was shocked by lowpressure pulse, all melts separated into submillimeter fragments, and nearly half of them manifested a size below 300 lm. By contrast, when the steel drop was shocked by high-pressure pulse, most of the melts were maintained but in the porous configuration. Chen et al. (1997) further compared the breakup scenario of a tin drop with a steel drop. The fragment sizes of tin and steel melts showed different distributions. The melt temperature exhibited a negligible effect on steel drop fragmentation, but the effect was significant on tin drops. Ciccarelli and Frost (1994) investigated the fragmentation process of a single molten metal drop immersed in water during vapor explosion. Their results on the low ambient velocity of water flow (<5 m/s) indicated that the cold drops (cerrolow alloy, T = 65 °C) broke up due to melt stripping in the relative flow, whereas the fragmentation of hot drops (tin, T = 700 °C) was governed by the growth and collapse of vapor bubbles at the drop surface. Hansson et al. (2009) visualized the vapor bubble dynamics and fine fragmentation process of a 0.5 g molten tin drop at 1000 °C during vapor explosion. Their analysis of simultaneous X-ray radiography and high-speed photography showed that instabilityinduced deformation and disintegration, such as swelling, occurred before droplet fine fragmentation. The experimental observations differed from the drop distortion and melt fingering phenomena observed by Ciccarelli, who conducted the experiment using an external pressure pulse trigger of two orders of magnitude larger. Increased pressure pulse enabled the melt drop surface to be further convoluted; subsequently, the explosive vaporization of the entrained coolant disintegrated the melt drop, leading to fine fragmentation. Park et al. (Park et al., 2005) observed a small-scale stratified explosion at the melt drop local region triggered by an external shock pulse of approximately 1.0 MPa; then, the explosion propagated along the melt surface, and thereby approximately 20% of the total mass of the melt drop was disintegrated. Uršic et al. (2014) investigated the hydrodynamic fine fragmentation process of partly solidified melt droplets during vapor explosion. A modified Weber number was proposed by replacing surface tension with crust stiffness as the stabilizing force. Their analysis indicated that crust constrained melt fine fragmentation and hence significantly influenced the intensity of vapor explosion. Haraldsson et al. (2001) evaluated the effect of surface solidification on melt drop fragmentation in liquid–liquid media. The experiments were performed with different melt materials, melt temperatures, and coolant temperatures. Results showed that the coolant temperature significantly affected melt fragmentation, and the fragmentation-controlled and freezing-controlled regimes were distinguished on the basis of different time scales of solidification and fragmentation processes. Several phenomenological modes to depict the melt fine fragmentation process have been proposed through experimental investigations and are categorized into thermal and hydrodynamic fragmentations. The thermal fragmentation mode is characterized by intense coolant evaporation that has penetrated melt drops at vapor film collapse; it mainly occurs in the triggering and early propagation phases of FCI. The hydrodynamic fragmentation mode is characterized by a strong shear flow at the melt drop surface. Interface instabilities incurred by relative flows and density gradients are the dominant factors of melt drop hydrodynamic deformation and fragmentation. Hydrodynamic fragmentation mainly occurs in the premixing, late propagation, and expansion phases. The time scale of hydrodynamic fragmentation is larger compared with that of thermal fragmentation, and it controls the global behavior of melt drop breakup.

The main determining parameter of droplet hydrodynamic deformation and fragmentation is the Weber number, which reflects the disturbance of relative motion between coolant flow and a melt drop. When the Weber number is less than a critical value, the melt drop simply vibrates in shape; otherwise, an irreversible drop disintegration occurs. Duan et al. (2003) numerically investigated this critical Weber number under different melt–coolant density ratios. Escobar et al. (2015) analysis showed that drop fragmentation in liquid–liquid configuration occurred at a Weber number of over 30, which was higher than that in liquid–gas configuration. Gelfand (1996) and Nomura et al. (2001) investigated drop breakup modes in conditions of Weber numbers higher than the critical value. Thakre et al. (2013) modeled the hydrodynamic deformation and the coalescence process of melt drops in a water pool by using the VOF method. Their results indicated that the initial Weber number significantly affected drop deformation, and the deformation rate increased with the rise of Weber number. Kim et al. (1983) experimentally investigated gallium drop fragmentation in water flow over a wide range of Weber numbers (30 – 3519). The gallium drop broke up into large daughter droplets at a low Weber number and into small fragments at a large Weber number; they formed fragment clouds when the Weber number further increased. Burger et al. (1984) calculated drop fragmentation behavior by using two different hydrodynamic modes, one based on Taylor instability and deformation breakup and the other one to describe fragmentation by stripping the capillary waves produced by the shear flow. The models can calculate the disintegrated mass and evaluate the actual size of the fragments. Past experimental studies have found thermal and hydrodynamic fragmentation modes in the interaction of the melt drop and the coolant and investigated the fragment morphology. Meanwhile, numerical simulations were performed in drop deformation modes, such as bag breakup, with different Weber numbers. However, little attention has been paid on fragment size distribution and melt–coolant contact areas, especially for melt drop with steam film. In recent years, meshless particle methods that have advantages in capturing fluid interfaces and free surface in complex boundary conditions have gained remarkable progress. Two representative methods are the smoothed particle hydrodynamics method initially developed using a fully explicit algorithm for compressible non-viscous flow and the moving particle semi-implicit (MPS) method developed for incompressible fluid flow. In the present study, we modeled the melt drop breakup behavior in a coolant flow by using the MPS method, which was previously applied in melt stratification and phase transition simulations (Li et al., 2016, 2017a, 2017b). However, the present MPS modeling faces problems of large density ratio and high Weber number. We validated numerical models step by step by simulating K-H instability, R-T instability, and liquid drop breakup. Then, the fragment number, size distribution, and melt–coolant contact boundary were comparatively analyzed with different Weber numbers. 2. Numerical method 2.1. Governing equations and models Mass and momentum conservation equations are the governing equations. The convection term is not included in the mass conservation equation but calculated by particle motion. Gravity was ignored in the present modeling, but surface tension term was included in the momentum equation. To model turbulence flow at the high Reynolds number, the large eddy simulation model was coupled with the MPS method.

Dq ¼0 Dt

ð1Þ

3

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378 

 Du 1  ¼  r P þmr2 u þr  s þ f s Dt q

ð2Þ

where q is density, u is velocity vector, t is time, P is pressure, m is kinematic viscosity, s is sub-particle scale turbulence stress, and fs is surface tension force. The overlines denote the variables in particle scale, which are directly solved from the governing equations. The turbulence flow caused by sub-particle scale vortexes is calculated by the sub-particle scale turbulence stress term, which is solved using the Smagorinsky eddy viscosity model. 



 

2 3

sab ¼ ua ub  ua ub ¼ 2mt Sab  kSPS dab 2





ð3Þ

1=2

mt ¼ ðC s DÞ ð2Sab Sab Þ 



Sab ¼

ð4Þ



1 @ua @ub ð þ Þ 2 @xb @xa

ð5Þ 

where vt is turbulence eddy viscosity, Sab is strain rate of mean flow, kSPS is turbulence kinetic energy, dab is Kronecker’s delta, D is filter width, and C s is a Smagorinsky constant. Turbulence kinetic energy at the sub-particle scale can approximate the turbulence eddy viscosity as follows: 1=2 mt ¼ C v kSPS D

ð6Þ

By combining Eqs. (5) and (6), turbulence kinetic energy can be expressed as

kSPS ¼

C 4s





D2 ð2Sab Sab Þ

2

Cm

ð7Þ

where C v is a sub-particle scale turbulence constant. Particle motion in the MPS method is implemented by calculating particle interaction models. The contribution of neighboring particles to the target particle is evaluated using a weight function. The differential operators, such as gradient, Laplacian, and divergence, in the governing equations are often discretized with the following equations:

hr/ii ¼ D

d X /j - /i ðrj - ri Þwðjrj - ri jÞ n0 j–i jrj - ri j2

E 2d X r / ¼ 0 ð/ - /i Þwðjrj - ri jÞ i kn j–i j 2

hr  uii ¼

ðuj  ui Þðrj  ri Þ 1 X re ni j–i jrj  ri j jrj  ri j2

1

q





rP

i

(  )    d X ðPj  Pi Þðrj  ri Þ d   ¼ 0 C wð r  r Þ þ 0 2 ij j i q þq  n j–i n ð i 2 j Þrj  ri  8  9 0 X <ðPi  P i;min Þðrj  ri Þ  =  Cij wðrj  ri Þ   : ; q rj  ri 2 j–i i

0

Cij . P i;min in the particle stabilizing term is the minimal pressure searched among the neighboring particles with the same type as particle i. The density qi in the denominator represents the fluid density of particle i instead of the average density; this treatment can effectively avoid the computational instability to be incurred by the high fluid density ratio. In addition, the correction matrix Cij can improve pressure gradient accuracy and is written in 2D as follows:

0 P w x2 P wij xij yij ij ij V i r2 V ij B j–i i r2ij j–i Cij ¼ B @ P wij yij xij P wij y2ij V i r2 V i r2

C C A

ð12Þ

ij

j–i

xij ¼ xj  xi , yij ¼ yj  yi , r2ij ¼ x2ij þ y2ij , wij ¼ wðrij Þ, and P V i ¼ d= wij . Cij is approximate to a unit matrix when the particle j–i

distribution is uniform and it has no effect on computational stability. Weight function is a function of particle distance and effective radius; its value decreases as the particle distance increases in the confined region.

 re wðrÞ ¼

r

 1 0 6 r < re

ð13Þ

re 6 r

0

where r is the distance between two particles, and re is the cutoff radius. In the present study, the same cutoff radius re = 3.1 l0 was used for the particle interaction models. The particle number density at position ri is defined as the summation of weight function calculated between particle i and its neighboring particles.

X

ni ¼

wðjrj - ri jÞ

ð14Þ

j–i

The Laplacian model coefficient is determined to ensure that the increase in variance in the Laplacian operator is equal to that of the analytical solution.

P j–i

ð9Þ

wðjrj - ri jÞjrj - ri j2 P wðjrj - ri jÞ

ð15Þ

j–i

Surface tension is an important influencing factor in melt drop deformation and fragmentation. A surface tension model based on pairwise potential, proposed by Kondo et al. (Kondo et al., 2007), was applied in the present study. The potential P(r) is a defined function that enables the potential force (the differential of potential function with respect to particle coordinate ri) to be long-range attractive and short-range repulsive.

(

P ðr Þ ¼

1 C 3



 r  32 l0 þ 12 re ðr  r e Þ2

ðr < r e Þ ðr P re Þ

0

ð16Þ

where C is a fitting coefficient calculated from the fluid surface tension, and l0 is the initial particle distance. The potential force fs between the two particles is calculated as

 fs ¼

ð11Þ

11

where

k¼

ð10Þ

ij

j–i

ð8Þ

where / and u represent a vector and a scalar, respectively; w(r) is weight function, n0 is the standard value of particle number density computed with an initial particle configuration; k is the Laplacian model coefficient, d is the number of spatial dimensions; r is the particle position vector; and the subscripts i and j denote particles i and j, respectively. However, when solving the pressure gradient term, an improved gradient model was adopted in present application to achieve high stability and accuracy (Duan et al., 2017).



Different from the standard gradient model (Eq. (8)), the present pressure gradient model includes a particle stabilizing term (second term of Eq. (11)) and a dimensionless correction matrix

Cðr  r e Þðr  l0 Þ

ðr < re Þ

0

ðr P re Þ

ð17Þ

which indicates that particle potential force is repulsive and attractive when the particle distance is lesser and larger than l0, respectively. When particle distance exceeds the effective radius, the potential force becomes sufficiently small and has to be set to zero to limit the computational expense.

4

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

The surface tension coefficient r is the potential energy stored on a unit fluid surface. The work required to create a unit surface is equivalent to potential energy. Therefore, when creating two new 2

fluid surfaces with the same area of l0 by detaching fluid particles that belong to two different regions, as shown in Fig. 1, the required work can be estimated by

X

  2 Pðri  rj Þ ¼ 2rl0

ð18Þ

i2A;j2B

where r is the surface tension coefficient. The particle potential in Eq. (18) is replaced by Eqs. (16), and (18) is reformed. Coefficient C can be calculated as 2

C¼

2rl0

P i2A;j2B;jr ij j
1 ðr 3



3 l 2 0

ð19Þ

þ 12 re Þðr  re Þ2

Parameter C is a constant coefficient, which is calculated at the initialization step rather than recalculated at each time step. 2.2. Computational algorithm

where c is an adjusting coefficient lying between 0.01 and 0.05, and a is an artificial compressibility coefficient. Coefficient a can effectively enhance computational stability and efficiency, but using this coefficient can result in the loss of fluid incompressibility. In practical applications, the value of a is usually set from 107 to 109 to avoid unnatural fluid deformation. The source term of pressure Poisson equation (right side of Eq. (22)) is a hybrid form, including the velocity divergence, particle number density difference and artificial compressibility terms. This hybrid source term can result in much smoother pressure distribution and higher computation stability, compared with the pure particle number density difference term (Koshizuka and Oka, 1996). In the solving process, the artificial compressibility term is moved to the left side and combined with the pressure Laplacian model. With pressure distribution, particle velocity and position are updated by solving the pressure gradient term.   1  ukþ1 ¼ u þ Dtð r PÞ

ð23Þ

1  rkþ1 ¼ r þ Dt 2 ð r PÞ

ð24Þ

q

q

The MPS method uses a semi-implicit computational algorithm to solve mass and momentum conservation equations. In the first step of one time step, the surface tension, sub-particle scale stress and viscous terms are explicitly calculated, and the temporal particle velocity and position are obtained as follows: 





u ¼ uk þ ðmr2 u þr  s þ f s ÞDt

ð20Þ



r ¼ rk þ Dtu

ð21Þ

where superscripts k and * denote the computation time step and temporal sub-step in one time step, respectively. After particle temporal motion, the particle distribution does not satisfy the fluid incompressibility and the mass conservation. Thus, in the second step, a pressure Poisson equation is deduced from the pressure term of the momentum equation in combination with the mass conservation equation (Duan et al., 2017; Khayyer and Gotoh, 2011).



r



1

q



kþ1

rP

¼ ð1  cÞ i

þa



 1 1 n  n0 r  u  c 2 0 Dt n Dt

1 kþ1 P Dt 2 i

ð22Þ

2.3. Validation of numerical method As described in Section 1, interface instabilities are the dominant factors in melt drop hydrodynamic deformation and fragmentation. To validate the numerical method, we thus simulated K-H and R-T instability phenomena, and compared the results with the data from the literature. Fig. 2 shows the growth of spirals at the interface of brine and fresh water in relative motion. The density difference between the fluids is 15.6 kg/m3, and the fluid layer height is 3.0 cm. The fresh water at the upper part moved to the right, whereas the brine moved to left. The density difference and the relative motion resulted in K-H instability that manifested as the rolling up of spirals. The MPS simulation has efficiently reproduced the spirals at the fluid interface in comparison with the experimental results. R-T instability is caused by the density gradient at the fluid interface in which a fluid penetrates another fluid. A classic example is the instability that occurs when a heavy liquid is placed over a light liquid, as shown in Fig. 3. The densities of heavy and light fluids are 1.225 and 0.1694 kg/m3, respectively. The viscosity of both fluids is 0.00313 m2/s. The surface tension coefficients are 1.0  103 and 2.0  103N/m for heavy and light fluids, respectively. A 1 m  4 m tank was discretized with a particle arrangement of 100  400, and the free-slip boundary condition was applied. The simulated mushroom-shaped spikes were similar to the interfaces obtained using the front-tracking method. The main front positions of the heavy liquid also agreed well. In addition to surface instability simulation, we simulated the water drop breakup in gas flow and compared it with experimental observations in Ref. (Li et al., 2019). We further simulated gallium drop breakup in water flow, i.e., liquid–liquid configuration. Considering that no experimental result was available for comparison, we compared the MPS simulation with the mesh method results, as shown in Fig. 4. Water flowed into the computation domain from the left with a velocity of 10.6 m/s, and the gallium drop diameter was 2.8 mm. The Weber and Reynolds numbers calculated by Eqs. (25) and (26) were 450 and 29775, respectively.

We ¼ Re ¼ Fig. 1. The work needed to create two new surfaces with area

2 l0 .

qc u20 d0 r

qc u0 d0 lc

ð25Þ ð26Þ

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

5

Fig. 2. K-H instability (upper part: experiment (Thorpe, 1968); lower part: simulation).

where u0 is water flow velocity, d0 is drop diameter, and the subscript c denotes the parameters of coolant, which means water flow in this case. Under the action of relative velocity, the gallium drop experienced deformation and fragmentation, and the daughter droplets formed and aggregated behind the mother drop. In general, the MPS simulation agreed well with the mesh method results.

premixing phase, the melt drops that separated from the melt jets continuously move in the coolant under the drive of inertia. By contrast, in the steam explosion and pressure wave propagation phases, the coolant is driven by pressure pulse and move rapidly. Either the former or the latter form is the relative motion between the melt drops and the coolant; thus, their fragmentation mechanisms are the same. Fig. 5 shows the schematic of the melt drop breakup. 2D configuration was considered because 3D simulation requires a huge amount of computation cost, especially for the present study. Past studies reported that the melt debris size should be approximately 0.1 – 10 mm; thus, we chose a UO2 drop with a diameter of d0 = 10 mm as the object. The UO2 drop was wrapped by an initial steam film that had a thickness of 1 mm and 2 mm accordingly in the analysis cases. The modeling of such initial steam film could be applied to the conditions of low superheat melt drops and the small melt drops, which have no further steam generation. The melt drop was impulsively accelerated by water flow with a velocity of u0. Water flowed in at the right boundary and out at the left boundary. The width and length of the computational domain were 6d0 and 14d0, respectively; the top and bottom were free-slip boundaries. Table 1 presents the calculation cases with determined velocities that enabled the Weber numbers to be larger than the breakup critical value. Gravity was ignored in the present calculation. The density ratio of UO2 melt to water is 8.77. The viscosities of UO2 melt and water are 0.0263 and 0.001 Pas, respectively. The surface tension coefficients of melt and water when exposed to atmosphere are 0.55 and 0.07 N/m, respectively. For the interface particles during the interaction with particles of different types, the harmonic means of their viscosity and surface tension coefficient were adopted accordingly. Particle number independence should be proven before setting the initial particle distance. Fig. 6 shows the fragment–coolant contact perimeter simulated with three different particle numbers at We = 1818. The negligible difference indicated that the selected particle configurations were all applicable to the present study. Moreover, considering the computation cost and to present the miniscule fragments with high resolution, we set the initial particle distance to 0.1 mm and the total particle number to 866621.

3. Analysis of melt drop fine fragmentation

3.2. Melt drop breakup behavior

3.1. Analysis conditions

Figs. 7 and 8 present the melt drop breakup behavior at We = 655 and 1818, respectively. An investigation of the cases with initial steam film showed that the steam film rapidly moved away from the melt drop under the effect of shear flow. The shear forces

Fig. 3. R-T instability (lines represent the interfaces simulated by the front-tracking method (Popinet and Zaleski, 1999).

Hydrodynamic fragmentation of melt drops mainly occurs in the premixing and propagation phases of steam explosion. In the

6

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

Fig. 4. Gallium drop breakup in the water flow at We = 450 (the top part is the simulation by using the Gerris flow solver proposed by Popinet (Escobar et al., 2014; Popinet, 2009; Popinet, 2003); the bottom part is the present simulation).

Fig. 5. Schematic of UO2 drop breakup.

acting on the melt drop and steam by the coolant flow were equivalent, but the consequent velocities significantly differed due to the large variances in their densities, as shown in Fig. 9a. Collapsed steam films aggregated at the low-pressure region behind the melt drop and then separated into two bubbles as the formation of vortexes, as shown in Fig. 9b. Since coolant evaporation was not considered in the present numerical modeling, the coolant directly interacted with the melt without additional steam generation. The interfacial instabilities presented at the coolant and melt interfaces due to the density and velocity differences. Consequently, part of the melt was stripped away and separated into many daughter droplets. The separation time of the initial steam was short at approximately 1–2 s in the present investigated cases. The deformation of the melt drop was insignificant when the initial steam film completely separated from the melt drop. Therefore, in the succeeding process, the melt drop deformed and fragmented in a similar manner as that without the steam film. The melt drop

Fig. 6. Particle number independence validation (fragment–coolant contact perimeter simulated with three different particle numbers).

was squeezed by the upstream and downstream flows into the shape of a crescent and then expanded upward and downward. Finally, the crescent-shaped melt disintegrated into many daughter ligaments and droplets. The fragment number was large at the high Weber number. Moreover, the time of the melt drop to complete the disintegration was approximately 20 s at We = 655 and approximately 12 s at We = 1818.

Table 1 Calculation conditions. Case

d0/mm

u0/ms1

qm/qc

We

Re

I II III IV V

10 10 10 10 10

6 8 10 12 14

8.77 8.77 8.77 8.77 8.77

655 1164 1818 2618 3564

6  104 8  104 10  104 12  104 14  104

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

Fig. 7. Melt drop breakup behavior at We = 655.

7

8

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

Fig. 8. Melt drop breakup behavior at We = 1818.

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

9

Fig. 9. Schematic of the steam separation mechanism.

3.3. Melt fragment number and size distribution Fig. 10 shows the total number of fragments at We = 655 and 2618. The statistics of the fragment number was performed using

Fig. 10. Fragment number of melt drops with and without initial steam film.

the image processing software Image-Pro Plus 6.0. The summation of all melt droplet areas was sustained the same as the initial melt drop area during post-processing to obtain an accurate fragment number. The fragment number increased at the initial moments because melt was stripped away from mother drop, thus forming many daughter droplets, as presented in Figs. 7 and 8. The fragment swarm moved to the downstream area, and the relative velocity between the melt droplets and the coolant flow decreased. Consequently, the effect of surface tension surpassed the initial flow inertia force, the melt fragments started to coalesce, and the number decreased. Therefore, the melt drop breakup process could be divided into two phases: initial deformation and fragmentation phase governed by flow inertial force and the subsequent fragment coalescence phase under the effect of fluid surface tension. The fragment numbers of melt drops with and without initial steam film had similar variation tendencies. Fig. 11 presents a comparison of the fragment numbers with different Weber numbers for the melt drop with 2 mm-thick steam film. The number increased as the Weber number increased, ascribing to the intensity of shear flow. Fig. 12 shows the fragment size distribution of the melt drop with 2 mm of initial steam film. The size was normalized with the fragment equivalent diameter and the initial melt drop diameter as d/d0. A large proportion of the fragments had a diameter of 0.01–0.02, accounting for nearly 70% of the total number. The number of small fragments presented increase and decrease variations

Fig. 11. Fragment number of melt drop with 2 mm initial steam film.

10

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

Fig. 12. Fragment size distribution of melt drop with 2 mm initial steam film.

(represented by filled symbols), whereas the relatively large fragments continuously increased in the entire duration (marked by empty symbols). The increase of large fragments was caused by the breakup of melt arms and the coalescence of daughter droplets, and these phenomena coincided with the image descriptions in Figs. 7 and 8. A comparison of size distributions at We = 655 and 2618 showed that the number of fragments with sizes of 0.01– 0.02 and 0.02–0.03 had a ratio of 4–5, but the numbers of large fragments (d/d0 > 0.05) were equivalent. 3.4. Fragment–coolant contact perimeter Melt–coolant contact area is an important parameter in the evaluation of heat transfer and melt oxidation. Here, the fragment surface area cannot be directly calculated because the present simulation was conducted in 2D. Instead, we used the increment of the normalized fragment–coolant contact perimeter to determine the contact area variation (Fig. 13) because the area usually has an order of magnitude of squared perimeter. The perimeter increased at the initial period and then tended to decrease. The cause of such tendency is similar to that of fragment number variation. The

Fig. 13. Fragment–coolant contact perimeter of melt drops with and without initial steam film.

change in perimeter further proved that no significant difference existed between melt drops with and without initial steam film. Fig. 14 presents a comparison of the fragment–coolant contact perimeters of the melt drop with 2 mm initial steam film and with different Weber numbers. The perimeter increased as the Weber number increased, and the maximum value at We = 3564 was approximately 23 times that of the initial melt drop perimeter. Hence, the melt–coolant contact area could be 529 times larger than the initial melt drop surface area. The maximum perimeter value at We = 3564 was approximately 1.9 times larger than that at We = 655. Fig. 15 shows the contributions of the fragments with different sizes to the melt–coolant contact perimeter. Only fragments with the size of d/d0 < 0.1 were compared to correspond with the preceding statistics of size distribution, whereas large fragments, including the mother droplet, were not considered. For We = 655, the contribution of fragments with d/d0 = 0.01–0.03 was larger than that of other fragments in the first 8 ms, but it decreased and was surpassed by fragments with d/d0 > 0.03. At

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

11

Fig. 14. Fragment–coolant contact perimeter of melt drops with 2 mm initial steam film.

the period of 10 to 20 ms, the largest perimeter contribution was from the fragments with d/d0 = 0.04–0.08 instead of the smaller or larger fragments. The case of We = 2618 yielded a different result where the perimeter contributions by fragments with d/d0 = 0.01 – 0.02 and 0.02 – 0.03 were constantly larger than those by other groups. The contributions of fragment groups marked with empty symbols in Fig. 15b increased with time. In addition, the summations of perimeter contributions of all fragments in the present statistics were approximately 50% and 60% – 80% for We = 655 and 2618, respectively. 4. Conclusions We modeled the breakup behavior of a melt drop with an initial steam film by using the MPS method. The numerical models were validated step by step by simulating K-H instability, R-T instability, and liquid drop breakup. Then, the breakup behavior of a UO2 melt drop with and without initial steam film was comparatively analyzed, and fragment number, size distribution, and fragment–coolant contact perimeter were investigated. The specific conclusions were drawn as follows: 1) No significant difference was observed in the melt drop deformation, fragment number and contact perimeter between the melt drops with and without initial steam film. 2) The initial steam film rapidly moved away from the melt drop under the effect of shear flow and separated into two main bubbles distributed at the central regions of the coolant flow vortexes. 3) The total number of fragments increased at the initial moments of the breakup process due to the flow inertia force and then decreased under the effect of surface tension. The number of small fragments, i.e., d/d0 = 0.01 – 0.03 in the present study, increased and decreased in variations, but the number of other relatively large fragments continuously increased in the entire duration. 4) The total melt fragment–coolant contact perimeter also presented increase and decrease variations with time. In the initial duration of the melt drop breakup at We = 655, the perimeter contribution by fragments with d/d0 = 0.01 – 0.03 was larger than that with d/d0 > 0.03, but it reversed subsequently. By contrast, the perimeter contribution by fragments with d/d0 = 0.01 – 0.02 and 0.02 – 0.03 at

Fig. 15. Fragment–coolant contact perimeter as contributed by different sized fragments.

We = 2618 was constantly larger than that by other groups. In addition, the summations of perimeter contributions of all fragments in the present statistics were approximately 50% and 60% – 80% for We = 655 and 2618, respectively. In future research, 3D simulation may be preferable in further investigating the melt–coolant contact area and in understanding the mechanisms of the physical phenomena during melt drop breakup. However, a huge challenge of 3D simulation is the computation cost. The estimated particle number is approximately 500 million when the uniform particle size adopted in the present study is eventually used. The particle size of a melt drop must be sufficiently small to be able to depict the miniscule melt fragments, but the coolant can be discretized using large particles. Therefore, a model to deal with particle interactions with different particle size under complex flow needs to be developed for particle method. Moreover, the phase transition of the melt and the coolant at the interface were not considered in the present hydrodynamic fragmentation mode, but they are important influencing factors in melt

12

G. Li et al. / Annals of Nuclear Energy 142 (2020) 107378

thermal fragmentation. The solid–liquid phase transition model can be implemented by changing particle types and properties when the particle liquid fraction reaches a certain value. However, for the coolant evaporation and steam condensation, the difficulty is the change in material volume because the coolant and the steam have considerable difference in density. As each particle has a constant volume and density, a potential approach to handle this problem is to vary the particle number as a means to maintain the mass conservation. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements This work was supported by the National Natural Science Foundation of China (Grant No. 11605129, 11975180) and the ‘‘Fundamental Research Funds for the Central Universities”. References Burger, M., Kim, D.S., Schwalbe, W., Unger, H., et al., 1984. Two-phase description of hydrodynamic fragmentation processes within thermal detonation waves. J. Heat Transfer. 106 (4), 728–734. Chen, X., Yuen, W.W., Theofanous, T.G., 1997. On the constitutive description of the micro interactions concept in steam explosions. Nucl. Eng. Des. 177, 303–319. Chen, X., Luo, R., Yuen, W.W., Theofanous, T.G., 1999. Experimental simulation of micro interactions in large scale explosions. Nucl. Eng. Des. 189, 163–178. Ciccarelli, G., Frost, D.L., 1994. Fragmentation mechanisms based on single drop steam explosion experiments using flash X-ray radiography. Nucl. Eng. Des. 146, 109–132. Duan, G., Chen, B., Koshizuka, S., Xiang, H., 2017. Stable multiphase moving particle semi-implicit method for incompressible interfacial flow. Comput. Methods Appl. Mech. Eng. 318, 636–666. Duan, R.Q., Koshizuka, S., Oka, Y., 2003. Numerical and theoretical investigation of effect of density ratio on the critical Weber number of droplet breakup. J. Nucl. Sci. Technol. 40 (7), 501–508. Escobar, S.C., Meignen, R., Rimbert, N., Gradeck, M., Secondary breakup in the context of fuel coolant interactions (FCI). Workshop Soci´et´e Franßcaise de Thermique, 2014, Paris, France. Escobar, S.C., Meignen, R., Picchi, S., A two-scale approach for modelling the corium melt fragmentation during fuel-coolant interaction. 16th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, 2015, Chicago, IL. Gelfand, B.E., 1996. Droplet breakup phenomena in flows with velocity lag. Prog. Energy Combust. Sci. 22, 201–265. Guildenbecher, D.R., Lopez-Rivera, C., Sojka, P.E., 2009. Secondary atomization. Exp. Fluids 46, 371–402. Han, J., Tryggvason, G., 2001. Secondary breakup of axisymmetric liquid drops. II. Impulsive acceleration. Phys. Fluids 13 (6), 1554–1565.

Hansson, R.C., Park, H.S., Dinh, T.N., 2009. Simultaneous high speed digital cinematographic and X-ray radiographic imaging of an intense multi-fluid interaction with rapid phase changes. Exp. Therm Fluid Sci. 33, 754–763. Haraldsson, H.O., Li, H.X., Yang, Z.L., Dinh, T.N., Sehgal, B.R., 2001. Effect of solidification on drop fragmentation in liquid-liquid media. Heat Mass Transf. 37, 417–426. Haraldsson, H.O., Breakup of jet and drops during premixing phase of fuel coolant interactions. Doctoral Thesis, Royal Institute of Technology, Stockholm, Sweden, 2000. Hsiang, L.P., Faeth, G.M., 1992. Near-limit drop deformation and secondary breakup. Int. J. Multiph. Flow 18 (5), 635–652. Kékesi, T., Amberg, G., Prahl Wittberg, P., 2014. Drop deformation and breakup. Int. J. Multiph. Flow 66, 1–10. Khayyer, A., Gotoh, H., 2011. Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230, 3093–3118. Kim, D.S., Burger, M., Frohlich, G., Unger, H., Experimental investigation of hydrodynamic fragmentation of gallium drops in water flows. International Meeting on Light-Water Reactor Severe Accident Evaluation, 1983, Cambridge, MA. Kondo, M., Koshizuka, S., Suzuki, K., Takimoto, M., Surface tension model using inter-particle force in particle method. 5th Joint ASME/JSME Fluids Engineering Conference, 2007, San Diego, California USA. Koshizuka, S., Oka, Y., 1996. Moving particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123, 421–434. Li, G., Liu, M., Duan, G., Chong, D., Yan, J., 2016. Numerical investigation of erosion and heat transfer characteristics of molten jet impinging onto solid plate with MPS-LES method. Int. J. Heat Mass Transf. 99, 44–52. Li, G., Liu, M., Wang, J., Chong, D., Yan, J., 2017. Crust behavior and erosion rate prediction of EPR sacrificial material impinged by core melt jet. Nucl. Eng. Des. 314, 44–55. Li, G., Liu, M., Wang, J., Chong, D., Yan, J., 2017. Numerical study of thermal erosion behavior of RPV lower head wall impinged by molten corium jet with particle method. Int. J. Heat Mass Transf. 104, 1060–1068. Li, Gen, Zhang, Jun, Yang, Qingzhong, Duan, Guangtao, Yan, Junjie, 2019. Numerical analysis of hydrodynamic fine fragmentation of corium melt drop during fuelcoolant interaction. Int. J. Heat Mass Transf. 137, 579–584. https://doi.org/ 10.1016/j.ijheatmasstransfer.2019.03.159. Nomura, K., Koshizuka, S., Oka, Y., Obata, H., 2001. Numerical analysis of droplet breakup behavior using particle method. J. Nucl. Sci. Technol. 38 (12), 1057– 1064. Park, H.S., Hansson, R.C., Sehgal, B.R., 2005. Fine fragmentation of molten droplet in highly subcooled water due to vapor explosion observed by X-ray radiography. Exp. Therm Fluid Sci. 29, 351–361. Popinet, S., 2003. Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries. J. Comput. Phys. 190 (2), 572–600. Popinet, S., 2009. An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228, 5838–5866. Popinet, S., Zaleski, S., 1999. A front-tracking algorithm for accurate representation of surface tension. Int. J. Numer. Meth. Fluids 30, 775–793. Strotos, G., Malgarinos, I., Nikolopoulos, N., Gavaises, M., 2016. Predicting droplet deformation and breakup for moderate Weber numbers. Int. J. Multiph. Flow 85, 96–109. Thakre, S., Ma, W., Li, L., 2013. A numerical analysis on hydrodynamic deformation of molten droplets in a water pool. Ann. Nucl. Energy 53, 228–237. Thorpe, S.A., 1968. A method of producing a shear flow in a stratified fluid. J. Fluid Mech. 32, 693–704. Uršic, M., Leskovar, M., Burger, M., Buck, M., 2014. Hydrodynamic fine fragmentation of partly solidified melt droplets during a vapour explosion. Int. J. Heat Mass Transf. 76, 90–98.