<figure><img src="image_1.png" /></figure>

HAL
open science

Evaluation of Sodium Boiling Models Using KNS-37
Loss of Flow Experiments

Sara Perez-Martin, Marine Anderhuber, Laurent Laborde, Nathalie Girault,
Calogera Lombardo, Luca Ammirabile, Konstantin Mikityuk, Stephane
Mimouni, Christophe Péniguel, Werner Pfrang

► To cite this version:

Sara Perez-Martin, Marine Anderhuber, Laurent Laborde, Nathalie Girault, Calogera Lombardo, et al.. Evaluation of Sodium Boiling Models Using KNS-37 Loss of Flow Experiments. Journal of Nuclear Engineering and Radiation Science, 2021, 8 (1), 10.1115/1.4050769. irsn-04112113

HAL Id: irsn-04112113

https://irsn.hal.science/irsn-04112113

Submitted on 8 Sep 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A review of models for the sodium boiling phenomena in sodium-cooled fast reactor subassemblies

**Haileyesus Tsige-Tamirat**

European Commission / Joint Research Centre Petten
Westerduinweg 3. 1755LE Petten
[Haileyesus.TSIGE-TAMIRAT@ec.europa.eu](mailto:Haileyesus.TSIGE-TAMIRAT@ec.europa.eu)

**Sara Perez-Martin**

Karlsruhe Institute of Technology
Institute of Neutron Physics and Reactor Technology
Hermann-von-Helmholtz-Platz 1
76344 Eggenstein-Leopoldshafen, Germany
[sara.perez@kit.edu](mailto:sara.perez@kit.edu)

**Werner Pfrang**

Karlsruhe Institute of Technology
Institute of Neutron Physics and Reactor Technology
Hermann-von-Helmholtz-Platz 1
76344 Eggenstein-Leopoldshafen, Germany
[werner.pfrang@kit.edu](mailto:werner.pfrang@kit.edu)

**Marine Anderhuber**

Atomic and Alternative Energies Commission,
Nuclear Energy Division, CEA Saclay
F-91191 Gif-sur-Yvette, France
[Marine.anderhuber@cea.fr](mailto:Marine.anderhuber@cea.fr)

**Antoine Gerschenfeld**

Atomic and Alternative Energies Commission,
Nuclear Energy Division, CEA Saclay
F-91191 Gif-sur-Yvette, France
[antoine.gerschenfeld@cea.fr](mailto:antoine.gerschenfeld@cea.fr)

**Laurent Laborde**

Institut de Radioprotection et de Sûreté Nucléaire (IRSN),
13115 Saint Paul Lez Durance Cedex, France
[laurent.laborde@irsn.fr](mailto:laurent.laborde@irsn.fr)

**Konstantin MIKITYUK**

Paul Scherrer Institut (PSI)
5232 Villigen, Switzerland
[konstantin.mikityuk@psi.ch](mailto:konstantin.mikityuk@psi.ch)

**Christophe Peniguel**

EDF-R&D, EDF Lab Chatou
6 Quai Watier, 78400 Chatou, France
[christophe peniguel@edf.fr](mailto:christophe peniguel@edf.fr)

**Stephane Mimouni**

EDF-R&D, EDF Lab Chatou
6 Quai Watier, 78400 Chatou, France
[stephane.mimouni@edf.fr](mailto:stephane mimouni@edf.fr)
ABSTRACT

The Euratom Horizon-2020 project ESFR-SMART aims to enhance the safety and performance of the European Sodium-cooled Fast Reactor (ESFR) considering safety objectives envisaged for Generation-IV reactors and the update of European and international safety frameworks, taking into account the Fukushima accident. Further, the project aims to support the development and validation of the computational tools for the situations to be considered in each defense-in-depth level in order to support the safety assessments using data produced in the project as well as selected legacy data. Within this activity, the focus is on the assessment and when needed further development of computer codes for the analysis of the sodium thermal-hydraulics phenomena in SFR subassemblies under operational and accidental conditions including sodium boiling and transitional convection from forced to natural/mixed convection. In support of this activity, a review is performed for the sodium boiling models used in the codes participating in benchmark activities within the project. The objective of the present paper is to summarize the result of the review, which encompasses both the phenomenological and mathematical models implemented in the codes. In particular, the review addresses the physical bases of sodium boiling models, the analytical models in connection with the numerical implementations in the codes and the geometry representations ranging from one-dimensional single channel and finite volume to full three-dimensional CFD and porous media representations. Finally the needs for further developments are discussed.
# 1 INTRODUCTION

The concept of the European Sodium-cooled Fast Reactor (ESFR) has been proposed within the FP7 Euratom project CP-ESFR [1] and is currently being improved in the Euratom Horizon-2020 project ESFR-SMART [2]. Within the ESFR-SMART project, new safety provisions are proposed for the ESFR concept considering safety objectives envisaged for Generation-IV reactors and the update of European and international safety frameworks [3] taking into account the Fukushima accident. Furthermore, the project aims to support the development and validation of the computational tools for the situations to be considered in each defense-in-depth level in order to support the safety assessments using data produced in the project as well as selected legacy data. Within this activity, the focus is on the assessment and when needed further development of computer codes for the analysis of the sodium thermal-hydraulics phenomena in SFR subassemblies under operational and accidental conditions including sodium boiling and transitional convection from forced to natural/mixed convection. In support of this activity, the sodium boiling models used in the codes participating in benchmark exercises within the project are reviewed. The ability of the computer codes participating in the current study to predict sodium boiling in separate effect and integral experiments is of fundamental importance in order to confirm the adequacy of the codes for use in safety analysis of full size reactor systems.

The objective of the present paper is to summarize the result of the review which encompasses both the phenomenological and the mathematical models implemented in the codes. In particular, the review addresses (i) the physical bases of sodium boiling models; (ii) the geometry representations ranging from one-dimensional single channel and finite volume to full three-dimensional CFD and porous media; and (iii) the analytical models in connection with the numerical implementations in the codes. Moreover, the modelling approaches to the various phenomena are discussed and the capability of the codes to model specific phenomenon such as boiling onset, nucleation, superheat, bubble growth, two-phase flow regime, liquid vapor interface, bubble collapse, re-entry and dry-
out is indicated. In the final part of the paper, the needs for further developments are discussed.

## 2 Physical model of Sodium boiling

### 2.1 Sodium boiling in SFRs

The sodium boiling is an important phenomenon in the safety assessment of SFRs. As consequence of the sodium boiling, the voiding of the reactor core could result in significant reactivity insertion which could in turn lead to an increase in the reactor power. If not mitigated, the boiling process could deteriorate to film boiling resulting in dry-out and overheating of the fuel pins ending in pin failure. However, in normal reactor operation, the sodium coolant is highly subcooled with a bulk outlet temperature 300 – 400°C below the saturation temperature. This provides a large safety margin to coolant boiling.

Despite the considerable safety margin, the possibility of sodium boiling is investigated in the safety assessment of SFRs. In particular, it plays a key role in the assessment of postulated hypothetical accidents such as subassembly blockage, unprotected loss of flow (ULOF), unprotected transient overpower (UTOP) and unprotected loss of ultimate heat sink (ULOHS). Although the probability of occurrence of these accidents is extremely low, their analysis is main part of the SFR safety assessment.

### 2.2 Sodium boiling in pin bundles

In case of subassembly blockage and in the initiation phase of the hypothetical accidents such as ULOF, UTOP, and ULOHS, transient sodium boiling within intact subassembly coolant channel could occur and if not mitigated could result in channel voiding. The key concerns in the safety assessment are [4]:

* the reactivity variations rate,
* the deterioration of the heat transfer conditions,

* the pressure surges due to boiling and

* subsequent condensation of the vapor in the subcooled region of the subassembly.

The focus of the present review is the phenomena related to transient boiling within intact subassembly coolant channel. Although, the basic phenomena involved in transient sodium boiling have been thoroughly investigated in the past, it remains a significant challenge to develop a mechanistic representation of the sodium boiling phenomena in an SFR subassembly for accident analysis computer codes.

The process of sodium boiling involves several complex sub-phenomena such as nucleation, boiling inception, superheat, vapor growth, vapor collapse in subcooled region, liquid vapor interface, liquid film characteristics, two-phase flow regime, etc. [4]. Further complication arises due to coolant channel geometry, flow blockage, flow reversal, rewetting, dry-out of the fuel pin, etc. As consequence of sodium boiling, transient two-phase flow regime with severe increase of the friction pressure drop with potential for flow excursion could establish. The two-phase flow regime is characterized by the existence of distinct single phase regions which are separated by moving interfaces [5]. In connection with the complex subassembly geometry, the accurate modelling of the transient two-phase flow is among the main challenges of the accident analysis.

In SFR subassembly (rod bundles), following sodium boiling processes may be distinguished [6]:

* Hot spot boiling: is localized boiling in a subchannel at nucleation points;

* Local boiling: is distributed boiling spreading over several subchannels; and

* Generalized or saturated boiling is propagated boiling which extends over a whole subassembly and beyond.

Hot spot boiling is triggered by local temperature rise resulting from flow disturbance initiated by the distortion of the flow channel geometry. It is mainly prevalent
under nominal condition where the heat fluxes are large. Local boiling may appear as a result of radial temperature gradients in a subassembly following the occurrence of a transient. If not controlled, local boiling may develop to saturated boiling.

Local boiling and hot spot boiling are phenomena producing small local bubbles which immediately collapse in the cooler surrounding. The impact on cross section averaged coolant density and axial coolant flow is very small and in consequence the associated coolant reactivity feedback is negligibly small.

Modeling of generalized boiling of sodium under transient boiling conditions requires consideration of following key phenomena [4]:

* nucleation, boiling inception, superheat,

* vapor growth and vapor collapse in subcooled region,

* phase interface, interfacial heat transfer and momentum exchange,

* liquid film characteristics, dry-out,

* two-phase flow regime and

* reentry and rewetting.

Many codes represent the subassembly cross section by an average pin with its associated coolant area and structure material (e.g. SAS-SFR, ASTEC-Na, CATHARE-3). This approach is assumed sufficient to calculate generalized boiling. To get more insight into local incoherencies and to detect also hot spot and local boiling it is necessary to model the complex coolant channel geometry in detail. With the concurrence of the transient two-phase flow simulation this represents a huge challenge. However, considerable progress has been made in the past to enhance the understanding of the sodium boiling in analytical model developments supported by dedicated separate effect and integral experiments.
## 2.3 Sodium boiling modeling

In what follows, the sodium boiling models in the computer codes participating in benchmark activities of the ESFR-SMART project are reviewed encompassing both the implemented phenomenological and mathematical models. The involved computer codes and the participating institutions are:

* SAS-SFR, Karlsruhe Institute of Technology (KIT);

* NEPTUNE_CFD, Électricité de France (EDF);

* ASTEC-Na, Institut de radioprotection et de sûreté nucléaire (IRSN);

* CATHARE-3, Commissariat à l'énergie atomique et aux énergies alternative (CEA); and

* TRACE, Paul Scherrer Institute (PSI);

# 3 Physical models of sodium boiling

For the modeling of the transient two-phase flow resulting from sodium boiling, the reviewed computer codes use essentially two approaches: multi-bubble and two fluid models. In what follows, these physical bases for the models in the codes are discussed. The limitation of the models is highlighted in section 5.

## 3.1 Overview of physical models

The sodium boiling in the SAS-SFR code [7] is modeled using a one-dimensional multi-bubble slug ejection model [8] which allows a finite number of bubbles separated by liquid slugs in a channel at any time. A representation of the model is depicted in **Figure 1**. Assuming a certain degree of superheat, sodium bubbles are generated when the
sodium temperature exceeds the saturation temperature by a user-supplied superheat.
The new bubble both expands and moves away from the heated region in the direction
of sodium flow, thereby allowing collapsing and forming of additional bubbles in the
heated region. The coolant channel is thus modelled as a series of liquid slugs and vapor
bubbles moving along the coolant channel. The model assumes two types of bubbles,
small ones where the pressure is uniform and larger ones where the pressure varies
axially. Both types of bubbles extend radially in the whole cross section of the coolant
channel leaving only a thin liquid film with a user defined initial thickness on the surfaces
(cladding, wires and hexcan surfaces). The axial extent of the bubbles and their dynamic
behaviour is determined mainly by bubble-slug interfaces and the liquid slug motion.
Therefore, the voiding calculation couples vapor pressure calculations for the bubbles
with momentum equations for the liquid slugs. Vapor flow of sodium between the slugs
is modelled as annular flow with a static liquid film of variable thickness wetting cladding
and structural surfaces. For small bubbles, the uniform vapor pressure throughout the
bubble corresponds to the actual saturation temperature where the flow rate is governed
by the momentum balance for the liquid slugs motion. Mass is computed through the
mass balance for the vapor phase. For axially extended bubbles, the vapor production on
hot clad surfaces and condensation on cold clad and structure surfaces can result in large
vapor velocities and appreciable axial pressure gradients. The axial distributions of the
vapor pressure, temperature and mass flow rate are calculated using the continuity
equation, the momentum equation and the mass conservation equation assuming local
saturation conditions. The vapor production and condensation rates relative to clad
surface and structure temperatures are determined by heat flow through liquid films on
clad and structure i.e. the liquid-vapor interfaces. The dominating heat transfer
coefficient for vapor production is determined by the ratio of the thermal conductivity
divided by the actual film thickness. The condensation heat transfer coefficient considers
the inert gas content in the vapor phase. The initial film thicknesses are specified by the
user based on coolant cross section and wetted surfaces. The recommended value is that
corresponding to a total liquid volume fraction of about 10% in the voided region. The
film thickness changes due to vaporization, condensation and redistribution of the liquid inventories of slugs in case of too small slugs leading to slug disruption. Film dry-out on clad and structures is calculated on basis of criteria related to minimum film thickness or flooding criteria. Heat transfer of dried-out regions is calculated according to the Dittus-Boelter correlation. Radiation heat transfer from the cladding to structures is also considered after dry-out respecting the viewing factor to structures.

SAS-SFR multi-bubble slug ejection model belongs to the group of separated flow approaches of computational methods developed so far to describe sodium two-phase flow behaviour. Under most conditions, the separated flow model provides a better representation of the pressure drop in a pipe flow. Although computational results suggest that the initial boiling phase is predicted by the slip model with sufficient accuracy, the phase with larger void fraction would require the application of the model of separated phases.

<figure><img src="image_3.png" /><figcaption>**Figure 1** Representation of the multiple-bubble boiling model</figcaption></figure>
NEPTUNE_CFD is a computational fluid dynamic code which uses the two-fluid model of Ishii et al. [9] for the modeling of two-phase flow. In this model, density, viscosity, volume fraction and local velocity are defined for each field in each cell. The equations solved are mass, momentum and energy balance equations for each field (continuous gas, dispersed gas, continuous liquid fields). Further, two separate models are used: one for the small and spherical bubbles and one for the larger and distorted ones as depicted in Figure 2. This three-field simulation consists in treating the small bubbles as a dispersed field where the closure laws are well known terms [9] and locating the interface of the larger bubbles which deformation cannot be taken correctly into account by non-dimensional numbers and empirical correlations. Although the achievement of the proposed approach is quite new, the concept of such three-field model has already been described in the literature. The idea of a so-called LES Hybrid approach is discussed in [10]. It is based on a spatial scale that separates the located interfaces scales from the smaller ones that are filtered and treated with a subgrid model.

In the standard approach, the gas phase is modelled by dispersed bubbles. For example, the forces exerted on dispersed bubbles are: drag force, lift force, wall lubrication, virtual mass and turbulent dispersion force. The closure laws used in the NEPTUNE_CFD code and more details on the dispersed modeling of bubbly flows can be found in [11]. These models are also implemented in commercial codes such as CFX and STAR-CCM with a formulation slightly different. Moreover, the modelling of dispersed bubbly flows has been extensively validated in previous works in the NEPTUNE_CFD code and give a reasonable agreement with experimental data.

Regarding the continuous vapor field, the same mass balance and momentum balance equations are solved but the interfacial momentum transfer term is different and represents the surface tension force and the drag force needed for the coupling between the liquid field and the gas field. This set of models is called Large Bubble Model.

Regarding the interface between the continuous liquid field and the continuous vapor field, an interface sharpening equation is solved in order to control the interface
thickness. In most CFD codes, this equation leads to the non-conservation of the mass flux
across the interface. In NEPTUNE_CFD, this equation is solved carefully in order to ensure
the mass conservation [27],[28].

<figure><img src="image_4.png" /><figcaption>**Figure 2** Illustration of the multifield approach used to simulate a two-phase flow</figcaption></figure>

The ASTEC-Na thermalhydraulic modelling is based on the two-fluid model with
separate treatment of non-condensable gases transport [12],[13]. In this model, separate
mass balances are setup for the vapor phase, the non-condensable gases and for the
liquid phase. The energy balance involves separate equations for the gas mixture phase
and the liquid phase. The momentum balance associated to the drift flux equations deals
with the mixture of liquid and gas phases. In addition, the approach implements models
for the interfacial drag between the liquid phase and the gas phase. The heat and mass
exchanges between the liquid and the vapor are modeled using the kinetic theory of
gases. In this model, the mass flux between sodium phases and the vapor-interface heat
flux are computed first and the liquid-interface is deduced by means of the interface jump
condition.
CATHARE-3 is a system code originally devoted to best estimate calculations of thermohydraulic transients in Water-Cooled Reactors [14]. It has been extended to SFR domain [15] and recent improvements have been carried out on two-phase models [16]. For the modeling of transient sodium boiling, the two-fluid model with non-condensable gases transport equations is used. A proper set of two-phase closure laws are implemented for the energy transfer (interface to liquid heat flux, interface to gas heat flux, wall to liquid heat flux, wall to gas heat flux and wall to interface heat flux for direct vaporization), the momentum transfer (entrainment fraction, interfacial friction and two-phase wall friction) and droplet diameter calculation. It did not appear necessary to model explicitly all flow pattern transitions in the CATHARE-3 code. Only two transitions are explicitly modelled and used in several closure terms: (i) transition between stratified and non-stratified flows and (ii) onset of droplet entrainment involving transition between annular and annular-mist flow. These two transitions are probably the most important as they describe the passage from a separated flow to a non-separated flow. On one side, the interfacial area is small and phases have a weak thermal and mechanical coupling. On the other side, the interfacial area is large and phases have a strong thermal and mechanical coupling.

According to the entrainment fraction, CATHARE-3 switches between different models for heat transfer, for wall friction and for interfacial friction. The wall friction term includes a coefficient which is the product of a single-phase friction coefficient and a two-phase multiplier. The standard law for the single-phase coefficient is derived from Blasius correlation but other bundle-specific laws can be used, in particular, the Rehme-correlation [18]. For the two-phase multiplier, a specific correlation, named after SENSAS experiment, has been implemented for sodium applications [16]. It is a function of the void fraction and densities and it is supposed to describe a high liquid occupation ratio on the wall up to rather high void fraction: its role is to correct the effect of an overestimated liquid film velocity actually in contact with the wall, compared to the mean liquid one. Concerning interfacial friction models for sodium applications, it has been considered that one could take benefit from the sound developments for PWRs since they cover many
flow patterns, depressurized situation as well as a large scope of hydraulic diameters.
Then the models developed for pressurized water reactors applications [17] are kept for
sodium application. Concerning heat transfer, the liquid to wall convective heat exchange
coefficient is calculated using Skupinski correlation [26]. The heat flux in case of nucleate
boiling is also computed using the Skupinski correlation. Indeed, it is considered for
sodium that the heat exchange efficiency at boiling inception is quite near the liquid
convective one (the liquid convective exchange being already very efficient and the
bubble expansion phase being relatively short in sodium compared to other fluids). The
fraction of nucleate boiling flux is kept quite low to be consistent with sodium physics
where most of the vapor is created in the bulk and very few from nucleate boiling
(because of very high sodium conductivity). For flashing, the superheat is limited to 1 or
2°C since it is supposed to be small in a SFR core considering the amount of impurities it
contains. The dry-out of the wall can be reached through two conditions: (i) The flow is in
pure gas conditions: void fraction $\alpha > \alpha_{max} = 0.99999$ and (ii) The wall temperature is too
high: $T_w > T_{sat} (P) + 200°C$; where $\alpha$ is the void fraction; $T_w$ is the wall temperature and $T_{sat}$
is the saturation temperature. The Dittus-Boelter correlation has been implemented for
the forced turbulent model of gas convection. The CATHARE-3 code also includes a point
kinetics model for reactor power calculation.

The TRACE (TRAC-RELAP Advanced Computational Engine) code is the latest in a
series of best-estimate system codes developed by the U.S. Nuclear Regulatory
Commission (NRC) for analyzing steady-state and transient thermohydraulic-neutronic
behaviour in light water reactors [19] as well as in advanced reactor systems cooled by
helium, sodium or lead-bismuth eutectic [20]. Models used include multidimensional
two-phase flow (single-phase flow for non-water fluids), non-equilibrium
thermodynamics, generalized heat transfer, re-flooding, level tracking and reactor
kinetics, using either the point kinetics model or the PARCS three-dimensional reactor
kinetics solver integrated in TRACE.
Several modifications were implemented in TRACE to improve the code's
capabilities to simulate advanced fast reactors [20]. For the modeling of transient sodium
boiling, the two-fluid model with new closure relations is used. The approach is based on
the same main conservation equations for two-phase water flow with modified closure
relations and equations-of-state. In particular, thermo-physical two-phase properties of
sodium are calculated according to [21]. The simplified two-phase sodium flow regime
map used in the modified TRACE distinguishes single-phase liquid and single-phase
vapour convection regimes with annular film boiling regime in-between as proposed in
[22]. The fraction of the perimeter æ wetted by liquid sodium is calculated by the
following correlation for all regimes: æ = $\frac{1}{1+500e^{250(\alpha-1)}}$ where α is the void fraction. This
correlation predicts æ equal to 1 for α~0.95 and in the interval from 0.95 to 1 æ is
smoothly reducing to 0. The value of 0.957 was recommended as a dryout criterion in
[22]. The æ value therefore describes a transition from the annular film boiling regime to
the single-phase vapour convection regime. Liquid-to-wall friction factor is calculated by
the Rehme model [18] for the wire-wrapped bundles and by the Churchill model [23] for
all other configurations, while vapour-to-wall friction factor is calculated by the Churchill
model [23] for all configurations. In all cases the friction factor is multiplied by the value
of æ for the liquid phase and by the value of (1 - æ) for the vapour phase. Liquid-to-wall
and vapour-to-wall heat exchange coefficients are calculated according to the model
proposed in [24] for the liquid phase and by the Dittus-Boelter model for the vapour
phase. Interfacial area density is based on the geometrical representation of the annular
film boiling pattern taking into account the value of æ as used in the SABENA code [25]
and recommended in [22]. Likewise, interfacial heat transfer and friction coefficients are
calculated as in the SABENA code [25] and as recommended in [22].
## 3.2 Remarks on the physical models

As mentioned above, two approaches are used for modeling of the transient two-phase flow resulting from sodium boiling: multi-bubble and two-fluid models. However, differences exist among the codes concerning the physical phenomena taken into account. An overview of the considered key phenomena in the codes is given in **Table 1**.

<table><caption>**Table 1** Main physical phenomena considered</caption><thead><tr><td>Code name</td><td>SAS-SFR</td><td>NEPTUNE_CFD</td><td>ASTEC-Na</td><td>CATHARE-3</td><td>TRACE</td></tr></thead><tbody><tr><td>Superheat</td><td>yes</td><td>yes</td><td>yes</td><td>yes</td><td>no</td></tr><tr><td>Nucleate boiling</td><td>no</td><td>yes</td><td>yes</td><td>yes</td><td>no</td></tr><tr><td>Voiding scheme</td><td>bubble/slug</td><td>large bubble</td><td>bubble / droplet</td><td>bubble</td><td>Porous medium</td></tr><tr><td>Re-entry rewetting</td><td>yes</td><td>yes</td><td>Yes</td><td>yes</td><td>yes</td></tr><tr><td>Film boiling</td><td>no</td><td>no</td><td>yes</td><td>yes</td><td>Yes</td></tr></tbody></table>

Despite the difference in the considered physical phenomena, both the multi-bubble and the two-fluid approach are able to satisfactorily model the sodium boiling and the related two-phase flow. However, both approaches have advantages and limitations. The multi-bubble slug ejection approach enables efficient analysis of sodium boiling both in subassembly and whole core scenario. Its limitation is mainly due to the type of the two-phase flow patterns it is able to represent as it is limited to slug flow and annular flow. Bubbly flow and dispersed flow are not reproduced. The two-fluid models take into account the two-phase flow pattern where for each phase the mass, energy and momentum conservation is separately considered. However, implementing and solving
those equations require significantly more effort. Further, the modeling of the heat, mass
and momentum exchange across the liquid-vapor interface remain a challenging task.

## 4 Model implementation

The reviewed codes implement analytical models for the sodium boiling in SFR subassembly using geometry representations which enable practical numerical solutions. The methods used for geometry representations in the codes are: one-dimensional channel, CFD, finite volume and porous medium representations. The following subsection discusses in detail the approaches for the implementation of the analytical models including the geometry representations with applied discretization schemes. In section 5, some of the limitations of the codes are highlighted.

### 4.1 Overview of analytical models

The SAS-SFR code uses one-dimensional channel geometry which can be assembled to a multi-channel to represent a SFR whole core. For each channel, the analytical model for the multi-bubble slug ejection predicts temperatures and sodium mass flow using the energy, mass, and momentum equations. It is assumed that heat transfer occurs in radial direction only. In addition, since the sodium flow is predominantly in axial direction, the flow is assumed to be one-dimensional. The mass and momentum equations are solved for the mass flow rate and pressure along a variable mesh spacing axial grid, and at each axial mesh segment, the energy equation is solved along a radial mesh for the component temperatures. The radial mesh structure can be made fine enough to model a radial temperature profile in the fuel pin, cladding, and structure as depicted in Figure 3.
<figure><img src="image_5.png" /><figcaption>**Figure 3** Axial Node Structure used in SAS-SFR code</figcaption></figure>

The description of the momentum conservation equation for the liquid flow is expressed with the flow rate instead of mass flux in order to take into account of the flow area variation for a node in the numerical computation. The frictional pressure drop is calculated using correlations.

In the NEPTUNE_CFD code, the flow geometry is represented using fine CFD meshes resolving the sub-assembly structure. The sodium boiling flows are assumed to involve predominantly large bubbles with minor component of dispersed vapor flow (small bubbles). As consequence, the analytical model focuses on the representation of large interfaces. The key components of the analytical model are models for surface tension drag force and interface sharpening.
The ASTEC-Na code uses one-dimensional finite volume method for flow geometry representation. It has been extended to multi-one-dimensional geometry in the core including radial cross-flows. The employed two-fluid model for the sodium boiling uses an interface to phase heat transfer theory to compute later the net evaporation flow. It uses internally the heat fluxes to and from the phases. The mass, energy and momentum conservation equations for each phase are solved separately. The model predicts for each mesh the coolant thermal-hydraulics data including temperature, pressure, void fraction, mass flow rate and velocity. Core structures are represented by two-dimensional meshed objects which may be cylinders for fuel pins and claddings, plates, grids or hexagonal box for the wrapper. Their geometry, temperature and thermal properties are computed at each time step, according to the implicitly computed heat fluxes exchanged with the fluid.

CATHARE-3 code uses a one-dimensional channel representation of the fuel sub-assembly for sodium applications. All fuel pins of a subassembly or a group of subassemblies are considered as a single average representative pin and the coolant channels around the fuel pins inside the hexcan wrapper are considered as a single average channel. For each one-dimensional average channel, temperatures in the structure (fuel, cladding, hexcan wrapper, etc.), and temperature, pressure and velocity in the sodium, are computed using the energy, mass, and momentum equations. It is assumed that heat transfer between the sodium and the structures and inside the structures occurs in radial direction only and the sodium flow is predominantly in axial direction. Thus, a three-dimensional calculation is reduced to two one-dimensional calculations coupled through the sodium mass flow. The wall radial mesh can be made fine enough to model a radial temperature profile in the fuel pin, cladding, and other structures. The inter-subchannel mixing is not taken into account in the one-dimensional model.
The TRACE code uses the porous medium approach where the two-fluid model is used for modelling of sodium boiling. The conservation equations for mass, momentum and energy are solved for each phase employing closure relations related to the interaction between the two phases at the liquid-gas interface and to the interaction of liquid or gas phase with the wall. The mass transfer rate is estimated using the interfacial heat transfer rates per unit area for the liquid-to-interface and gas-to-interface. The wetted fraction of the wall is estimated using empirical formula obtained from experimental data. Further models included concern the interfacial area density, interfacial friction factor, interfacial heat transfer rate, wall friction factors and wall heat transfer rate.

## 4.2 Numerical implementation

The numerical implementation of the sodium boiling models in the reviewed codes is summarized in **Table 2**. In SAS-SFR code, the equations are linearized about the values at the beginning of the time step, and fully implicit time differencing is used, leading to the requirement to solve all of the pressures and flow rates in a bubble, simultaneously. Since the non-zero elements in the matrix that must be solved are all near the diagonal, the solution by Gaussian elimination is not time-consuming. Because the equations are linearized and solved directly, no iteration is involved in the solution.

In the NEPTUNE_CFD code, a common pressure is assumed for all fields and the Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) solver is used. An iterative coupling of these equations is applied to ensure mass conservation.

In the ASTEC-Na code, the numerical method follows the finite volume technique. The space is discretized using a staggered grid with the use of the donor cell principle. The time integration is performed using a Newton's method and a fully implicit scheme is used. The Jacobian matrix inversion is based on a highly optimized Lower Upper algorithm resulting in a fast running code [12].
<table><caption>**Table 2** Numerical implementation details</caption><thead><tr><td>Code name</td><td>SAS-SFR</td><td>NEPTUNE_CFD</td><td>ASTEC-Na</td><td>CATHARE-3</td><td>TRACE</td></tr></thead><tbody><tr><td>Geometry</td><td>1-D, channel</td><td>3-D, CFD</td><td>1-D, finite volume</td><td>1-D, channel</td><td>1-3 D</td></tr><tr><td>Two-phase flow</td><td>slug / annular</td><td>large bubble, Plug / slug</td><td>annular</td><td>Annular, annular mist</td><td>annular</td></tr><tr><td>Numerical scheme</td><td>semi-implicit</td><td>SIMPLE</td><td>Newton-Raphson</td><td>fully implicit</td><td>Semi-Implicit</td></tr><tr><td>Interfacial transfer</td><td>no</td><td>tracking</td><td>tracking</td><td>no</td><td>tracking</td></tr></tbody></table>

The CATHARE-3 code discretization scheme employs the staggered mesh with scalar nodes and vector nodes and a first order upwind differencing method. This enhances the stability of the system. Mass and energy balance equations (scalar equations) are integrated on a scalar mesh from a vector node to the next one. Transport terms are evaluated at vector cell faces of a scalar mesh using the donor cell values. Velocities at scalar nodes are weighted averages of velocities at the vector nodes. Momentum balance equations (vector equations) are integrated from a scalar node to the next one and centered on a vector node. The CATHARE-3 code uses a fully implicit numerical scheme for the solution of the system of equations. The numerical scheme is not centered in time. The nonlinear system is solved using a Newton-Raphson iterative method.

The TRACE code implements sodium boiling model using the original numerical implementation which is based on the finite volume approximation [19]. For three-
dimensional application, the porous medium approach is considered. A semi-implicit
method is used for the numerical implementation.

Concerning calculated quantities, all the codes predict the coolant thermal-hydraulics data including temperature, pressure, void fraction, mass flow rate and velocity as well as temperature distributions in the structures. In addition, relevant parameters are calculated including boiling inception and onset, film boiling onset, two-phase flow duration, void fraction and pressure drops.

The SAS-SFR code calculates time and location of bulk boiling onset. After boiling onset the code provides the number of bubbles created their size, location, velocity and pressure. The volume of vapor created at any time is not directly printed by the code, and it should be calculated by the user based on the size and number of bubbles created and on the film thickness. As for the liquid film, the code provides the thickness at every time step and once the criteria selected by the user in the input are fulfilled, dry-out onset and axial front evolution is provided as well.

The NEPTUNE_CFD code is a CFD code and all variables are available at any
location and at any moment during the calculation such as void fraction, liquid and vapor
velocities, liquid and vapor temperatures, mass transfer, wall temperature.

The ASTEC-Na code provides for each time step and each mesh the coolant thermal-hydraulics data (temperature, pressure, void fraction, mass flow rate, velocity, etc.), as well as temperature of structures (fuel pin, cladding, grids, shroud, etc.). Data is available for all meshes according to the radial and axial meshing provided by the user.

The CATHARE-3 code provides for each scalar node of the hydraulic module, the pressure, the liquid enthalpy, the gas enthalpy and the void fraction of the coolant. Besides, it provides for each vector node of the hydraulic module, the liquid velocity and
the gas velocity. Finally, it provides the temperature at each radial mesh border of wall structures.

The TRACE code has the possibility to output the complete list of the calculated values (liquid, vapor and solid temperature, pressures, internal energies, densities, heat fluxes, etc.) Moreover it is possible to use the control system to calculate the needed values on the fly.

## 5 Concluding remarks

The object of the present review was to summarize the models used for transient sodium boiling phenomena in fully intact SFR subassemblies in computer codes used in benchmark activities within the ESFR-SMART project. The codes present various possibilities to model the sodium boiling and the two-phase flow in SFR subassemblies. Various differences exist between the codes with respect to physical models, assumptions, boundary conditions, geometry representations and numerical schemes. With regard to the physical model for the two-phase flow resulting from sodium boiling, essentially the multi-bubble slug ejection model and the two-fluid model are used in the codes. The geometry representations used in the codes range from one-dimensional single channel and finite volume to full three-dimensional CFD and porous media representations. In addition, the numerical implementations differ according to the chosen physical model and geometry representation.

The assessment of the accuracy, efficiency and reliability of the methods require dedicated benchmarking exercises which are being carried out within the ESFR-SMART project. Previous benchmark exercises have already shown that the codes obtain satisfactory results. However, there is a need for validations of the codes using experiments in order to improve their reliability for safety assessment.
With regard to the need for improvements, following points have been identified:

**SAS-SFR:** it computes a radially averaged sodium temperature across the fuel subassembly, it predicts boiling onset significantly later than the time observed in a fuel subassembly experiment using a limited number of heated pins. In consequence, during the time that the initially local void is expanding radially and axially in the experiment, SAS-SFR extends the prescribed flow coast down until the time that the radially averaged sodium temperature exceeds the saturation temperature (plus a few degrees of superheat) and a bubble is formed. In a reactor subassembly the difference between peak and average sodium temperature is much less pronounced. Due to the one-dimensional approach SAS-SFR can only calculate bulk boiling onset and no local boiling. To represent both axial and radial behavior of sodium would require to have access to a multi-channel representation of the fuel pin subassembly, but this is out of scope of SAS-SFR possibilities. The limitation of the multi-bubble slug ejection model concerns the type of two-phase flow patterns able to represent, it is limited to slug flow and annular flow. Bubbly flow and dispersed flow are not reproduced. Although it is a limitation, it was observed experimentally that in sodium fast reactors those are mostly the two-phase patterns taking place.

**NEPTUNE_CFD:** The models for dispersed phase are dedicated for water at high pressure and are not validated for sodium boiling less than 1 bar. But, the models for large bubbles come from theoretical considerations and can be considered as valid for sodium boiling flows. In sodium boiling flows, one can reasonably assume that large bubbles will be of first order when compared to disperse vapor flow (small bubbles). Nevertheless, the calculation of heat and mass exchange across the liquid-vapor interface needs to control efficiently the interface thickness. As consequence, the interface sharpening equation has to be solved. Unfortunately, the code is currently unable to solve this kind of equation whatever the type of cells: the cells have to be as similar as possible to cubic shape. This point is particularly uncomfortable for industrial geometries and is currently the main
difficulty for the simulation experiments. In addition, the numerical treatment of the interface sharpening equation requires further improvement.

ASTEC-Na: the code is able to correctly reproduce the apparition of the two-phase regime. However, after boiling onset, the expansion of the gas zone looks very fast compared to experimental observations in multi-pin devices, leading to an early dryout and a late contraction of the bubble after power shutdown. Several modifications of the code have been tested to correct this behaviour (modification of the flashing time constant, introduction of a wall condensation model), but without improvement of the result. It has to be mentioned that running the code with the 6 equations model (two momentum equations) leads to the same result. The reason for this too fast vaporization seems to be related to the interface area, but further work would be needed to confirm that and to correct this behaviour.

CATHARE-3: Due to the one-dimensional approach chosen for sodium applications, CATHARE-3 can only calculate bulk boiling onset and no local boiling. To represent both axial and radial behavior of sodium would require having access to a multi-channel representation of the fuel pin subassembly, but this is currently out of scope of CATHARE-3 possibilities for sodium applications.

Vaporization, condensation and flashing models are today the same as the one used for water applications with minor modifications. Those modifications are only made of common sense judgment. A further validation work has to be carried out based on separate effect tests (for example, depressurization tests). Besides, the formulation of dry-out criteria for sodium application in the CATHARE-3 code has not been evaluated in detail, and needs to be further validated. Moreover, the film boiling model, which is based on water needs, should be removed for sodium applications as film boiling is considered as unlikely (or very limited) in SFR.

Finally, the droplet entrainment formulation should be revised according to sodium specificities along with the friction model for high quality range.
TRACE: the main work on the current sodium boiling model focuses on improvement of the numerical stability for different transient conditions. In a number of cases especially characterized by the high void fraction and by the vapour condensation, the current model is not able to converge. Apart from that, the main direction of the code's improvement is linked to the validation procedure. The increase of the validation database will allow improving the physical models described above.

The performance of the codes is not discussed in the present paper. A dedicated upcoming paper will address this issue including compilations of benchmarking of the codes using experimental data.

## ACKNOWLEDGMENTS

The research leading to these results has received funding from the Euratom research and training program 2014-2018 under grant agreement No 754501.

## Acronyms and Abbreviations

ASTEC-Na = Accident Source Term Evaluation Code-Na

CATHARE = Code for Analysis of THERMALhydraulics during an Accident of Reactor and safety Evaluation

CEA = Commissariat de l'Energie Atomique

CFD = Computational Fluid Dynamics

CP-ESFR = Collaborative Project European Sodium Fast Reactor

EDF = Electricité de France

ESFR = European Sodium Fast Reactor
ESFR- = European Sodium Fast Reactor Safety Measures Assessment and Research

SMART Tools

FP7 = Seventh Framework Program

IRSN = Institut de Radioprotection et de Sûreté Nucléaire

KIT = Karlsruhe Institute of Technology

PSI = Paul Scherrer Institute

SABENA = SubAssembly Boiling Evolution Numerical Analysis

SFR = Sodium Cooled Fast Reactor

Tsat = Saturation Temperature

Tw = Wall Temperature

TRACE = TRAC/RELAP Advanced Computational Engine

ULOF = Unprotected Loss of Flow

ULOHS = Unprotected Loss of Heat Sink

UTOP = Unprotected Transient Over Power

Tsat = Saturation Temperature

Tw = Wall Temperature

1-D = One-dimensional

3-D = Three-dimensional

α = Void fraction
REFERENCES

[1] G.L. Fiorini, A. Vasile, 2011. European Commission – 7th Framework Programme: The Collaborative Project on European Sodium Fast Reactor (CP ESFR), Nuclear Engineering and Design. Vol. 241, pp. 3461–3469.
https://doi.org/10.1016/j.nucengdes.2011.01.052

[2] K. Mikityuk, E. Girardi, J. Krepel, E. Bubelis, E. Fridman, A. Rineiski, N. Girault, F. Payot, L. Buligins, G. Gerbeth, N. Chauvin, C. Latge, J.-C. Garnier, 2017. ESFR-SMART: new Horizon-2020 project on SFR safety, in proc. IAEA FR2017.

[3] U.S. Department of Energy Nuclear Energy Research Advisory Committee and the Generation IV International Forum, 2002. A Technology Roadmap for Generation IV Nuclear Energy Systems, OECD/NEA.

[4] A. E. Waltar and A. B. Reynolds, 1981. Fast Breeder Reactors, Pergamon Press.

[5] H. Fauske, 1972-1973. Nucleation of Liquid Sodium in Fast Reactors, Reactor Technol., 15, 4, 278-302.

[6] J.M. Seiler, D. Juhel and P. Dufour, 2010. Sodium Boiling Stabilization in a Fast Breeder Subassembly during an Unprotected Loss of Flow Accident, Nucl. Eng. Des., 240, 3329 – 3335, https://doi.org/10.1016/j.nucengdes.2010.07.001

[7] R. Kruessmann, A. Ponomarev, W. Pfrang, D. Struwe, J.Champigny, B.Carluec, D.Schmitt, D.Verwaerde, 2015. Assessment of SFR reactor safety issues: Part II: Analysis results of ULOF transients imposed on a variety of different innovative core designs with SAS-SFR, Nucl. Eng. Des, 285, pp. 263-28, https://doi.org/10.1016/j.nucengdes.2014.11.037

[8] A. M. Tentner, G. Birgersson, J. E. Cahalan, F. E. Dunn, Kalimullah, K. J. Miles, 1987. SAS4A A Computer Model for the Analysis of Hypothetical Core Disruptive Accidents in Liquid Metal Reactors. Paper submitted for presentation at the 1987 Eastern Computer Simulation Conference, Orlando, FLA.
[9] M. Ishii and T. Hibiki, 2006. Thermo-fluid dynamics of two-phase flows, Springer Science & Business Media.

[10] D. Bestion, 2012. Applicability of two-phase CFD to nuclear reactor thermalhydraulics and elaboration of best practice guidelines, Nucl. Eng. Des. 253, 311-321, https://doi.org/10.1016/j.nucengdes.2011.08.068

[11] S. Mimouni, J. Lavieville, N. Mérigoux, 2017. Combined evaluation of bubble dynamics, polydispersion model and turbulence modeling for adiabatic two-phase flow, Nucl. Eng. Des., 321, 57-68, https://doi.org/10.1016/j.nucengdes.2017.03.041

[12] P. Chatelard, S. Belon, L. Bosland, L. Carénini, O. Coindreau, F. Cousin, C.Marchetto, H.Nowack, L. Piara, L. Chailana, et al., 2016, Main modelling features of the ASTEC V2.1 major version, Ann. Nucl. Ener. 93, 83-93, https://doi.org/10.1016/j.anucene.2015.12.026.

[13] Gómez-García-Toraño and L. Laborde, 2019. Validation of the CESAR friction models of the ASTECV21 code based on Moby Dick experiments. J. Nucl. Eng. Radiat. Sci. 5, 1-9, https://doi.org/10.1115/1.4042119.

[14] D. Bestion and G. Geffraye, 2001. The CATHARE code. Operational Practice of Nuclear power Plants, Budapest University, Hungary.

[15] M. Robert, M. Farvacque, M. Parent, and B. Faydide, 2003. CATHARE 2 V2.5: A fully validated CATHARE version for various applications. In 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-10) Seoul, Korea.

[16] M. Anderhuber, J. Perez-Manes, N. Alpyet, 2017, Validation of CATHARE-3 system code for sodium two-phase flow application: R&D path and simulation of key-tests from the SIENA program. NURETH-17, Xi'an, China.

[17] D. Bestion, 1990. The physical closure laws in the CATHARE code. Nuclear Engineering and Design, 124, 229-245, https://doi.org/10.1016/0029-5493(90)90294-8
[18] K. Rehme, 1972. Pressure Drop Performance of Rod Bundles in Hexagonal Arrangements, International Journal of Heat and Mass Transfer. Vol 15, pp. 2499-2517, https://doi.org/10.1016/0017-9310(72)90143-3

[19] Nuclear Regulatory Commission, 2013. TRACE V5.0. Theory Manual. Field Equations, Solution Methods, and Physical Models, Office for Nuclear Regulatory Research, Washington.

[20] K. Mikityuk, S. Pelloni, P. Coddington, E. Bubelis, R. Chawla, 2005. FAST: An advanced code system for fast reactor transient analysis, Annals of Nuclear Energy 32, 1613-1631, https://doi.org/10.1016/j.anucene.2005.06.002

[21] J. K. Fink and L. Leibowitz, 1995. Thermodynamic and transport properties of sodium liquid and vapor, Argonne National Laboratory, ANL/RE-95/2.

[22] A. Chenu, 2011. Single- and Two-Phase Flow Modelling for Coupled Neutronics/Thermal-Hydraulics Transient Analysis of Advanced Sodium-Cooled Fast Reactors, EPFL Thesis, p.42, http://dx.doi.org/10.5075/epfl-thesis-5172.

[23] S.W. Churchill, 1977. Friction Factor Equations Spans All Fluid-Flow Regimes, Chemical Eng., November, 91-92.

[24] K. Mikityuk, 2009. Heat transfer to liquid metal: Review of data and correlations for tube bundles, Nucl. Eng. Des., vol. 239, 680-687, https://doi.org/10.1016/j.nucengdes.2008.12.014

[25] H. Ninokata and T. Okano, 1990. SABENA: Subassembly boiling evolution numerical analysis, Nucl. Eng. Des., vol. 120:pp. 349 – 367, https://doi.org/10.1016/0029-5493(90)90386-C

[26] Skupinski, E., Tortel, J. et Vautrey L., 1964. Détermination des coefficients de convection d'un alliage sodium-potassium dans un tube circulaire. Int. J. Heat Mass Transf., vol 8, 937-951, https://doi.org/10.1016/0017-9310(65)90077-3
[27] Mimouni, S., Guingo, M. and Lavieville, J., 2017. Assessment of RANS at low Prandtl number and simulation of sodium boiling flows with a CMFD code, Nucl. Eng. Des., vol. 312, 294–302, https://doi.org/10.1016/j.nucengdes.2016.07.006

[28] Mimouni, S., Fleau, S. and Vincent, S., 2017. CFD calculations of flow pattern maps and LES of multiphase flows. Nucl. Eng. Des., vol 321, 118–131, https://doi.org/10.1016/j.nucengdes.2016.12.009
Figure Captions List

Fig. 1 Representation of the multiple-bubble boiling model

Fig. 2 Illustration of the multifield approach used to simulate a two-phase flow

Fig. 3 Axial Node Structure used in SAS-SFR code
Table Caption List

Table 1 Main physical phenomena considered

Table 2 Numerical implementation details