Improvement of solidification model and analysis of 3D channel blockage with MPS method

Abstract

In a severe accident of a nuclear power reactor, coolant channel blockage by solidified molten core debris may significantly influence the core degradations that follow. The moving particle semi-implicit (MPS) method is one of the Lagrangian-based particle methods for analyzing incompressible flows. In the study described in this paper, a novel solidification model for analyzing melt flowing channel blockage with the MPS method has been developed, which is suitable to attain a sufficient numerical accuracy with a reasonable calculation cost. The prompt velocity diffusion by viscosity is prioritized over the prompt velocity correction by the pressure term (for assuring incompressibility) within each time step over the “mushy zone” (between the solidus and liquidus temperature) for accurate modeling of solidification before fixing the coordinates of the completely solidified particles. To sustain the numerical accuracy and stability, the corrective matrix and particle shifting techniques have been applied to correct the discretization errors from irregular particle arrangements and to recover the regular particle arrangements, respectively. To validate the newly developed algorithm, 2-D benchmark analyses are conducted for steady-state freezing of the water in a laminar flow between two parallel plates. Furthermore, 3-D channel blockage analyses of a boiling water reactor (BWR) fuel support piece have been performed. The results show that a partial channel blockage develops from the vicinity of the speed limiter, which does not fully develop into a complete channel blockage, but still diverts the incoming melt flow that follows to the orifice region.

Key words： boiling water reactor (BWR); severe accident; channel blockage; moving particle semi-implicit (MPS) method; solidification

1 Introduction

An accurate prediction of flow with computation fluid dynamics has been one of the major issues of nuclear engineering [1]. In a postulated severe accident of a nuclear power reactor, molten fuel and core structures, or mixtures thereof (corium) may flow down to the lower part of the core. During this process, corium may be cooled, which may solidify and narrow down or block the flow path. Such channel blockages may have a large influence on the core degradation process that follows as it may block the water/steam supply from the lower plenum of the reactor pressure vessel to the damaged core region. It may also develop to form a crust layer, which can hold up a molten corium pool above it [2]. Alternatively, the partial formation of a blockage may allow the continuous steam supply to the damaged core region while allowing relocations of the melt to the lower part of the core through the unblocked regions as observed in some experiments, such as XR-2 [3], and as expected to some extent in Unit-2 and Unit-3 of Fukushima Daiichi Nuclear Power Plant [4].

In the meantime, most severe accident analysis codes have limited capability in considering channel blockages with simple geometries and models. For example, MELCOR-2.2 assumes simple vertical channels, which are assumed to be blocked when refrozen materials completely fill the available sub-nodes based on energy balance calculations [5]. In reality, the channel geometries are complex as found in the core support structures and the complexity develops as core degradation proceeds. In general, modeling uncertainty for analyzing core degradation is large and is expected to be one of the key uncertainties, which is responsible for the diverged predictions regarding accident progressions of the Fukushima reactors [6] as also indicated by MAAP-MELCOR Crosswalk [7].

The moving particle semi-implicit (MPS) method is one of the particle methods for analyzing incompressible flows. It is a fully Lagrangian method, in which the fluid as well as the solid is discretized with calculation points (particles) without using calculation meshes [8]. The MPS method has been extensively developed and validated for analyzing the melt behavior relevant to severe accidents of nuclear power reactors, such as melt freezing in a penetration tube and a control rod guide tube [9–10], ablation of vessel wall with melt pool [11], melt spreading, and molten core-concrete interaction (MCCI) [12]. In general, the Lagrangian nature of the MPS method has advantages over the Eulerian based methods for analyzing melt behavior, when the analysis is focused on the evaluation of the melt/solid interface and interactions involved in dynamic flows, whose flow patterns are difficult to characterize prior to the analysis. These advantages arise from such natures of the MPS method as there is no need to explicitly evaluate the interface and free surfaces; there is no need to solve the advection term (convective term); and there are no calculation meshes.

In the meantime, the issues need to be addressed to further develop the MPS method for analyzing channel blockage in complex geometries include assuring incompressibility (keeping particle number density and particle arrangement uniform), which is important for assuring numerical accuracy of the Gradient, Divergence and Laplacian models; accurate modeling of solidification by solving the melt flow up to a sufficiently high level of viscosity before fixing the coordinates of fully solidified particles; and reducing the calculation cost, especially associated with pressure and viscosity calculations, which is the dominant calculation cost for the solidification analysis with MPS method.

In the preceding studies, solidification of melt was modeled by fixing the coordinates of solidified particles [9–10]. Such method is suitable for reducing the calculation cost, where a large part of the calculation domain is occupied by solidified particles, because both pressure and viscosity calculations can be skipped for solidified particles. However, the sudden fixing of solidifying particle coordinates tends to violate incompressibility (regular particle arrangement) and may lead to severe deterioration of numerical accuracy, such as melt particles penetrating through wall particles [10].

More specifically, Fig. 1 shows an example of the problem encountered in the previous research. It demonstrates the difficulty to attain qualitatively reasonable channel blockage analysis with the simple solidification model.

Fig.1 Unsuccessful calculation of fuel support piece channel blockage with melt freezing.

Full size|PPT slide

Alternatively, the solidification of the melt in a dynamic flow has been well modeled with increasing the viscosity of the melt for analyzing melt spreading [13]. Such method is superior to the former method with respect to assuring incompressibility and numerical accuracy. However, such method requires a heavy calculation cost associated with pressure and viscosity calculations of solidified particles. In the meantime, various techniques such as particle shifting have been developed to keep the regular arrangement of particles for improving the numerical accuracy of particle methods [14–15]. Corrective matrix has also been developed to correct discretizing errors associated with irregularity of particle arrangement [16].

With the above explained understanding, the study described in this paper aims to develop a novel solidification model (modifications to the original MPS algorithm) for analyzing melt flowing channel blockage using the MPS method. The novel method aims to attain sufficient numerical accuracy (at least, accurate enough to prevent unphysical behavior such as melt particles penetrating through the wall particles) with a reasonable calculation cost. In another word, in the preceding studies, the analysis of complex channel flow blockage was impossible, either because the solidification model was too simple, or the solidification model was too calculation cost demanding. This study aims to provide a new method with an appropriate simplification in modeling solidification (skipping some of the unnecessary calculations) so that practical channel blockage, such as the melt flow blockage problem in the boiling water reator (BWR) fuel support piece, can be analyzed with an affordable calculation cost. The basic approach is to divide solidification process to the following two steps to implement particle shifting and increasing the viscosity to a sufficiently high level, so that incompressibility is well guaranteed as the particles are being solidified to a partially solidified state through the “mushy zone” (between the solidus and liquidus temperature), and to fix the coordinates of the particles, when they are regarded as completely solidified and skipping unnecessary calculations for solidified particles, which would otherwise greatly increase the calculation cost (e.g., viscosity term and pressure term).

To validate the newly developed algorithm, 2-D benchmark analyses are conducted for the steady-state freezing of water in a laminar flow between two parallel plates [17]. Furthermore, the 3-D channel blockage analyses of BWR fuel support piece are also performed with a focus on how the flow path is partially blocked and how the incoming flow that follows is diverted and the flow pattern develops.

2 Calculation method

With regard to numerical accuracy, momentum conservation of the MPS method is often discussed with respect to other particle methods. For example, it is often discussed in comparison to the smoothed particle hydrodynamics (SPH) method, where the conservative formulations are usually adopted. As also discussed in the recent review for the grand challenges in SPH [18], it is well known that the conservative formulations usually perform well for stability but are not truly convergent on random particle distributions. Meanwhile, as argued by Jeong et al. [19], the conservative formation can cause additional numerical diffusion/dissipation especially for multiphase simulations. On the other hand, it is noted that the high-order MPS formulations using corrective matrix are more accurate but easily cause instability when the neighbor support is not complete (e.g., near free surface). In this paper, to achieve a good balance between numerical accuracy and stability, the non-conservative formulations with corrective matrix for internal particles to improve accuracy and the conservative formulations near free surface to enhance overall stability are employed. Therefore, such a method can enhance the accuracy for internal particles and guarantee the overall stability especially near the free surface. This idea has been adopted and validated in the MPS study in Ref. [16] and in a similar SPH study by Sun et al. [20]. Therefore, even though the employed formulations are not exactly conservative, satisfactory and reliable results can still be obtained. The other way is to apply the Hamiltonian MPS method, in which the momentum and mechanical energy of the system are conserved [8]. In the current paper, the accuracy of the pressure term is improved by the corrective matrix and the stability is guaranteed by introducing the particle stabilizing term [9]. More details regarding the employed discretization and stabilization models in this paper are explained in Section 2.1.

2.1 Governing equations and discretization models

The governing equations of the MPS method are conservation equations of mass, momentum, and energy, expressed as [8,21]

(1){\displaystyle \frac{\mathrm{d}\rho }{\mathrm{d}t}}=\mathrm{0},

(2){\displaystyle \frac{\mathrm{d}\mu }{\mathrm{d}t}}=-{\displaystyle \frac{\mathrm{1}}{\rho }}\nabla P+v{\nabla }^{\mathrm{2}}\mathit{u}+g+\mathit{f},

where ρ is the density, kg/m³; u is the velocity vector, m/s; P is the pressure, Pa; ν is the kinematic viscosity, m²/s; g is the gravity, m/s²; h is the enthalpy, kJ/kg; k is the thermal conductivity, W/(m∙K); T is the temperature; and Q is the heat source, J/m³, f is the force.

In the MPS method, the differential operators of the governing equations are discretized with particle interaction models such as gradient, divergence, and Laplacian models [8]. In this paper, the following models with corrective matrix are utilized to improve the numerical accuracy compared to the original MPS method [16].

(4){\langle \nabla P\rangle }_i={\displaystyle \frac{d}{n^{\mathrm{0}}}}{\displaystyle \sum _{j\ne i}\left({\displaystyle \frac{(p_j-P_i)}{\Vert r_j-r_i\Vert }}\left({\mathit{C}}_i{\displaystyle \frac{r_j-r_i}{\Vert r_j-r_i\Vert }}\right)w(\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert )\right\}},

(5){\langle \nabla ·\mathit{u}\rangle }_i={\displaystyle \frac{d}{n^{\mathrm{0}}}}{\displaystyle \sum _{j\ne i}\left\{{\displaystyle \frac{({\mathit{u}}_j-{\mathit{u}}_i)}{\Vert {\mathit{r}}_j-r_i\Vert }}\left({\mathit{C}}_i{\displaystyle \frac{r_j-r_i}{\Vert {\mathit{r}}_j-r_i\Vert }}\right)w(\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert )\right\}},

(6){\langle {\nabla }^{\mathrm{2}}\varphi \rangle }_i={\displaystyle \frac{d}{n^{\mathrm{0}}}}{\displaystyle \sum _{j\ne i}\left\{({\varphi }_j-{\varphi }_i)\left({\displaystyle \frac{\mathrm{2}}{\lambda }}-{\displaystyle \frac{{\mathit{L}}_i{\mathit{c}}_i(r_j-r_i)}{{\Vert r_j-r_i\Vert }^{\mathrm{2}}}}\right)w(\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert )\right\}},

where d is the dimension number, n⁰ is the initial particle number density, \varphiis a scalar quantity, and λ is the correction factor, which is defined as \lambda ={\displaystyle \frac{{\displaystyle {\sum }_{j\ne i}{{\displaystyle r}}_{ij}^{\mathrm{2}}{{\displaystyle w}}_{ij}}}{{\displaystyle {\sum }_{j\ne i}{{\displaystyle w}}_{ij}}}}.

The corrective matrix C_i and the row vector L_i are calculated by using

(7){\mathit{C}}_i^{-\mathrm{1}}={\displaystyle \frac{d}{n^{\mathrm{0}}}}\left(\begin{array}{ll}{\displaystyle \sum _{j\ne i}\frac{x_{ij}^{\mathrm{2}}}{r_{ij}^{\mathrm{2}}}}{w}_{ij}& {\displaystyle \sum _{j\ne i}\frac{x_{ij}{y}_{ij}}{r_{ij}^{\mathrm{2}}}}{w}_{ij}\\ {\displaystyle \sum _{j\ne i}\frac{y_{ij}{x}_{ij}}{r_{ij}^{\mathrm{2}}}}{w}_{ij}& {\displaystyle \sum _{j\ne i}\frac{y_{ij}^{\mathrm{2}}}{r_{ij}^{\mathrm{2}}}}{w}_{ij}\end{array}\right),

(8){\mathit{L}}_i=\left({\displaystyle \frac{\mathrm{2}d}{n^{\mathrm{0}}\lambda }}{\displaystyle \sum _{j\ne i}{x}_{ij}}{w}_{ij}\hspace{1em}{\displaystyle \frac{\mathrm{2}d}{n^{\mathrm{0}}\lambda }}{\displaystyle \sum _{j\ne i}{y}_{ij}}{w}_{ij}\right),

Where x_{ij}=x_j-x_i,y_{ij}=y_j-y_i,\hspace{1em}{r}_{ij}^{\mathrm{2}}=x_{ij}^{\mathrm{2}}+y_{ij}^{\mathrm{2}}, and w_{ij}=w(\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert ).

These models are applied using a weight function, as also utilized in the original MPS method. The weight function w(r) is represented as [8]

(9)w(r)=\{\begin{array}{l}{{\displaystyle \left(\mathrm{1}-{\displaystyle \frac{{{\displaystyle r}}_e}{r}}\right)}}^{\mathrm{2}}(\mathrm{0}\le r\le {{\displaystyle r}}_e),\\ \mathrm{0}({{\displaystyle r}}_e\le r),\end{array}

where ris the position vector [m], and r_e is the effective interaction radius [m].

The technique drives particles from a particle-dense area to a particle-diluted area by adjusting the particle coordinates. The displacement by particle shifting is determined as

(10)\delta {\mathit{r}}_i={\displaystyle \frac{\Delta r_i}{n^{\mathrm{0}}}}{\displaystyle \sum _{j\ne i}\frac{l_{\mathrm{0}}}{\Vert r_j-r_i\Vert }}\left({\mathit{C}}_i{\displaystyle \frac{{\mathit{r}}_i-{\mathit{r}}_j}{\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert }}\right)w(\Vert {\mathit{r}}_j-{\mathit{r}}_i\Vert ),

where \delta {{\displaystyle r}}_i is the displacement, m; l₀ is the particle size, m; and is the distance variable whose value depends on the value of r_i,min, which is the minimum distance between particles i and j. Specifically, it is represented as

(11)

where a is the threshold parameter to activate particle shifting. Some trial tests have shown that a = 0.9 is appropriate [16]. Therefore, when the minimum distance between particles i and j is smaller than 0.9 times the particle size, particle shifting is activated.

2.2 Phase change and viscosity models

In this paper, the temperature of each particle is calculated from its enthalpy, expressed as

(12)

where T_s is the solidus temperature, K; C_ps is the specific heat capacity of the solid, J/(kg∙K); C_pl is the specific heat capacity of the liquid, J/(kg∙K); h_s0 is the enthalpy of the material at solidus temperature, kJ/kg; and h_s1 is the enthalpy of the material at liquidus temperature, kJ/kg.

The solid fraction, which represents the volumetric fraction of the solid phase in one particle, is evaluated with the enthalpy at solidus and liquidus temperature as shown in Eq. (13). In this way, the solidification status of the particle between the solidus and liquidus temperatures is defined. When the solid fraction is 0, the state of particles is fully liquid. When the solid fraction is 1, the state of particles is fully solid. When the solid fraction is between 0 and 1, the state of particles is a mixture of solid and liquid, as is known in a mushy zone.

(13)

The viscosity of the mixture in the mushy zone has been experimentally measured for various melts and various correlations have been proposed [22–24]. Among those correlations, this study adopts the Ramacciotti model [24]. The correlation calculates the effective viscosity of the mixture as a function of solid fraction, expressed as

(14)\mu ={\mu }_l\exp [\mathrm{2.5}\cdot \gamma \cdot CR],

where μ is the dynamic viscosity, Pa∙s; μ_l is the dynamic viscosity of liquidus, Pa∙s; and CR is the parameter of Ramacciotti model (CR typically takes the value of 4.0–8.0 [24] and it is tentatively determined to be 7.0 in this paper).

2.3 Novel algorithm and solidification model

In Ref. [13], some modifications to the MPS algorithm were made to extend the applicable range of the method to the extremely highly viscous flow with an acceptable numerical accuracy so that a complete solidification can be practically represented by extremely viscous fluid. Basically, the order of velocity diffusion calculation by the viscosity and velocity correction for assuring incompressibility by the pressure gradient was switched in each time step. This modification was intended to avoid the so-called “numerical creep” of the particles, whose velocity should be practically zero due to extremely high viscosity. In another word, the velocity corrections from the pressure calculation is promptly diffused by viscosity in the same time-step (with a drawback of incompressibility correction delayed by one time-step). In addition, the particle shifting technique was implemented to improve incompressibility and numerical accuracy and stability. With these modifications, the developed algorithm could successfully model the solidification of melt with a high viscosity [13].

However, one problem that arises from the preceding study is the significant calculation cost increase as the calculation domain is occupied with practically “solidified particles,” which are still treated as “fluid type” particles in the MPS algorithm. In MPS method, different calculation algorithms can be applied for different “particle types.” Hence, in this paper, a new particle type, “immobilized type” is defined so that some unnecessary calculations can be skipped for the practically solidified particles to reduce the calculation cost.

Specifically, when the solid fraction of a fluid particle is below 0.5 (or any other appropriate fraction between 0 and 1.0), its viscosity through the mushy zone is increased with the solid fraction described in Eq. (14). When the solid fraction of a fluid particle reaches 0.5, its viscosity is further increased by a factor of 10 (or by any other appropriate factor that is large enough so that the movement of the particle is almost terminated). This strategy is to simulate the highly-viscous solidifying melt and will be kept until the solid fraction reaches 1.0. Once the solid fraction of a particle reaches 1.0, it is judged as being fully solidified and is converted from the “fluid type” to the “immobilized type” and its coordinate is fixed. It should be noted that the current immobilization method is not applicable to the free surface flow, which involves a significant solidification from the free surface, because solidified particles on the free surface should still be transported by the bulk melt flow. However, the current immobilization method (fixing the coordinates of fully solidified particles) seems reasonable when analyzing channel blockage by fluid solidification from the wall boundary. This algorithm is shown in Fig. 2(b) as Algorithm-1.

In addition, for assuring a more robust numerical stability, Algorithm-2, as shown in Fig. 2(c), has been developed. In this algorithm, when the velocity of the particle is below a threshold, some costly calculations such as viscosity calculations are skipped, while the particle type is still kept as fluid-type to assure a more numerical stability. In this manner, a considerable calculation cost can be saved while retaining the numerical stability and accuracy depending on the nature of calculations. More specifically, the omitted calculations for the particles at a solid fraction of 1.0 and a velocity below the threshold are described as follows: external force calculations including gravity, temporary movement of particles due to external forces, viscosity calculation and final movement of particles due to pressure and viscosity calculations. Such treatment seems reasonable for the current application, because freezing is expected to occur always from the wall boundary and any particles adjacent to the boundary with a solid fraction of 1.0 and a velocity below the threshold may be expected to anchor to the solidified wall.

Thus, the features and advantages of the developed method relative to the methods of the preceding studies may be summarized as follows:

Preceding study [13]: Accurate and generalized modeling of solidification is adopted, which can be applied to wide applications involving the free surface and multi-interface flow. However, it is accompanied with high calculation cost, because all solidified particles need to be treated as fluid particles with extremely large viscosity.

Algorithm-1: Calculation cost is substantially reduced, because almost all calculations are skipped for solidified particles. However, this method is applicable only to the cases where solidification builds-up from the wall boundary in simple geometry (e.g., 2-D parallel plate channel).

Algorithm-2: Calculation cost is moderate as some costly calculations (e.g., viscosity calculation) are skipped for solidified particles. This method is still only applicable to cases where solidification develops from the wall boundary. However, the applicability of the method is expected to extend to channels with complex geometries, because pressure calculations are still performed for the solidified particles to assure stable fluid – solid interface calculations.

Fig.2 MPS algorithms for modeling melt solidification.

Full size|PPT slide

3 Benchmark analysis of laminar water flow and freezing between two parallel plates

3.1 Experimental conditions

In this section, the experimental data of a laminar flow of freezing water between the cooled parallel plates [17] have been used to benchmark the developed MPS method. The experimental geometry, as depicted in Fig. 3, is composed of two sections. The acrylic plastic section in the upstream is for forming a laminar inflow of water under the adiabatic condition. The test section in the downstream consists of two parallel copper plates, which are kept at a constant temperature so that an ice layer can be developed from the inner plate surface with time until an equilibrium state is reached. In both sections, the side walls are made of acrylic plastic plates and are thermally insulated.

This experiment was performed in six cases at different temperatures and inflow rates. The equilibrium thicknesses of the developed ice layers were measured for each case. Table 1 summarizes the six test conditions with different inflow Reynolds number Re and dimensionless wall temperature θ_W, which represents the ratio of the wall temperature to the fluid temperature.

Fig.3 Experimental geometry.

Full size|PPT slide

Tab.1 Experimental cases [15]

Case	A	B	C	D	E	F
Re	700				1200	2300
θW	2.5	1.1	0.6	0.4	1.1

3.2 Analysis conditions

Figure 4 displays the analysis geometry of the MPS method in 2-D. The particle size is 0.4 mm in the reference case. In the simulation, the inflow boundary is simulated by giving an initial velocity to the inflow particles at the inlet. Furthermore, the test section length has been reduced from 400 mm to 300 mm to reduce the calculation cost. In the experiment, a constant flow velocity is given at the inlet. At the outlet boundary, the fluid particles are converted into ghost particles when they pass the given outlet line. In this situation, a free surface boundary condition is maintained for the remaining fluid particles, which is similar to the true outlet with a free surface in reality. It is well known that the gravity does not influence the steady solution of the Poiseuille flow in the channel. Therefore, the gravity is neglected in the simulation and thus the free surface outlet does not influence the upstream. The effectiveness of this kind of outlet boundary is also verified in Ref. [25]. As observed in the experiment [17], the last 100 mm did not show changes in the flow pattern and ice layer thickness. Therefore, it is reasonable to omit the calculation of the last 100 mm in experiment and consider that the outlet boundary does not have much effect to the upstream.

Fig.4 Analysis geometry in 2D.

Full size|PPT slide

The inflow rate {{\displaystyle \bar{u}}}_{mo} (mo represents average of channel width), wall temperature T_W, and inflow temperature T₀ are converted from the experimental conditions (Table 1) with the following relationships as explained in the experimental report and summarized in Table 2.

(15)Re={\displaystyle \frac{\mathrm{4}{\bar{u}}_{mo}H}{v}},

(16){\theta }_W={\displaystyle \frac{-T_W}{T_{\mathrm{0}}}},

where H is half of the flow channel width (10 mm), m and ν is the kinematic viscosity, m²/s.

Tab.2 Analysis cases

Case	A	B	C	D	E	F
Re	700				1200	2300
{{\displaystyle \overline{u}}}_{mo}/(\mathrm{m}\cdot {\mathrm{s}}^{-\mathrm{1}})	0.02625				0.045	0.08625
{{\displaystyle \theta }}_{\mathrm{W}}	2.5	1.1	0.6	0.4	1.1
{{\displaystyle T}}_{\mathrm{0}}/°\mathrm{C}	2.0	4.0	4.0	6.0	4.0
{{\displaystyle T}}_{\mathrm{W}}/°\mathrm{C}	–5.0	–4.4	–2.4	–2.4	–4.4

To investigate the applicability of the new algorithms, reasonable particle size should be determined, so that differences due to the fact that different algorithms can also be captured with considerations of the experimental results. The major difference between Algorithm-1 and Algorithm-2 is the consideration of the minimum velocity threshold. Therefore, velocity field should be reasonably considered in the simulations. Figure 5 exhibits the velocity scale in the direction perpendicular to the flow direction near the inlet for the three particles sizes of 0.2 mm, 0.4 mm, and 0.8 mm. As can be seen, the particle size of 0.8 mm cannot capture the velocity field well. The results indicate that a particle size of 0.4 mm or even smaller is necessary. Moreover, the ice thickness obtained by the experiment ranges from about 3.0–8.0 mm. With considerations of these factors, this study adopts a particle size of 0.4 mm from hereafter.

Fig.5 Velocity scale perpendicular to the flow direction near the inlet at different particle sizes.

Full size|PPT slide

3.3 Analysis results of different solidification algorithms

In this section, the results of the solidification algorithm of the preceding study in Ref. [13] are compared with that of the developed Algorithm-1 and those of the experiment. Figure 6 presents the final equilibrium state with the developed ice layer distribution as water is continuously flowing from the left to the right. As confirmed by the particles with zero pressure and initial viscosity values ( = 1.5 × 10⁻³ Pa∙s), the pressure and viscosity calculations of the immobilized type particles are skipped in Algorithm-1. In contrast, the ice particles have nonzero pressure values and high viscosity values to model solidification with the algorithm of the preceding study in Ref. [13].

Fig.6 Pressure and viscosity distributions in equilibrium state.

Full size|PPT slide

Consequently, the calculation cost of Algorithm-1 has been reduced from that of the preceding study. By using Intel Xenon processor E5-2650v2 8 CPU cores machine with 16 threads, the calculation time has been decreased from about 14 days to about 4 days.

Using the algorithm of the preceding study in Ref. [13], different cases of the developed Algorithm-1 and Algorithm-2 summarized in Table 2 are simulated. Figure 7 shows the final ice layer thickness distribution with different dimensionless wall temperatures (all with the same Re = 700). Similarly, Fig. 8 illustrates the results with different Re values (all with the same θ_W = 1.1). Considering the calculation resolution (with a particle size of 0.4 mm), the quantitative results shown in Fig. 7 for Re = 700 agree well with the measurements. On the other hand, the results shown in Fig. 8 indicate that the error between the algorithms becomes larger as Re increases. Qualitatively, it may be understood that Algorithm-1 tends to overestimate buildup of the ice layer relative to other algorithms, because only the solid fraction is considered as a condition to immobilize the solidified particles. In reality, even if the water is frozen by the ice layer surface, it may not always be immobilized on the ice layer surface as it may slip by the ice layer and flow out from the channel. The algorithm of the preceding study may be more accurate, although it depends on accurate modeling of interaction of the ice layer surface and the water flowing by the ice layer surface (i.e., viscosity interaction model of the water and the ice particles). Algorithm-2 may be regarded as in between the other two algorithms as the minimum velocity threshold can influence the ice layer buildup in addition to the water-ice viscosity interaction model. In summary, all algorithms inherit parameters, which need to be tuned for a given resolution and flow regime to quantitatively match the experimental results to a high degree of accuracy. Such investigation may be considered for future study. Moreover, all simulations are expected to underestimate the ice layer buildup as flow pattern becomes complex and involves significant flow velocity perpendicular to the ice layer surface, which cannot be captured with the current simulation resolution. When the channel wall surfaces are smooth, Re<2000 usually indicates that the flow is laminar. However, in the freezing flow, the freezing surface may not be very smooth, and the flow may turn into a turbulent flow even if Re<2000. It is also noted by the original experimental report that the ice layer development is different from the theoretical prediction, which is assumed to be a perfect laminar flow [17].

Fig.7 Ice thickness for cases at different dimensionless wall temperatures (Re = 700).

Full size|PPT slide

Fig.8 Ice thickness for cases at different Re values (θ_W = 1.1).

Full size|PPT slide

4 Analysis of channel blockage in BWR fuel support piece

4.1 Analysis configurations and conditions

MPS simulations of metallic melt flowing through a BWR fuel support piece are conducted utilizing the developed Algorithm-2, which is more robust for modeling solidification involving complex free-surfaces and liquid-solid interfaces than Algorithm-1. The 1/4 symmetric calculation geometry covers one of the four orifice flow paths and 1/4 of the central flow channel and the cruciform control blade as shown in Fig. 9. The 1/4 symmetry is represented by the dummy wall particles, which acts as the no-slip reflection boundary. The particle size is tentatively determined as 2.0 mm and the inflow particles are arranged as 20 mm × 20 mm at the center of the cruciform channel. The fuel support piece and the velocity limiter are assumed to be made of stainless steel (SS).

Different analysis cases have been determined by referring to the XR2-1 experiment in Ref. [3]. However, the current simulations do not directly correspond to XR2-1, because the current simulation covers only the vicinity of the fuel support piece, while the XR2-1 is an integral experiment, which covers much a larger section of a reactor core. It should be noted that although there are past integral tests on the melting of BWR fuel assembly such as XR2 [3], the experimental data especially designed for the BWR fuel support piece is scarce and the data that can be directly used for validation of numerical simulation are not available to the best knowledge of the authors. This study shows that it is possible to reproduce the flow blockage as a mechanism. The comparison to an experiment devoted to the melting or damage of the BWR support piece can be reserved for the future study.

Tables 3 and 4 summarize the physical properties of the melt and the different simulation cases, respectively. All structures are initially assumed to be at 600 K as indicated by the XR2-1 experimental report. For the cases with Zr melt, sensitivity cases are prepared with a higher flow rate (Zr-FlowRate_H) and a higher superheat (Zr-Superheat_H). Assuming the representative length to be half of the channel width (which is the assumption in Section 3), the Re for the SS/B4C case, the Zr case, and the Zr-Superheat_H case are all about 600 and that of the Zr-FlowRate_H is about 1200, which seems to be reasonably covered by the conditions validated in Section 3.

Tab.3 Physical properties [10,26]

Property	SS/B4C	Zr	SS
Density/(kg·m–3)	6646	6520	7930
Specific heat/(J·(kg·K)–1)	452	377	840
Latent heat/(kJ·kg–1)	289	230	268
Melting temperature/K	1420	2100	1700
Thermal conductivity/(W·mK–1)	30.8	36.0	21.0

Tab.4 Analysis cases

	Different simulation cases
	SS/B4C	Zr	Zr-FlowRate_H	Zr-Superheat_H
Melt	SS/B4C	Zr	Zr	Zr
Inflow mass/kg	18	29	29	29
Inflow rate/(kg·s–1)	0.12	0.12	0.24	0.12
Initial melt temperature and superheat/K	1430/10	2110/10	2110/10	2310/210

Fig.9 Calculation geometry and boundary particles.

Full size|PPT slide

4.2 Results and discussions on partial channel blockage and melt flow distributions

Figures 10 and 11 show the development of the internal melt/solidified debris distributions with colors indicating the solid fraction at some representative moments until total inflow masses have reached 18 kg for both cases. For both cases, the melt solid fraction increases downstream as it loses its heat to the fuel support piece structure, and channel blockage develops from vicinity of the velocity limiter structure upstream. For both cases, the control blade channel is never fully blocked, but the inflow melt is diverted to the adjacent orifice channel as the partial blockage develops toward the melt injection point. Compared with SS/B₄C, Zr has a higher initial temperature, a smaller specific heat, a larger conductivity, and a smaller latent heat. Hence, the blockage development is more significant with Zr and there is no outflow of Zr through the fuel support piece, whereas some SS/B₄C melt has flowed through without solidifying in the fuel support piece.

Fig.10 MPS simulation results of SS/B₄C flowing down BWR fuel support piece.

Full size|PPT slide

Fig.11 MPS simulation results of Zr flowing down BWR fuel support piece.

Full size|PPT slide

The calculated outflow masses through the fuel support piece for the three cases, SS/B₄C, Zr-FlowRate_H, and Zr-Superheat_H, are shown in Fig. 12. The total outflow masses in the SS/ B₄C, Zr-FlowRate_H, and Zr-Superheat_H are 10.0 kg, 4.8 kg, and 3.9 kg, with an outflow/inflow mass percentage (volume ratio) of 55.6 %, 16.6 %, and 13.4 %, indicating that more significant channel blockage by melt solidification has been expected in the Zr-FlowRate_H, and Zr-Superheat_H cases. For the SS/B4C case, the blockage in the control blade channel did not develop enough to divert the inflow melt to the orifice channel and about half of the total inflow mass flowed through. The volume fraction occupied by the solidified melt relative to the original channel space (control blade channel+ orifice channel) was merely about 3%.

For the two Zr cases, the outflow through the control blade channel was temporarily stopped as the blockage developed in the control blade channel until the diverted melt started to flow out of the orifice channel. As can be seen from the outflow masses before the “plateaus” of Fig. 12, in both the Zr-FlowRate_H and Zr-Superheat_H cases, the outflow masses of the melts from the control blade were about 2.8 kg. It could be understood that the outflow masses from the control blade channel is the melt mass that can flow out of the channel, before the walls around the melt flow path are fully covered by the solidified melt, which does not depend much on the flow rate or superheat as any Zr in direct contact with the wall is effectively solidified almost instantaneously. Note that the injected melt is diverted to the orifice channel as soon as the blockage develops from the downstream to the upstream injected point. Then, the diverted melt exited from the orifice outlet. The total melt masses, which exited from the orifice channel for the Zr-FlowRate_H and Zr-Superheat_H cases were 2.0 kg and 1.1 kg, respectively. The outflow masses from the orifice channel is different among different cases, because the melt can still flow out of the orifice channel until the channel is fully blocked by the solidified melt.

Eventually, the orifice channels were completely blocked for the two Zr cases, while the control blade channel remained partially blocked (but no melt went through the control blade channel after the melt had been diverted to the orifice channels). The volume fractions occupied by the solidified melt relative to the original channel space for the Zr, Zr-FlowRate_H and Zr-Superheat_H cases were 18%, 20%, and 21%, respectively.

Thus, the current simulation results indicate that if the channel box has been lost and melt inflow mainly takes place near the center of the control blade channel, the channel may be only partially blocked and the inflow melt can be diverted to the adjacent orifice channel, which may eventually be fully blocked, while the control blade channel remains partially unblocked with relatively large space still remaining for any steam to flow through.

Fig.12 Comparison of outflow mass evaluation.

Full size|PPT slide

5 Conclusions and future issues

A novel solidification model for analyzing melt flowing channel blockage utilizing the MPS method has been developed, which can attain a sufficient numerical accuracy with a reasonable calculation cost. The solidification process is modeled by increasing the melt viscosity to sufficiently high level before solidification, by converting the particle type from fluid-type to immobilized-type at solidification, and by fixing coordinates and skipping some of the calculations unnecessary for the immobilized-type particles. The developed method has a complex melt flow/channel blockage behavior such as diverting inflow melt to adjacent flow channel when some of the channel area is partially blocked.

For future studies, appropriate modeling for melt-structure interface heat transfer may be conducted for more quantitative analyses of melt flowing through structures and flow channels during severe accident conditions of light water reactors. Considerations of entectic reactions may also be necessary. Such analyses may also help identifying possible core degradation progressions and the current damaged core status of Fukushima Daiichi Nuclear Power Plant.

Acknowledgments

This work was supported by the Nuclear Energy Science & Technology and Human Resource Development Project (through concentrating wisdom), the Advanced Theoretical and Experimental Physics, Waseda University, and the TCU priority promotion research supported by Tokyo City University.

References

1 Yu Y, Yang Y. URANS simulation of the turbulent flow in tight lattice bundle. Frontiers in Energy, 2011, 5(4): 404–411

2 Tolman E L, Kuan P, Broughton J M. TMI-2 accident scenario update. Nuclear Engineering and Design, 1988, 108(1–2): 45–54

3 Gauntt R O, Humphries L L. Final results of the XR2–1 BWR metallic melt relocation experiment. Office of Scientific and Technical Information (OSTI), 1997

4 Tokyo Electric Power Company. Evaluation of the situation of cores and containment vessels of Fukushima Daiichi Nuclear Power Station Units-1 to 3 and examination into unsolved issues in the accident progression progress. Tokyo Electric Power Company, Report No. 4, 2015

5 Sandia National Laboratories. MELCOR Computer Code Manuals Vol. 2: Reference Manual Version 2.2.9541 2017. SAND2017–0876 O, 2017

6 Organisation for Economic Co-operation and Development, Nuclear Energy Agency. Benchmark Study of the Accident at the Fukushima Daiichi Nuclear Power Plant (BSAF Project) Phase 1 Summary Report. NEA/CSNI/R18, 2015

7 Luxat D. Modular Accident Analysis Program (MAAP)-MELCOR Crosswalk Phase 1 Study. Technical Update, 2014

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

9 Chen R, Oka Y. Numerical analysis of freezing controlled penetration behavior of the molten core debris in an instrument tube with MPS. Annals of Nuclear Energy, 2014, 71: 322–332

10 Goto Y, Yamaji A. Analysis of eutectic and metallic melt flow and blockage in BWR control rod guide tube by MPS method. In: Proceedings of TopSafe 2017, Vienna, Austria, 2017

11 Takahashi N, Duan G, Furuya M, Analysis of hemispherical vessel ablation failure involving natural convection by MPS method with corrective matrix. International Journal of Advanced Nuclear Reactor Design and Technology, 2019, 1: 19–29

12 Yamaji A, Li X. Development of MPS method for analyzing melt spreading behavior and MCCI in severe accidents. Journal of Physics: Conference Series, 2016, 739(1): 012002(in Chinese)

13 Duan G, Yamaji A, Koshizuka S. A novel multiphase MPS algorithm for modeling crust formation by highly viscous fluid for simulating corium spreading. Nuclear Engineering and Design, 2019, 343: 218–231

14 Xu R, Stansby P, Laurence D. Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach. Journal of Computational Physics, 2009, 228(18): 6703–6725

15 Khayyer A, Gotoh H, Shimizu Y. Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context. Journal of Computational Physics, 2017, 332: 236–256

16 Duan G, Koshizuka S, Yamaji A, An accurate and stable multiphase moving particle semi-implicit method based on a corrective matrix for all particle interaction models. International Journal for Numerical Methods in Engineering, 2018, 115(10): 1287–1314

17 Kikuchi Y, Shigemasa Y, Oe A, Steady-state freezing of liquids in laminar flow between two parallel plates. Journal of Nuclear Science and Technology, 1986, 23(11): 979–991

18 Vacondio R, Altomare C, De Leffe M, Grand challenges for smoothed particle hydrodynamics numerical schemes. Computational Particle Mechanics, 2021, 8: 575–588

19 Jeong S M, Nam J W, Hwang S C, Numerical prediction of oil amount leaked from a damaged tank using two-dimensional moving particle simulation method. Ocean Engineering, 2013, 69: 70–78

20 Sun P N, Colagrossi A, Le Touzé D, Extension of the δ-Plus-SPH model for simulating Vortex-Induced-Vibration problems. Journal of Fluids and Structures, 2019, 90: 19–42

21 Koshizuka S, Shibata K. MPS-SW-MAIN-Ver.2.0. P 8827–1. 2006

22 Spindler B, Veteau J M. Simulation of spreading with solidification: assessment synthesis of THEMA code. Rapport CEA-R-6053, 2004

23 Alsmeyer H, Cron T, Foit J J, Test report of the Melt Spreasding Tests ECOKATS-V1 and ECOKATS-1. 2004

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

25 Duan G, Chen B. Large Eddy Simulation by particle method coupled with Sub-Particle-Scale model and application to mixing layer flow. Applied Mathematical Modelling, 2015, 39(10–11): 3135–3149

26 Kirillov P L, Terentieva M I, Deniskina N B. Thermophysical Properties of Materials for Nuclear Engineering. 2nd ed. Moscow: IzdAT Publishing House, 2007