# Assessment of CTF Boiling Transition and Critical Heat Flux Modeling Capabilities Using the OECD/NRC BFBT and PSBT Benchmark Databases

Authors: Maria Avramova; Diana Cuervo
Journal: Science and Technology of Nuclear Installations
 (2013)

---

## Abstract

Over the last few years, the Pennsylvania State University (PSU) under the sponsorship of the US Nuclear Regulatory Commission (NRC) has prepared, organized, conducted, and summarized two international benchmarks based on the NUPEC data—the OECD/NRC Full-Size Fine-Mesh Bundle Test (BFBT) Benchmark and the OECD/NRC PWR Sub-Channel and Bundle Test (PSBT) Benchmark. The benchmarks’ activities have been conducted in cooperation with the Nuclear Energy Agency/Organization for Economic Co-operation and Development (NEA/OECD) and the Japan Nuclear Energy Safety (JNES) Organization.
This paper presents an application of the joint Penn State University/Technical University of Madrid (UPM) version of the well-known sub-channel code COBRA-TF (Coolant Boiling in Rod Array-Two Fluid), namely, CTF, to the steady state critical power and departure from nucleate boiling (DNB) exercises of the OECD/NRC BFBT and PSBT benchmarks. The goal is two-fold: firstly, to assess these models and to examine their strengths and weaknesses; and secondly, to identify the areas for improvement.


---

## Body

## 1. Introduction

### 1.1. Motivation for the Work

The increased use and importance of detailed reactor core descriptions for light water reactor (LWR) safety analysis and coupled local neutronics/thermal-hydraulics evaluations requires the use of advanced two-phase thermal-hydraulic codes. These codes must be extensively validated against full-scale high-quality experimental data. In that sense, the international OECD/NRC Boiling Water Reactor Full-Size-Fine-Mesh Bundle Test Benchmark [1] and the OECD/NRC PWR Sub-Channel and Bundle Tests Benchmark [2] provide an excellent opportunity for validation of models for critical power and departure from nucleate boiling.The OECD/NRC BFBT and PSBT benchmarks were established to provide test beds for assessing the capabilities of various thermal-hydraulic subchannel, system, and computational fluid dynamics (CFD) codes and to encourage advancements in the analysis of fluid flow in rod bundles. The aim was to improve the reliability of the nuclear reactor safety margin evaluations. The benchmarks are based on one of the most valuable databases identified for the thermal-hydraulics modelling, which was developed by the Nuclear Power Engineering Corporation in Japan.This paper presents the results obtained with the thermal-hydraulic code CTF [3] for the Exercise II-1 (steady-state critical power) of the OECD/NRC BFBT benchmark and Exercise II-1 (steady-state departure from nucleate boiling) of the OECD/NRC PSBT benchmark. Transient CHF simulations are not considered in this study. Assessments of the code capabilities for prediction of the steady-state and transient void distribution along with uncertainty and sensitivity analysis were previously performed and published [4–6].Although CTF was already subjected to an extensive verification and validation program and applied to a variety of LWR steady-state and transient simulations [7, 8], the code assessment to CHF experiments was limited to single-tube geometries. Moreover, the CTF heat transfer package has not been further developed and it is essentially the same as in the original COBRA-TF code from the early 1980s. A brief description of the CTF flow regimes and heat transfer package is given in the following section.The goal of the presented work is twofold: firstly, to assess these models and to examine their strengths and weaknesses; and secondly, to identify the areas for improvement.


### 1.2. Background on the Thermal-Hydraulic Subchannel Code CTF

CTF, a version of the well-known subchannel code COBRA-TF, is being maintained by the Reactor Dynamics and Fuel Management Group (RDFMG) at the Pennsylvania State University (PSU) in cooperation with several partners including the Technical University of Madrid (UPM) in Spain. The original version of COBRA-TF was developed at the Pacific Northwest Laboratory as a part of the COBRA/TRAC thermal-hydraulic code under the sponsorship of US NRC [9]. Since then, various academic and industrial organizations have adapted, developed, and modified the code in many directions. The code is worldwide used for academic and general research purposes. The code version used at PSU originates from the COBRA-TF version modified during the FLECHT SEASET program [10]. In parallel to the code utilization to teach and train students in the area of nuclear reactor thermal-hydraulic safety analyses at PSU and UPM, the theoretical models and numerics of CTF were substantially improved [11, 12].CTF is a transient code based on a separated flow representation of the two-phase flow. The two-fluid formulation, often used in thermal-hydraulic codes, separates the conservation equations of mass, energy, and momentum to vapor and liquid. CTF extends this treatment to three fields: vapor, continuous liquid, and entrained liquid droplets, which results in a set of nine time-averaged conservation equations. The conservation equations for each of the three fields and for heat transfer from and within the solid structure in contact with the fluid are solved using a semi-implicit, finite-difference numerical technique on an Eulerian mesh, where time intervals are assumed to be long enough to smooth out the random fluctuations in the multiphase flow, but short enough to preserve any gross flow unsteadiness. The code is able to handle both pre- and post-CHF flow regimes and is capable of calculating reverse flow, counter flow, and cross-flow situations. The code is developed for use with either three-dimensional (3D) Cartesian or subchannel coordinates and, therefore, it features extremely flexible nodding for both the thermal-hydraulic and the heat-transfer solution. This flexibility allows a fully 3D treatment in geometries amenable to description in a Cartesian coordinate system.It is worth mentioning here that CTF is being used at both universities for coupling with different 3D neutron-kinetics codes. At UPM, the code is part of the COBAYA3 [13] system of codes for multiphysics and multiscale core calculations. The code has been coupled with the ANDES nodal scale diffusion code [14] for nodal calculations and with the COBAYA3K pin-by-pin diffusion code [15] for fine-mesh calculations. Both systems of coupled codes are part of a multiscale calculation methodology based on a subdomain decomposition of the core for fast pin-by-pin diffusion calculations of the whole core [16]. Validation of this system is being carried out [17]. At PSU, a 3D neutron kinetics module was implemented into CTF by a serial integration coupling to the PSU NEM code. The new PSU coupled code system was named CTF/NEM [18].


## 1.1. Motivation for the Work

The increased use and importance of detailed reactor core descriptions for light water reactor (LWR) safety analysis and coupled local neutronics/thermal-hydraulics evaluations requires the use of advanced two-phase thermal-hydraulic codes. These codes must be extensively validated against full-scale high-quality experimental data. In that sense, the international OECD/NRC Boiling Water Reactor Full-Size-Fine-Mesh Bundle Test Benchmark [1] and the OECD/NRC PWR Sub-Channel and Bundle Tests Benchmark [2] provide an excellent opportunity for validation of models for critical power and departure from nucleate boiling.The OECD/NRC BFBT and PSBT benchmarks were established to provide test beds for assessing the capabilities of various thermal-hydraulic subchannel, system, and computational fluid dynamics (CFD) codes and to encourage advancements in the analysis of fluid flow in rod bundles. The aim was to improve the reliability of the nuclear reactor safety margin evaluations. The benchmarks are based on one of the most valuable databases identified for the thermal-hydraulics modelling, which was developed by the Nuclear Power Engineering Corporation in Japan.This paper presents the results obtained with the thermal-hydraulic code CTF [3] for the Exercise II-1 (steady-state critical power) of the OECD/NRC BFBT benchmark and Exercise II-1 (steady-state departure from nucleate boiling) of the OECD/NRC PSBT benchmark. Transient CHF simulations are not considered in this study. Assessments of the code capabilities for prediction of the steady-state and transient void distribution along with uncertainty and sensitivity analysis were previously performed and published [4–6].Although CTF was already subjected to an extensive verification and validation program and applied to a variety of LWR steady-state and transient simulations [7, 8], the code assessment to CHF experiments was limited to single-tube geometries. Moreover, the CTF heat transfer package has not been further developed and it is essentially the same as in the original COBRA-TF code from the early 1980s. A brief description of the CTF flow regimes and heat transfer package is given in the following section.The goal of the presented work is twofold: firstly, to assess these models and to examine their strengths and weaknesses; and secondly, to identify the areas for improvement.


## 1.2. Background on the Thermal-Hydraulic Subchannel Code CTF

CTF, a version of the well-known subchannel code COBRA-TF, is being maintained by the Reactor Dynamics and Fuel Management Group (RDFMG) at the Pennsylvania State University (PSU) in cooperation with several partners including the Technical University of Madrid (UPM) in Spain. The original version of COBRA-TF was developed at the Pacific Northwest Laboratory as a part of the COBRA/TRAC thermal-hydraulic code under the sponsorship of US NRC [9]. Since then, various academic and industrial organizations have adapted, developed, and modified the code in many directions. The code is worldwide used for academic and general research purposes. The code version used at PSU originates from the COBRA-TF version modified during the FLECHT SEASET program [10]. In parallel to the code utilization to teach and train students in the area of nuclear reactor thermal-hydraulic safety analyses at PSU and UPM, the theoretical models and numerics of CTF were substantially improved [11, 12].CTF is a transient code based on a separated flow representation of the two-phase flow. The two-fluid formulation, often used in thermal-hydraulic codes, separates the conservation equations of mass, energy, and momentum to vapor and liquid. CTF extends this treatment to three fields: vapor, continuous liquid, and entrained liquid droplets, which results in a set of nine time-averaged conservation equations. The conservation equations for each of the three fields and for heat transfer from and within the solid structure in contact with the fluid are solved using a semi-implicit, finite-difference numerical technique on an Eulerian mesh, where time intervals are assumed to be long enough to smooth out the random fluctuations in the multiphase flow, but short enough to preserve any gross flow unsteadiness. The code is able to handle both pre- and post-CHF flow regimes and is capable of calculating reverse flow, counter flow, and cross-flow situations. The code is developed for use with either three-dimensional (3D) Cartesian or subchannel coordinates and, therefore, it features extremely flexible nodding for both the thermal-hydraulic and the heat-transfer solution. This flexibility allows a fully 3D treatment in geometries amenable to description in a Cartesian coordinate system.It is worth mentioning here that CTF is being used at both universities for coupling with different 3D neutron-kinetics codes. At UPM, the code is part of the COBAYA3 [13] system of codes for multiphysics and multiscale core calculations. The code has been coupled with the ANDES nodal scale diffusion code [14] for nodal calculations and with the COBAYA3K pin-by-pin diffusion code [15] for fine-mesh calculations. Both systems of coupled codes are part of a multiscale calculation methodology based on a subdomain decomposition of the core for fast pin-by-pin diffusion calculations of the whole core [16]. Validation of this system is being carried out [17]. At PSU, a 3D neutron kinetics module was implemented into CTF by a serial integration coupling to the PSU NEM code. The new PSU coupled code system was named CTF/NEM [18].


## 2. Overview of the CTF Flow Regimes and Heat Transfer Package

### 2.1. CTF Flow Regime Maps

As in the earlier COBRA-TF versions [9, 10], CTF contains two different types of flow regime maps: “normal wall” map and “hot wall” map. The normal wall map is used when the maximum wall surface temperature in a given computational mesh cell is below the critical heat flux temperature and is, thus, expected to be fully wetted. The hot wall map, on the other hand, is selected when the maximum wall surface temperature exceeds the critical heat flux temperature. The critical heat flux temperature is assumed to be well approximated by a wall superheat of 41.7°C or higher. There is a CHF temperature upper limit of 374.1°C, which corresponds to the critical temperature of water.The normal wall flow regime map includes the following flow regimes: small bubble; small-to-large bubble (slug); churn/turbulent; and annular/mist. If the maximum wall temperature exceeds the CHF temperature, a whole new range of significantly different flow regimes become possible since the liquid can only partially wet the wall. This occurs in PWRs during accident conditions, like the blowdown phase of a large-break loss of coolant accident (LOCA). The flow regimes recognized in the hot wall map are inverted annular flow; inverted slug flow; dispersed droplet flow; falling flow; and top deluge flow.A detailed description of the CTF flow regime maps and transition logic can be found in [19]. Since the physical models used in the numerical solution must be defined for each mesh cell, the flow regime must be determined from fluid properties and flow conditions within each cell or in the immediate surrounding cells. Once the flow regime of the mesh cell is correctly identified, the appropriate models can be chosen for calculation of the closure terms such as the interfacial heat transfer, the interfacial drag, and the wall drag.It has to be mentioned here that the code was developed for vertical two-phase flow in rod bundle geometries and, therefore, horizontal flow regimes were not considered; however, an implementation of a horizontal flow regime map is currently being carried out at PSU [20].


### 2.2. CTF Heat Transfer Package

The heat transfer models in CTF determine the material heat release rates and the temperature response of the fuel rod and structural components of LWRs during operating and transient conditions. All heat transfer calculations are performed at the beginning of each time step before the hydraulic solution. Heat transfer coefficients based on the liquid conditions in the previous time step are used to advance the material conduction solution. The resultant heat release rates are explicitly coupled to the hydrodynamic solution as source terms in the fluid energy equations. The CTF heat transfer package consists of a library of heat transfer coefficients and a selection logic. Together these produce a boiling curve that is used to determine the phasic heat fluxes. CTF recognizes the following heat transfer regimes: single-phase liquid convection; single-phase vapor convection; subcooled nucleate boiling; saturated nucleate boiling; transition boiling; inverted annular film boiling; dispersed droplet film boiling; dispersed droplet deposition heat transfer.CTF selects the appropriate heat transfer regime after evaluating several criteria. If the heated surface temperature is 0.1°C less than the critical heat flux temperature, a pre-CHF heat transfer regime will be selected (i.e., single-phase liquid convection, subcooled nucleate boiling, or saturated nucleate boiling). The single-phase vapor convection regime will be selected if the void fraction is above 0.999, regardless of the heated surface temperature. On the other hand, if the heated surface temperature is greater than or equal to 0.1°C above the critical heat flux temperature, then one of the post-CHF heat transfer regimes will be selected (i.e., transition boiling, inverted annular film boiling, dispersed droplet film boiling, or dispersed droplet deposition heat transfer). A further distinction is made by selecting the transition boiling regime if the heated wall temperature is less than the minimum film boiling temperature and one of the other post-CHF regimes if the heated wall temperature is larger than the minimum film boiling temperature.Prior to determining the heat transfer regime in the aforementioned manner, it is necessary to first determine the CHF temperature and the minimum film boiling temperature. The CHF temperature is calculated iteratively using the previously calculated CHF. The minimum and maximum boundaries for CHF temperature are 11°C over the liquid saturation temperature (minimum) and the maximum of 111°C over the saturation temperature or the critical temperature of water (maximum). For unheated conductors, the minimum film boiling temperature is set to a constant value of 482°C. For heated structures, the minimum film boiling temperature is set to a minimum value of 482°C for void fractions less than 80% and to a minimum value of 371°C for void fractions equal to or greater than 80%. The minimum film boiling temperature may be much higher, though, and is evaluated using two different methods: first, it is calculated assuming it equals the wall temperature that results in an instantaneous contact temperature equal to the homogeneous nucleation temperature; and second, the Henry’s modification of the Berenson correlation is used [19].If the mesh cell contains a hot wall, CHF is calculated for the annular film dryout region. If the mesh cell does not contain a hot wall, the void fraction is checked because the boiling regime could still be in annular film dryout. If the void fraction is higher than 90%, the annular film dryout CHF is calculated; and if the void fraction is higher than 99%, the annular film dryout CHF is multiplied by a ramping factor to reduce CHF down to 20% of its calculated value as the liquid fraction approaches 0.005. If the void fraction is less than 90%, then CHF is calculated for the forced-convection boiling regime.The heat transfer regime selection logic and the correlations used in each regime are briefly discussed next.Single-Phase Vapor. The maximum of the Dittus-Boelter turbulent convection correlation [21]; the FLECHT SEASET 161-rod steam cooling correlation [22]; and a laminar flow Nusselt number is used. For single-phase convection to vapor, all vapor properties are evaluated at the film temperature.Single-Phase Liquid. Convection to single-phase liquid is computed as the larger of either the Dittus-Boelter turbulent convection correlation or the laminar flow with a limit Nusselt number equal to 7.86 [23].Nucleate Boiling. When the temperature is greater than saturation but less than the critical heat flux temperature and liquid is present on the wall, the Chen nucleate boiling correlation [24] is used. The Chen correlation applies to both the saturated nucleate boiling region and the two-phase forced convection evaporation region. It automatically makes the transition to single-phase convection at low wall superheat and pool boiling at low flow rate. The Chen correlation assumes a superposition of a forced-convection correlation (Dittus-Boelter type) and a pool boiling equation (Forster-Zuber).Subcooled Nucleate Boiling. An extension of the Chen nucleate boiling correlation into the subcooled region is used for subcooled nucleate boiling. During the subcooled boiling, vapor generation occurs and a significant void fraction may exist despite the presence of subcooled water. The processes of interest in this regime are forced convection to liquid, vapor generation at the wall, condensation near the wall, and bulk condensation (subcooled liquid core).Critical Heat Flux and Transition Boiling Regime. Three critical heat flux regimes are considered—pool boiling, forced convection DNB, and annular film dryout. Pool boiling DNB is selected when the mass flux is low (below 30 g/cm2-sec) and the flow regime is not annular film flow. The pool boiling heat flux is given by Griffith’s [25] modification of the Zuber [26] equation. The critical heat flux in this region is chosen as the larger of the Griffith’s modification and the forced convection DNB heat flux at a mass flux of 30 g/cm2-sec. Forced-convection DNB is considered when the mass flux is greater than 30 g/cm2-sec and the flow regime is not annular film flow. The critical heat flux is given by the Biasi correlation [27], which consists of two equations: one for low-quality CHF and one for high-quality CHF. The critical heat flux is defined as the maximum of the two equations. If annular flow exists, the departure from nucleate boiling is caused by annular film dryout. In this regime, the critical heat flux is not limited by a correlation, but rather forced convection vaporization exists until the film dries out. Film dryout is a complex function of the film flow rate, the applied heat flux, and the entrainment-deentrainment rate. Film dryout is determined by the solution of the hydrodynamic equations. A value of 41.7°C wall superheat is selected to be a CHF point for annular film dryout and the CHF is set to that given by the Zuber equation. The critical heat flux temperature is defined using an iterative procedure to determine the wall temperature at which the heat flux from the Chen nucleate boiling correlation is equal to the CHF.The transition boiling regime is bounded by the CHF point (below which the wall is continuously wetted and nucleate boiling exists) and the minimum stable film boiling point (above which the liquid cannot wet the wall and film boiling exists). It is assumed that the minimum film boiling temperature is the wall temperature that results in an instantaneous contact temperature equal to the homogeneous nucleation temperature. In addition, the minimum film boiling temperature is restricted to varies between 427°C and 649°C.CTF employs a simple additive scheme for heat transfer beyond the critical heat flux temperature. It is assumed that the transition boiling heat transfer is composed of both liquid contact (wet wall) and film boiling (dry wall).Heat transfer in the film boiling region is assumed to result either from dispersed flow film boiling or from inverted annular film boiling.Dispersed Flow Film Boiling. Dispersed flow film boiling is selected if the void fraction is greater than 80%. It is treated by a “two-step” method where the dominant heat transfer mode is forced convection to superheated steam. The steam superheat is determined by the interfacial heat transfer rate to the entrained droplets as part of the hydrodynamic solution. Heat fluxes due to wall-droplet radiation and droplet impingement are superimposed upon the vapor convective heat flux.Inverted Annular Film Boiling. When the void fraction is less than 60%, inverted annular film boiling is assumed to occur. The heat flux for this regime is computed from the larger of either dispersed film boiling heat flux as defined above or the value from the modified Bromley correlation [28]. At intermediate void fractions (60%–80%), the heat flux is interpolated between the value for inverted annular and dispersed flow film boiling.


## 2.1. CTF Flow Regime Maps

As in the earlier COBRA-TF versions [9, 10], CTF contains two different types of flow regime maps: “normal wall” map and “hot wall” map. The normal wall map is used when the maximum wall surface temperature in a given computational mesh cell is below the critical heat flux temperature and is, thus, expected to be fully wetted. The hot wall map, on the other hand, is selected when the maximum wall surface temperature exceeds the critical heat flux temperature. The critical heat flux temperature is assumed to be well approximated by a wall superheat of 41.7°C or higher. There is a CHF temperature upper limit of 374.1°C, which corresponds to the critical temperature of water.The normal wall flow regime map includes the following flow regimes: small bubble; small-to-large bubble (slug); churn/turbulent; and annular/mist. If the maximum wall temperature exceeds the CHF temperature, a whole new range of significantly different flow regimes become possible since the liquid can only partially wet the wall. This occurs in PWRs during accident conditions, like the blowdown phase of a large-break loss of coolant accident (LOCA). The flow regimes recognized in the hot wall map are inverted annular flow; inverted slug flow; dispersed droplet flow; falling flow; and top deluge flow.A detailed description of the CTF flow regime maps and transition logic can be found in [19]. Since the physical models used in the numerical solution must be defined for each mesh cell, the flow regime must be determined from fluid properties and flow conditions within each cell or in the immediate surrounding cells. Once the flow regime of the mesh cell is correctly identified, the appropriate models can be chosen for calculation of the closure terms such as the interfacial heat transfer, the interfacial drag, and the wall drag.It has to be mentioned here that the code was developed for vertical two-phase flow in rod bundle geometries and, therefore, horizontal flow regimes were not considered; however, an implementation of a horizontal flow regime map is currently being carried out at PSU [20].


## 2.2. CTF Heat Transfer Package

The heat transfer models in CTF determine the material heat release rates and the temperature response of the fuel rod and structural components of LWRs during operating and transient conditions. All heat transfer calculations are performed at the beginning of each time step before the hydraulic solution. Heat transfer coefficients based on the liquid conditions in the previous time step are used to advance the material conduction solution. The resultant heat release rates are explicitly coupled to the hydrodynamic solution as source terms in the fluid energy equations. The CTF heat transfer package consists of a library of heat transfer coefficients and a selection logic. Together these produce a boiling curve that is used to determine the phasic heat fluxes. CTF recognizes the following heat transfer regimes: single-phase liquid convection; single-phase vapor convection; subcooled nucleate boiling; saturated nucleate boiling; transition boiling; inverted annular film boiling; dispersed droplet film boiling; dispersed droplet deposition heat transfer.CTF selects the appropriate heat transfer regime after evaluating several criteria. If the heated surface temperature is 0.1°C less than the critical heat flux temperature, a pre-CHF heat transfer regime will be selected (i.e., single-phase liquid convection, subcooled nucleate boiling, or saturated nucleate boiling). The single-phase vapor convection regime will be selected if the void fraction is above 0.999, regardless of the heated surface temperature. On the other hand, if the heated surface temperature is greater than or equal to 0.1°C above the critical heat flux temperature, then one of the post-CHF heat transfer regimes will be selected (i.e., transition boiling, inverted annular film boiling, dispersed droplet film boiling, or dispersed droplet deposition heat transfer). A further distinction is made by selecting the transition boiling regime if the heated wall temperature is less than the minimum film boiling temperature and one of the other post-CHF regimes if the heated wall temperature is larger than the minimum film boiling temperature.Prior to determining the heat transfer regime in the aforementioned manner, it is necessary to first determine the CHF temperature and the minimum film boiling temperature. The CHF temperature is calculated iteratively using the previously calculated CHF. The minimum and maximum boundaries for CHF temperature are 11°C over the liquid saturation temperature (minimum) and the maximum of 111°C over the saturation temperature or the critical temperature of water (maximum). For unheated conductors, the minimum film boiling temperature is set to a constant value of 482°C. For heated structures, the minimum film boiling temperature is set to a minimum value of 482°C for void fractions less than 80% and to a minimum value of 371°C for void fractions equal to or greater than 80%. The minimum film boiling temperature may be much higher, though, and is evaluated using two different methods: first, it is calculated assuming it equals the wall temperature that results in an instantaneous contact temperature equal to the homogeneous nucleation temperature; and second, the Henry’s modification of the Berenson correlation is used [19].If the mesh cell contains a hot wall, CHF is calculated for the annular film dryout region. If the mesh cell does not contain a hot wall, the void fraction is checked because the boiling regime could still be in annular film dryout. If the void fraction is higher than 90%, the annular film dryout CHF is calculated; and if the void fraction is higher than 99%, the annular film dryout CHF is multiplied by a ramping factor to reduce CHF down to 20% of its calculated value as the liquid fraction approaches 0.005. If the void fraction is less than 90%, then CHF is calculated for the forced-convection boiling regime.The heat transfer regime selection logic and the correlations used in each regime are briefly discussed next.Single-Phase Vapor. The maximum of the Dittus-Boelter turbulent convection correlation [21]; the FLECHT SEASET 161-rod steam cooling correlation [22]; and a laminar flow Nusselt number is used. For single-phase convection to vapor, all vapor properties are evaluated at the film temperature.Single-Phase Liquid. Convection to single-phase liquid is computed as the larger of either the Dittus-Boelter turbulent convection correlation or the laminar flow with a limit Nusselt number equal to 7.86 [23].Nucleate Boiling. When the temperature is greater than saturation but less than the critical heat flux temperature and liquid is present on the wall, the Chen nucleate boiling correlation [24] is used. The Chen correlation applies to both the saturated nucleate boiling region and the two-phase forced convection evaporation region. It automatically makes the transition to single-phase convection at low wall superheat and pool boiling at low flow rate. The Chen correlation assumes a superposition of a forced-convection correlation (Dittus-Boelter type) and a pool boiling equation (Forster-Zuber).Subcooled Nucleate Boiling. An extension of the Chen nucleate boiling correlation into the subcooled region is used for subcooled nucleate boiling. During the subcooled boiling, vapor generation occurs and a significant void fraction may exist despite the presence of subcooled water. The processes of interest in this regime are forced convection to liquid, vapor generation at the wall, condensation near the wall, and bulk condensation (subcooled liquid core).Critical Heat Flux and Transition Boiling Regime. Three critical heat flux regimes are considered—pool boiling, forced convection DNB, and annular film dryout. Pool boiling DNB is selected when the mass flux is low (below 30 g/cm2-sec) and the flow regime is not annular film flow. The pool boiling heat flux is given by Griffith’s [25] modification of the Zuber [26] equation. The critical heat flux in this region is chosen as the larger of the Griffith’s modification and the forced convection DNB heat flux at a mass flux of 30 g/cm2-sec. Forced-convection DNB is considered when the mass flux is greater than 30 g/cm2-sec and the flow regime is not annular film flow. The critical heat flux is given by the Biasi correlation [27], which consists of two equations: one for low-quality CHF and one for high-quality CHF. The critical heat flux is defined as the maximum of the two equations. If annular flow exists, the departure from nucleate boiling is caused by annular film dryout. In this regime, the critical heat flux is not limited by a correlation, but rather forced convection vaporization exists until the film dries out. Film dryout is a complex function of the film flow rate, the applied heat flux, and the entrainment-deentrainment rate. Film dryout is determined by the solution of the hydrodynamic equations. A value of 41.7°C wall superheat is selected to be a CHF point for annular film dryout and the CHF is set to that given by the Zuber equation. The critical heat flux temperature is defined using an iterative procedure to determine the wall temperature at which the heat flux from the Chen nucleate boiling correlation is equal to the CHF.The transition boiling regime is bounded by the CHF point (below which the wall is continuously wetted and nucleate boiling exists) and the minimum stable film boiling point (above which the liquid cannot wet the wall and film boiling exists). It is assumed that the minimum film boiling temperature is the wall temperature that results in an instantaneous contact temperature equal to the homogeneous nucleation temperature. In addition, the minimum film boiling temperature is restricted to varies between 427°C and 649°C.CTF employs a simple additive scheme for heat transfer beyond the critical heat flux temperature. It is assumed that the transition boiling heat transfer is composed of both liquid contact (wet wall) and film boiling (dry wall).Heat transfer in the film boiling region is assumed to result either from dispersed flow film boiling or from inverted annular film boiling.Dispersed Flow Film Boiling. Dispersed flow film boiling is selected if the void fraction is greater than 80%. It is treated by a “two-step” method where the dominant heat transfer mode is forced convection to superheated steam. The steam superheat is determined by the interfacial heat transfer rate to the entrained droplets as part of the hydrodynamic solution. Heat fluxes due to wall-droplet radiation and droplet impingement are superimposed upon the vapor convective heat flux.Inverted Annular Film Boiling. When the void fraction is less than 60%, inverted annular film boiling is assumed to occur. The heat flux for this regime is computed from the larger of either dispersed film boiling heat flux as defined above or the value from the modified Bromley correlation [28]. At intermediate void fractions (60%–80%), the heat flux is interpolated between the value for inverted annular and dispersed flow film boiling.


## 3. CTF Application to the Steady-State Critical Power Exercise of the OECD/NRC BWR Full-Size-Fine-Mesh Bundle Test Benchmark

A full-scale bundle, simulating an 8 × 8 high burn-up fuel assembly, was installed in the NUPEC BWR test section for pressure drop, void distribution, and critical power measurements [1]. Three combinations of radial and axial power shapes were tested: (1) beginning of cycle (BOC) radial power pattern/cosine axial power shape; (2) end of cycle (EOC) radial power pattern/cosine axial power shape; and (3) beginning of cycle radial power pattern/inlet peaked axial power shape. The individual radial and axial power distributions for all three combinations are provided in Volume I of the BFBT benchmark specifications [1]. The steady-state test series consisted of three parts: void distribution tests, pressure drop tests, and critical power tests. The pressure drop was measured in both single-phase flow and two-phase flow conditions that cover the normal operational behavior.CTF has been previously applied to the steady-state and transient void distribution exercises and to the single- and two-phase pressure drop exercises of the BFBT benchmark [4]. These studies have indicated that the code reproduces the qualitative behavior of the steady-state and transient void fraction distributions, but quantitatively overpredicts the vapor content in the BFBT bundles. This coincided with the results of the two-phase pressure drop exercises where the total pressure drop was also slightly overpredicted. The cause of both phenomena is believed to be an overestimation of the interfacial drag forces leading to an overpredicted slip and, subsequently, to an underpredicted vapor velocity yielding a higher void fraction. On the other hand, the code was able to reproduce the radial void distribution, except for the regions next to unheated structures. Regarding the sensitivity to the turbulent mixing coefficient, the code performed better when smaller values were used; larger mixing coefficients resulted in less accurate in-bundle void distribution [4].In the BFBT tests, the critical power was measured by slowly increasing the bundle power while monitoring the individual heater rod thermocouple signals. The critical power was defined when the peak rod surface temperature became 14°C higher than the steady-state temperature level before dryout occurred. The dryout was observed in the peak power rod located at the peripheral row adjacent to the channel box. The boiling transition was always observed just upstream of spacers. The estimated accuracies of the major process parameters were 1% for the pressure and 1.5% for the power. Figure1 describes the definition of thermocouple position. Each thermocouple position was identified according to rod no., axial location, and rotational angle.Figure 1
Definition of thermocouple position for assembly C2A [1].In this paper, only results for assembly C2A are shown, because C2A is the bundle used in the above-mentioned assessment of the CTF turbulent mixing model. The supplied measured data includes critical power, axial location of boiling transition, and corresponding boundary conditions (pressure, flow, inlet subcooling, and power shapes). The radial and axial power profiles of assembly type C2A are given in Table1.Table 1
Steady-state critical power measurement conditions for assembly C2A [1].
Axial power profile
Radial power profile
Boundary conditions
Pressure (MPa)
5.5, 7.2, 8.6
Flow rate (t/h)
10, 20, 30, 45, 55, 60, 65
Inlet subcooling (KJ/kg)
25, 50, 84, 104, 126
No. of cases
Exercise cases
14A full C2A bundle model on a subchannel by subchannel resolution (no symmetry) was used in the calculations. The heated length was divided axially into forty (40) equidistant nodes. The number of axial nodes was selected based on the desired aspect ratio and the modelling of the spacer grids. Previous verification studies have shown that the number of axial nodes will not influence the stability of the CTF thermal-hydraulic solution, but will impact the magnitude of the local pressure drop at the spacer locations.The pressure losses due to spacer grids were calculated as velocity head losses with subchannel loss coefficient as calculated by the Shiralkar’s method [4]. The total cross-flow between two adjacent subchannels is due to diversion cross-flow by lateral pressure gradients and cross-flow by turbulent mixing and void drift. Turbulent mixing and void drift phenomena are modeled in CTF with the Lahey and Moody approach [29], where the net two-phase mixing (including void drift) is assumed to be proportional to the nonequilibrium void fraction gradient. The void drift is only assumed to occur in bubbly, slug, and churn flow, where liquid is the continuous phase and vapor is the dispersed phase. The single-phase mixing coefficient is either specified as an input value or calculated using an empirical correlation derived by Rogers and Rosehart [30]. The Beus model for two-phase turbulent mixing is utilized [31].It is important to mention that the code version used in this study does not consider a multifilm formulation within a subchannel.Following the experimental procedure, in the code calculations the power was increased gradually until the rod surface temperature became 14°C higher than the temperature at the previous steady-state level; then a dryout occurrence was considered and the critical power was determined. It was found that the criterion of 14°C temperature excursion agrees perfectly with the code prediction of complete film dryout.Sensitivity studies were performed on the effect of turbulent mixing modeling. The calculations were repeated using three different options: (1) without modelling of turbulent mixing and void drift; (2) Lahey and Moody model with a user-specified single-phase mixing coefficient; and (3) Lahey and Moody model with a single-phase mixing coefficient by Rogers and Rosehart’s correlation and Beus’ model for two-phase mixing. Results are summarized in Figures2, 3, 4, 5, and 6 in the form of predicted versus measured (P/M) values. As it can be seen, the best agreement is obtained when turbulent mixing and void drift were not modeled. For this case, the mean relative error in the code predictions was found to be 3.4%. These results should not be misunderstood—it is a fundamental fact that the lateral exchange of momentum, mass, and heat due to increased turbulence in the flow improve the heat transfer rates. However, for this particular bundle setup, both models available in CTF had an adverse effect on the code accuracy of dryout prediction resulting in an overestimation of the critical power. Stronger was the turbulent mixing, larger was the overprediction (for typical BWR bundles, Rogers and Rosehart’s correlation generally gives a single-phase mixing coefficient in the order of 10E-3). This is in an agreement with the findings in [4]—stronger mixing, less accurate in-bundle void distribution: larger mixing coefficients led to more homogeneous void distribution pushing the vapor phase away from the peripheral region toward the central channels and resulting in a significant overprediction of the void fraction in the channels next to the water rod. In Table 2, results for a representative test case, SA505500, are given as an example. It can be clearly seen that the current turbulent mixing model gives inaccurate, and in this case less conservative, estimates of the critical power. It should be kept in mind that the axial location of the dryout could be everywhere between thermocouples A and B in Figure 1.Table 2
Results for test case SA505500.
SA505500
Measured data
Without turbulent mixing and void drift
Lahey and Moody model [29]
Two-phase mixing enhancement by Beus
Mixing coefficient of 0.03
Mixing coefficient byRogers and Rosehart [30]
Critical power
6.13 MW
6.67 MW
8.25 MW
7.59 MW
Radial location
Peripheral rod
Peripheral rod
Central rod next to water channel
Central rod next to water channel
Axial location
3.521 m
3.335 m
3.612 m
3.612 mFigure 2
Predicted versus measured (P/M) critical power for assembly C2A with different turbulent mixing models.Figure 3
P/M Critical Power versus Power for Assembly Type C2A.Figure 4
P/M critical power versus flow rate for assembly type C2A.Figure 5
P/M critical power versus inlet sub-cooling for assembly type C2A.Figure 6
P/M critical power versus pressure for assembly type C2A.On the other hand, the spacers instrumented along the C2A bundle are ferrule type spacers which are not designed to enhance the flow turbulence and, therefore, any attempt to model enhanced turbulent mixing at the spacer locations (by large coefficients) would be meaningful. Also, in the CTF simulations, the lateral pressure gradient due to spacers was accounted for by applying subchannel-based loss coefficients in both axial and transverse directions. Cross-flow due to coolant temperature and density gradients was handled by the diversion cross-flow models.Another observation in the code predictions was the bias with the pressure—the code tends to overpredict the critical power at lower pressure (~5.5 MPa) and to underpredict it at higher pressure (~9 MPa). No bias was seen with the flow rate and the inlet subcooling.


## 4. CTF Application to the Steady-State DNB Exercise of the OECD/NRC PWR Subchannel and Bundle Tests Benchmark

In the NUPEC PWR DNB measurements, the test assembly configuration consisted of twenty-five (25) rods in a 5 × 5 square bundle or thirty-six rods (36) in a 6 × 6 square bundle [2]. The configuration of rods in this geometry approximated a typical 17 × 17 commercial power reactor fuel assembly. Each rod had a heated length of 3.658 m, an outer diameter of 9.5 mm, and a rod pitch of 12.6 mm. Between thirteen (13) and seventeen (17) spacers (both with and without mixing vanes) along the axial length supported the rods in a vertical grid. The rods were cylindrical in shape with a hollow insulator of alumina radially encircled in a heater made from Inconel 600. For the steady-state departure from nucleate boiling cases considered in this paper, a series of experiments were performed in five different configurations [2]. The NUPEC test series (numbered 0, 2, 3, 4, 8, and 13) were conducted at various pressures and temperatures where prior experience demonstrated that departure from nucleate boiling was likely to occur. The thermocouples were attached to the inner surface of the heater rods to determine the boiling transition. The bundle power was increased gradually by fine steps to the vicinity of DNB power, which was based on preliminary analysis. The occurrence of DNB was confirmed by a rod temperature rise of more than 11°C as measured by the thermocouples. The DNB power was defined as the power corresponding to the step just before the step where the temperature increased. Figure 7 shows the axial position of the thermocouples for each configuration. The various test configurations used several axial and radial power schemes, which provided an ample cross-section of calculation data. Five assemblies were utilized, denoted by A0, A2, A3, A4, and A8 (see Tables 3, 4, and 5). The estimated accuracies of different process parameters for the DNB measurements were 1% for the pressure; 1.5% for the flow; 1°C for the fluid temperature; and 1% for the power.Table 3
Geometry and power shape for test assembly A0 [2].
Item
Data
Assembly
Rods array
5 × 5
Number of heated rods
25
Number of thimble rods
0
Heated rod outer diameter (mm)
9.50
Thimble rod outer diameter (mm)
—
Heated rods pitch (mm)
12.60
Axial heated length (mm)
3658
Flow channel inner width (mm)
64.9
Radial power shape
A
Axial power shape
Uniform
Number of MV spacers
5
Number of NMV spacers
2
Number of simple spacers
6
MV spacer location (mm)
610, 1219, 1829, 2438, 3048
NMV spacer location (mm)
0, 3658
Simple spacer location (mm)
305, 914, 1524, 2134, 2743, 3353Table 4
Geometry and power shape for test assemblies A2 and A3 [2].
Item
Data
Assembly
Rods array
5 × 5
6 × 6
Number of heated rods
25
36
Number of thimble rods
0
0
Heated rod outer diameter (mm)
9.50
9.50
Thimble rod outer diameter (mm)
—
—
Heated rods pitch (mm)
12.60
12.60
Axial heated length (mm)
3658
3658
Flow channel inner width (mm)
64.9
77.5
Radial power shape
A
D
Axial power shape
Uniform
Uniform
Number of MV spacers
7
7
Number of NMV spacer
2
2
Number of simple spacers
8
8
MV spacer location (mm)
457, 914, 1372, 1829, 2286, 2743, 3200
NMV spacer location (mm)
0, 3658
Simple spacer location (mm)
229, 686, 1143, 1600, 2057, 2515, 2972, 3429Table 5
Geometry and power shape for test assemblies A4 and A8 [2].
Item
Data
Assembly
Rods array
5 × 5
5 × 5
Number of heated rods
25
24
Number of thimble rods
0
1
Heated rod outer diameter (mm)
9.50
9.50
Thimble rod outer diameter (mm)
—
12.24
Heated rods pitch (mm)
12.60
12.60
Axial heated length (mm)
3658
3658
Flow channel inner width (mm)
64.9
64.9
Radial power shape
A
B
Axial power shape
Cosine
Cosine
Number of MV spacers
7
7
Number of NMV spacer
2
2
Number of simple spacers
8
8
MV spacer location (mm)
471, 925, 1378, 1832, 2285, 2739, 3247
NMV spacer location (mm)
2.5, 3755
Simple spacer location (mm)
237, 698, 1151, 1605, 2059, 2512, 2993, 3501Figure 7
Axial thermocouple locations in the PSBT DNB measurements [2].Approximately, twenty-five (25) tests were chosen from configurations A4 (test series 4 and 13) and A8 (test series 8) as good candidates for the benchmark. Additionally, ten (10) tests were chosen from each of the remaining test series 0, 2, and 3. In total, one hundred (100) tests were modeled with CTF.Full-bundle models on a subchannel-by-subchannel scale were used in the calculations. The heated length was divided axially into seventy (70) equidistant nodes. The number of axial nodes was selected based on the modelling of the spacer grids. The spacers’ pressure losses were calculated with subchannel loss coefficients as specified in [2].The PSBT bundles were equipped with three (3) different spacer types: single support spacer, nonmixing vane spacers, and mixing vane spacers [2]. While the first two types of spacers mostly affect the pressure drop in the bundle, the third type increases the turbulence and creates strong cross-flows between the subchannels. The code version utilized in this work does not have the capabilities to simulate transverse flows locally created by the mixing vanes. Instead, a large overall single-phase mixing coefficient of 0.05 was used. The choice of the mixing coefficient value was based on previously performed sensitivity studies on void distribution tests [6].The default models, described in Section2, for the flow and heat transfer regimes were utilized in these CTF simulations. The calculations were performed in two sets: (1) without modelling of turbulent mixing and void drift; and (2) turbulent mixing and void drift by Lahey and Moody with a user-specified single-phase mixing coefficient of 0.05 and Beus’ model for two-phase mixing enhancement. The simulations without turbulent mixing and void drift were carried out to confirm or disclaim the observations from the C2A BFBT bundle tests, where the model introduces high inaccuracy in the code predictions of void distribution and critical power.Following the experimental procedure, in the simulations the power was increased gradually until the rod surface temperature became 11°C higher than the temperature at the previous steady-state level; then a dryout occurrence was considered and the critical power was determined.Code-to-data comparisons are given in Figures8 and 9. Unlike the C2A BFBT bundle results, the agreement is significantly improved when turbulent mixing and void drift are modelled and a large underprediction on the DNB power is seen if not. No bias with the power, the flow, rate and the subcooling was found (Figures 10, 11, and 12). Similarly to the C2A BFBT critical power calculations, a code bias with the pressure was seen (Figure 13)—the code tends to overpredict the critical power at lower pressure (~5 MPa) and to underpredict it at higher pressure (~15 MPa).Figure 8
Predicted versus Measured DNB without Turbulent Mixing and Void Drift Modeling.Figure 9
Predicted versus measured DNB with turbulent mixing and void drift modeling (single-phase mixing coefficient of 0.05 and two-phase mixing by Beus).Figure 10
P/M DNB power versus power for assembly types A0, A2, A3, A4, and A8.Figure 11
P/M DNB power versus flow rate for assembly types A0, A2, A3, A4, and A8.Figure 12
P/M DNB power versus subcooling for assembly types A0, A2, A3, A4, and A8.Figure 13
P/M DNB power versus pressure for assembly types A0, A2, A3, A4, and A8.


## 5. Statistical Analysis

Mean error and standard deviation are calculated for each test series. The mean error is represented asCP¯=∑n=1NCPn/N, where the critical power (CP) error for test case “n” is given as CPn=(CPcoden-CPexpn) and N is the total number of test cases.The standard deviation is given asσ=±∑n=1N(CPn-CP¯)2/N-1.Results are given in Tables6 and 7, respectively for the BFBT assembly C2A and the PSBT assemblies A0, A2, A3, A4, and A8. As previously discussed, the sensitivity to the turbulent mixing and void drift models showed very inconsistent behavior—the modeling of stronger turbulent mixing worsened the code predictions for critical power in the BWR-type bundle (C2A) and vice-versa improved the code predictions for the departure from nucleate boiling in the PWR-type bundles (A0, A2, A3, A4, and A8).Table 6
Statistical analysis for the BFBT test assembly C2A.
Bundle type
Turbulent mixing:β = 0.03
Turbulent mixing: Rogers and Rosehart [30]
Without turbulent mixing
Mean error
Standard deviation
Mean error
Standard deviation
Mean error
Standard deviation
C2A
−1.63
0.24
−1.03
0.23
0.09
0.32Table 7
Statistical analysis for the PSBT test assemblies A0, A2, A3, A4, and A8.
Bundle type
Turbulent mixing:β = 0.05 and Beus
Without turbulent mixing
Mean error
Standard deviation
Mean error
Standard deviation
A0
0.04
0.18
−0.39
0.12
A2
0.34
0.25
−0.91
0.55
A3
−0.10
0.40
−1.60
0.66
A4(TS4)
−0.09
0.31
−0.63
0.41
A8
0.09
0.29
−0.35
0.24
A4(TS13)
−0.08
0.25
−0.67
0.37
Overall
0.01
0.30
−0.62
0.47To recall, no empirical correlations were used to calculate the critical power/critical heat flux. In annular flow regime, the heat flux was not limited by a correlation, but rather forced convection vaporization existed until the film dries out. In CTF, film dryout is a complex function of the film flow rate, the applied heat flux, and the entrainment/de-entrainment rate, and is determined by the solution of the hydrodynamic equations. In these assessments, a rapid increase of the cladding temperature was used as a criterion for DNB or dryout occurrence.


## 6. Conclusions

To assess its accuracy of dryout and departure from nucleate boiling prediction, the subchannel thermal-hydraulic code CTF was applied to Exercise II-1 (steady-state critical power) of the OECD/NRC BFBT benchmark and Exercise II-1 (steady-state departure from nucleate boiling) of the OECD/NRC PSBT benchmark. The obtained results showed that the code predicts fairly well the critical power and departure from nucleate boiling power with no specific tendency of over- or underprediction. The boiling crisis location was well captured. However, some areas of potential improvement were identified.The turbulent mixing and void drift model, including spacer grid effects, has to be enhanced. The current approach of using a larger overall mixing coefficient to capture spacer induced turbulence had proven to be highly inaccurate.The bias with the pressure indicated possible inconsistencies in the fluid solution and a need of further improvements of flow regime transition logic in CTF.The code tendency to overpredict the void generation rates and two-phase pressure drop may contribute to an inaccurate boiling crisis prediction and, therefore, the interfacial friction models might need an improvement.In summary, this work indicated that some code models have to be further improved to address the new trends in nuclear reactor core designs.



---

*Source: 508485-2013-02-03.xml*