THEFIS Test Simulation to Validate a Freezing Model of ASTERIA-SFR Core Disruptive Accident Analysis Code

Abstract

The mechanical consequences of core disruptive accidents (CDAs) are a major safety concern in sodium-cooled fast reactors. Once core disruption occurs, liquefied core materials rapidly disperse vertically and radially. The dispersed materials penetrate the pin bundles and control rod guide tubes (CRGTs) before freezing at the edge of the penetration zone as heat is transferred to surrounding structures. Such freezing phenomena can suppress the negative reactivity feedback of fuel dispersion. The discharge of core materials can be impeded, resulting in a molten core pool formation when tight blockages occur inside CRGTs due to frozen material. Accordingly, freezing phenomena of core materials play a key role in governing the mechanical consequences of a CDA. To validate a freezing model implemented in our CDA analysis code, ASTERIA-SFR, a preliminary simulation of the THEFIS RUN#1 test, was performed. The calculation results show that freezing on the structural wall and crust formation were key phenomena affecting the penetration behavior, and the structural heat transfer is an important parameter. A remarkable reduction of the heat transfer coefficient was required to reproduce the penetration length observed in the experiment. This suggests that the momentum exchange and flow regime at the leading edge as well as heat transfer should be well modeled to predict the freezing phenomena in rapidly evolving CDAs.

1. Introduction

Sodium-cooled fast reactor (SFR) technology is one of the key candidates for future energy generation. Numerous studies have investigated SFRs because of their high core performance, such as the high-power density of the fuel, high burn-up over 100 GWd/t, the scope to reduce the amount of minor actinides and inherent safety features.

From the lessons learned from TEPCO’s Fukushima Daiichi power plants, countermeasures to prevent and/or mitigate the consequences of severe accidents are required for SFRs, as well as for light water reactors.

Despite the advantages of SFR technology, conventional SFR cores are not designed with an optimally reactive configuration in the nominal state. This means that in the event of a core disruptive accident (CDA), severe re-criticality governed by a prompt critical condition may be induced. This is accompanied by significant mechanical energy, which is sufficient to threaten the integrity of the reactor vessel.

To enhance the safety of SFRs, particular attention has been paid to the robustness of safety features in each level of the Defense-in-Depth (DiD) safety strategy [1]. With reference to the definition of DiD and plant state based on IAEA SSR2/1 [2] and INSAG-12 [3], the objective of the fourth level of DiD is to cope with accidents and control the plant condition, including the prevention of further accident progression and mitigation of the consequences. The plant design condition corresponding to the fourth DiD level, the so-called design extension condition (DEC) should be considered. An unprotected loss-of-flow (ULOF) event is considered a representative event for DEC evaluation because it can initiate CDA by the rapid evolution of the event to disrupt the entire core.

State-of-the-art SFR designs have sufficient flow coast-down and a well-designed system, and even without control rod scramming, SFR cores can maintain a fairly significant margin for fuel melting. Considerable efforts and resources have been devoted to developing countermeasures to avoid the consequences brought by severe re-criticality. Some examples include the use of a heterogeneous core with a negative void reactivity close to zero [4], and a FAIDUS-type fuel subassembly [5] for controlled material relocation before a significantly large molten fuel pool is formed. It is critical to confirm the effectiveness of the countermeasures against such severe accidents with reference to the regulation requirements. For example, the Nuclear Regulation Authority (NRA) of Japan has been developing a CDA analysis code, the ASTERIA-SFR code [6], that can be used to confirm the validity of such measures.

If molten core materials generated from failed fuel pins during an ULOF event are not largely discharged out of the active core region, a molten core pool may be created within a very short time. The molten core material from the pool penetrates the gaps in the remaining intact upper/lower pin bundles, gaps between inter-wrapper tubes and inside the control rod guide tubes. The penetration is enhanced by the fission product (FP) gas pressure and sodium vapor pressure developed by fuel–coolant interactions (FCIs). The extent of discharge of molten materials from the active core region is important in the following perspective. A previous study reported that discharging around 40% of the total core fuel inventory could help significantly reduce the re-criticality potential [7]. Conversely, freezing the leading edge of the penetrated core material could disturb the discharge of the following materials. When a molten core pool is created by residual materials in the core, radial pool motion, so-called “pool sloshing,” may be induced and could act as an initiator for severe re-criticality. The pool sloshing is driven by some triggers such as a fuel or steel vaporization by nuclear energy deposition in the pool [8]. Therefore, the penetration and freezing behavior of molten materials play key roles in re-criticality potential during a CDA.

An ULOF sequence is evaluated in multiple phases classified into initiating, transition, post-accident material relocation and post-accident heat removal phases, based on the dominant phenomenon in each phase [9]. Various calculation codes have been developed to evaluate the consequences of each phase. However, the event continuously proceeds through these phases and such an evaluation strategy may introduce additional uncertainties due to the discrepancy among the calculation models. For this problem, a combination of two different calculation codes is conventionally used to calculate the core behavior during ULOF, which is characterized by the radial motion of molten core materials that may induce re-criticality. One is a core behavior calculation code for the initiating phase, which calculates fuel pin and coolant behavior coupled with neutronics-based first-order perturbation theory, e.g., SAS4A [10]. The other is a calculation code to evaluate the core behavior during the transition phase, e.g., SIMMER-III [11] and SIMMER-IV [12]. The latter performs calculations based on space–time neutronics coupled with thermo-fluid dynamics in the transition phase using the calculated results from the initiating phase as the input. The combination of these two calculation codes is considered to have the best verification and validation results, reliability, and recognition among the current evaluation methods. However, it has limitations when integrated with different codes and models for each process. As a result of the integration, the evaluation results can include uncertainties caused by both the differences in the calculation models used in the codes and the selection of the time step when switching between the codes. Motivated by the need to avoid arbitrariness in the connection points and to provide an evaluation method with a unified model, new CDA analysis codes to conduct a unified calculation of initiating and transition phases are required.

In relation to the above, this study focuses on model validation for molten material penetration and freezing phenomena which can disturb the discharge of the molten core material from the active core region and may cause core pool formation during the initiating and transition phase of ULOF event.

Berthoud proposed an interfacial resistance model during molten fuel relocation on the structure surface [13]. Based on this interfacial resistance model, a fuel cap freezing model was developed to describe the freezing evolution by the formation of a super-cooled layer between the melt (crust) and the structural wall surface [14]. As the melt flows along the structural wall, it freezes and creates solid material caps. These caps then propagate a freezing zone comprising a gap and a super-cooled layer that merges with the crust component. In this study, the above-mentioned freezing model was implemented into CONCORD, a thermo-fluid dynamics calculation module of ASTERIA-SFR. This work was conducted as part of a study to validate ASTERIA-SFR, and the objective is to identify key phenomena and parameters in the freezing model. The validation was conducted by simulating THEFIS experiments [15] involving molten alumina injection into a quartz tube to observe the freezing behavior.

This paper describes CONCORD and the modeling of freezing phenomena in Section 2, model validation through THEFIS test analysis in Section 3, discussions in Section 4 and a summary in Section 5.

2. Overview of the Calculation Module, CONCORD

As shown in Figure 1, ASTERIA-SFR comprises three major modules: a thermo-fluid dynamics calculation module, CONCORD; a fuel behavior calculation module, FEMAXI-FBR; and a space–time neutronics calculation module, PARTISN/RKIN. Verification and validation studies have been performed over the past decade via test simulations and code-to-code comparisons [16,17]. The following section describes the physical models of CONCORD in more detail.

2.1. Modeling of CONCORD

2.1.1. Basic Design and Calculation Flow

CONCORD is a three-dimensional, multi-velocity field, multi-phase, multi-component and Eulerian fluid dynamics code. CONCORD comprises models such as flow regimes, interfacial area and phase change, e.g., vaporization/condensation and melting/freezing and equation-of-state of sodium coolant and core materials. The Fractional Step (FS) method is applied to the system with the phase change process taken into consideration. The FS method generally separates the convection term from other terms to solve conservation equations. As shown in Figure 2, the calculation flow of CONCORD can be roughly divided into three steps: A, B and C. The first step, A, constitutes the procedure for setting the initial state for the calculation cell of interest. By evaluating constitutive equations, the flow regime, interfacial area and binary contact area are each assessed in each calculation cell, whereupon momentum exchange, heat transfer between the components and phase change are calculated. At the end of Step A, basic conservative equations for mass, energy and momentum are solved. The next step, B, addresses iteration calculations for pressure. Mass and energy conservation equations are solved repeatedly until the residual error satisfies the criteria of the iteration calculation in CONCORD, then the iteration calculation is conducted to obtain consistent pressure between fluid-dynamically and thermally. In the final step, C, a final calculation is performed to calculate the advection flow under consistent pressure and update the mass, internal energy and momentum at the end of the time step.

2.1.2. Freezing Model

Freezing models are classified into two types: non-equilibrium phase change model and equilibrium phase change model. While the former is governed by limited heat transfer between two components via the interfacial area, the latter is a model based on the heat transfer between two components with bulk mass and internal energy. Mass exchange caused by phase change is calculated among vapor, liquid and solid phases. A set of heat transfer paths for vaporization/condensation (V/C) and for melting/freezing (M/F) are prepared for possible phase change of the components in an actual reactor and related experimental conditions.

A freezing model based on the non-equilibrium phase change model is prepared in CONCORD as follows.

Heat and mass transfer in a non-equilibrium state are calculated based on the balance of heat fluxes at the interface of the melt and structural wall. Net energy transfer at the melt–structure interface, $q_L$(W/m³), can be evaluated by the interfacial temperature $T_i$(K):

$$q_L=A_{bc}·\left[h_L\left(T_i-T_L\right)+h_s\left(T_i-T_s\right)\right],$$

(1)

where $A_{bc}$ (m²) is the binary contact area, h (W/m²/K) is the heat transfer coefficient, T (K) is temperature and subscripts L and s are the melt and structural components, respectively. Moreover, the mass transfer rate is determined by the net energy transfer $q_L$ divided by the specific enthalpy difference to solidus or liquidus points for the melt and structure. The interfacial temperature $T_i$ is evaluated based on the density ρ (kg/m³), specific heat $C_p$(J/kg/K) and thermal conductivity κ of each component by the following formulas:

$$T_i={\frac{h_L\sqrt{{\left(\rho C_p\kappa \right)}_L}{T}_L+h_S\sqrt{{\left(\rho C_p\kappa \right)}_S}{T}_S}{h_L\sqrt{{\left(\rho C_p\kappa \right)}_L}+h_S\sqrt{{\left(\rho C_p\kappa \right)}_S}}}_{}\mathrm{for}\mathrm{heat}\mathrm{transfer}\mathrm{between}\mathrm{the}\mathrm{melt}\mathrm{and}\mathrm{structural}\mathrm{wall}\mathrm{and}$$

(2)

$$T_i=\frac{h_L{T}_L+h_S{T}_S}{h_L+h_S}\mathrm{for}\mathrm{heat}\mathrm{transfer}\mathrm{between}\mathrm{the}\mathrm{melt}\mathrm{and}\mathrm{crust}.$$

(3)

Seban–Shimazaki’s formula [18], which is suitable for analyzing low-Prandtl-number fluids, is employed to calculate the heat transfer coefficient of melt, $h_L$ (W/m²/K), while the heat transfer coefficient of the structure, $h_S$ (W/m²/K), is calculated by thermal conductivity divided by half the thickness of the structural materials, as follows:

$$h_S=\frac{2k_S}{d}.$$

(4)

3. Model Validation through THEFIS Test Analysis

3.1. Overview of the THEFIS Experiment

The THEFIS test series comprises out-of-pile experiments performed by the Karlsruhe Institute of Technology (ex-KfK) in the 1980s to measure and observe the freezing behavior during molten alumina penetration in a thin quartz tube [15]. In this study, RUN#1 was selected for the test simulation as the most fundamental test, which assumes the injection of molten material unhindered by particle obstacles. The experimental setup comprises a pressurized volume occupied with molten alumina at 2300 K and a quartz tube 6 mm in diameter, 1 mm thick and 1.5 m high, filled with air at 293 K. The test was initiated from melt injection into the quartz tube from the bottom, driven by a pressure difference of approximately 0.1 MPa. Melt penetration behavior into the quartz tube was observed and the maximum penetration length was measured to be approximately 1350 mm in the RUN#1 test.

3.2. THEFIS Test Calculation

3.2.1. Calculation Conditions

The calculation system for the RUN#1 simulation is schematically depicted in Figure 3. In this study, a one-dimensional cylindrical calculation system is used. Since a single radial mesh is considered, material relocation is limited to the axial direction only. A total of 85 axial calculation cells are prepared: 10 for the pressurized melt zone and 75 for the gas zone. The tube wall is divided into three calculation nodes: inner, center and outer areas of the tube wall. Once the transient is initiated, melt is injected into the quartz tube at a pressure of 0.2 MPa and penetrates the tube. The thermal energy of the melt can be reduced due to heat transfer to a quartz tube, whereupon the melt is frozen onto the wall of the quartz tube. A crust is formed on the tube surface. As a result of heat transfer and momentum exchange among the melt, quartz tube and crust, the melt stops at a particular depth due to freezing.

3.2.2. Melt Freezing Behavior of the RUN#1 Test

(1)

Case 1

Case 1 is the calculation case using the original heat transfer model. The heat transfer coefficients for melt and structure are given by Seban–Shimazaki’s formula and Equation (4), respectively, as mentioned in Section 2.1.2. These coefficients can be calculated based on the geometry of the test equipment and the thermal properties of the materials. The energy transfer between melt and structure is calculated by using Equation (1). Furthermore, the interfacial temperature is calculated from Equations (2) and (3) based on the thermal properties of the materials. The calculation result using the original heat transfer model is shown in Figure 4. Note that VF10% and VF1% in the legends of Figure 4a refer to the uppermost location of the calculation cell, where molten alumina and crust occupy volume fractions of 10% and 1%. Figure 4b shows the material distribution, which expresses the volume fractions per calculation cell at 2.0 s after the onset of melt injection. The X axis indicates the volume fraction of each material component and the Y axis indicates the axial location. Because of the 1-D calculation, the left side of this figure indicates the centerline of the cylindrical calculation system. As shown in Figure 4a,b, the penetration did not terminate within 2 s in the calculation. The final penetration length was calculated to exceed 1800 mm and the structural quartz tube was melted in this case, although the test tube maintained its integrity during the test in reference [15]. Accordingly, it emerged that the interfacial temperature at the structure should be set to 2100 K, taking into account the measured temperature of another test with the non-melting wall structure in the THEFIS series [19].

(2)

Case 2

The result of Case 1 suggests that the calculation using these original heat transfer coefficients cannot reproduce the test results. This may be caused by insufficient contact (not perfect contact) between melt and structure, and it may induce thermal resistance of heat transfer between them. The purpose of Case 2 is to examine the effect of interfacial temperature. The calculation result assuming an interfacial temperature of 2100 K is shown in Figure 5. As shown in Figure 5a, the final penetration length was limited to 840 mm from the bottom of the test tube. This is approximately 510 mm shorter than the test result of RUN#1, namely 1350 mm. Figure 5b shows the material distribution during the initial 2.0 s after the melt injection onset. As shown in Figure 5b, molten alumina freezes on the surface of the structure due to contact with the structure and forms a crust. This crust then exerts thermal resistance during the heat transfer from the residual melt to the test tube, while crust growth impedes the penetration behavior due to the decrease in the cross-sectional area. In this case, the created crust occupied the majority of the hydraulic cross-section in the initial 2 s.

(3)

Case 3

The calculation result of Case 2 indicates imperfect contact between the melt and structural wall and thermal resistance at the interface. A calculation was conducted, assuming a 0.1 factor of structural heat transfer coefficient instead of considering the effect of thermal resistance. The 0.1 factor was determined based on a parametric calculation of the heat transfer coefficient for the structure described in Section 4.1 in more detail. The melt–structure interfacial temperature was deemed to be 2100 K. The transient of the penetration length and volume fractions of components are shown in Figure 6a,b, respectively. As shown in Figure 6a, the transient of the penetration length of Case 3 is enhanced compared with that of Case 2. The final penetration length was calculated to be 500 mm longer than Case 1, 1340 mm from the bottom of the test tube. Compared to the original case, this longer penetration length is a factor of the limited crust growth, due in turn to the thermal resistance of heat transfer between the melt and the test tube. The calculated final penetration length of 1340 mm, however, remains 10 mm shorter than the test result. The calculated penetration velocity into the tube is smaller than that of the test result, particularly for the initial 1 s after the melt injection onset.

4. Discussion

4.1. Effect of the Structural Heat Transfer Coefficient

When using an original freezing model with an interfacial temperature of 2100 K, as in Case 2, the penetration length is approximately 500 mm less than the experimental result. The underestimation indicated an excessive heat transfer from melt to structure. This resulted in a crust forming on the surface of the quartz tube, which was then occupied by frozen alumina, which impeded subsequent melt penetration.

Conversely, in Case 3, the calculated final penetration length was drastically improved due to limited evolution of crust formation compared with Case 2. This result is obtained assuming a 0.1 factor of structural heat transfer coefficient. This assumption was also made because it is consistent with the reference [13], namely that one order less of heat transfer from the melt leading edge to the structure would reproduce the penetration behavior. Figure 7 shows the relation between the factors of structural heat transfer coefficient, and the calculation results relate to the final penetration length. Clearly, the lower the coefficient, the longer the penetration. The case of a 0.1 factor of structural heat transfer coefficient (Case 3) comes closest to reproducing the depth of the test result. This indicates that the original heat transfer coefficient at the leading edge is one order of magnitude higher than the experimental condition. In the calculation, the volume fraction of the leading edge is sufficiently small, and the flow regime at the cell of the leading edge remains within the range of those dispersed (0–30% of the calculation cell) at all times. It is considered that there is a shortage of an appropriate flow regime for the reproduction of the experimental results.

On the other hand, the calculation could not reproduce the penetration velocity of the test even if any heat transfer at the leading edge of the melt decreased significantly. It is considered that less heat is transferred to the structural walls from the melt in the test case than in calculation case during penetration. The calculated small penetration velocity is considered to indicate a relatively low level of thermal resistance between crusts or the melt and structure in the calculation. The flow regime and related heat transfer model should be improved based on further investigation of the freezing mechanism.

4.2. Effect of the Initial Temperature

Reference [15] reported the test condition in which the melt was generated by ignition with aluminum–iron oxide thermite, and so the tests were performed using alumina of 2300 K ± 50 K with alumina superheated by 200 K. However, the equation of state may show variance between the test condition and the calculation, due to uncertainty over the melt composition in the initial state. Accordingly, the effect of initial temperature on penetration behavior was investigated by parametric calculations setting at 2524 K, considering a superheated temperature of 200 K and a 0.1 factor of structural heat transfer coefficient. Figure 8 shows the final penetration lengths of parametric cases for the initial temperature and the factor of the structural heat transfer coefficient. The X axis indicates the calculation conditions of each case, namely the initial temperature and the factors of the heat transfer coefficient of the structure (h_S). It emerged that adding 200 K to the initial temperature resulted in a final penetration length 200 mm longer. The penetration velocity of the calculation, however, remained smaller than that of the test result for all parametric calculation results in this study. This discrepancy between test and calculation results may be caused by momentum exchange between melt and gas. The calculated velocity of gas is greater than that of melt despite the fact that the calculated volume fraction of gas in the melt was small, 10⁻³~10⁻⁶. The assumed pressure difference of the boundary condition might not adequately contribute to the penetration of melt. Further investigation is needed, not only of the freezing mechanism but also of the momentum exchange and flow regime.

5. Summary

As part of the model validation study for ASTERIA-SFR, which is developed by NRA, a THEFIS test simulation was performed to identify the key phenomena and parameters of a freezing model of CONCORD, a thermo-fluid dynamics calculation module. It was identified that freezing on the structural wall and crust formation are key phenomena for penetration behavior, and the structural heat transfer coefficient is an important parameter in the freezing model.

Through the test simulation of RUN#1, the most fundamental test in the THEFIS test series, it emerged that freezing on the structural wall and crust formation governed the penetration length and the reduced heat transfer at the interface, as the thermal resistance proved decisive. In preparation for implementing the thermal resistance model into CONCORD, preliminary parametric calculations were performed for the structural heat transfer coefficient and initial melt temperature.

It emerged that the crust that formed on the quartz tube surface occupied most of the hydro-dynamic cross-section and impeded the following melt penetration. The calculated penetration length, which was 500 mm less than the experimental result, indicated that the heat transfer of the melt to the structure in the calculation was estimated to be larger than the experimental condition. Assuming a 0.1 factor of structural heat transfer coefficient, the penetration length improved, approaching the test result of 1350 mm. The freezing model can be improved by implementing the thermal resistance model taking account of outcomes from this study. In addition, as for another discrepancy that the penetration velocity of the calculation was smaller than that of the test result, further investigation and model improvement are planned, particularly on the momentum exchange between melt and gas. Model application to melt penetration behavior at lower energies shall be investigated through other THFIS tests with particle beds in the future.