About this Journal Submit a Manuscript Table of Contents
Science and Technology of Nuclear Installations
Volume 2012 (2012), Article ID 465059, 12 pages
Research Article

Validation of the Subchannel Code SUBCHANFLOW Using the NUPEC PWR Tests (PSBT)

Institute for Neutron Physics and Reactor Technology (INR), Karlsruhe Institute of Technology (KIT), Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany

Received 19 March 2012; Accepted 5 June 2012

Academic Editor: David Novog

Copyright © 2012 Uwe Imke and Victor Hugo Sanchez. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


SUBCHANFLOW is a computer code to analyze thermal-hydraulic phenomena in the core of pressurized water reactors, boiling water reactors, and innovative reactors operated with gas or liquid metal as coolant. As part of the ongoing assessment efforts, the code has been validated by using experimental data from the NUPEC PWR Subchannel and Bundle Tests (PSBT). The database includes single-phase flow bundle outlet temperature distributions, steady state and transient void distributions and critical power measurements. The performed validation work has demonstrated that the two-phase flow empirical knowledge base implemented in SUBCHANFLOW is appropriate to describe key mechanisms of the experimental investigations with acceptable accuracy.

1. Introduction

The requirements for computational resources of high-resolution computational fluid dynamics (CFD) are still large, so the thermal-hydraulic analysis of nuclear reactor cores is frequently performed using subchannel computer codes. The current development of subchannel codes concentrates on refined modeling of two-phase flow. The two-fluid formulation separates the conservation equations of mass, energy, and momentum to vapor and liquid. The COBRA-TF family of subchannel codes [1, 2] extends this treatment to a different description of continuous liquid and entrained liquid droplets, which results in a set of nine conservation equations. This kind of refinement needs additional constitutive relations which have to be derived from single-effect experiments. For example, the behavior of droplets colliding with grid spacers and the entrainment of droplets in the vapor flow has to be described by physical models, which are validated using such experimental data [3]. In addition, the computational requirements are growing strongly along with the number of equations to be solved.

For the design and safety assessment of nuclear power plants, the coupled multiphysics description of the core behavior becomes more and more important [4]. Fast running and extensively validated numerical tools are needed for industrial, regulatory and research purposes. Therefore a good performance of the thermal-hydraulic calculation is essential, if it is combined with numerical simulation of neutron physics or fuel pin mechanics.

An alternative to the simulation of the processes on a microscale level is to use empirical correlations related to pressure drop, heat transfer, void generation, and so forth collected over the last decades. These correlations are combined with liquid-vapor mixture equations for the conservation of mass, momentum, and energy as used by legacy codes. After a critical review of various legacy subchannel codes freely available, it was decided to develop SUBCHANFLOW as a modern and modular subchannel code starting from the COBRA-family [5, 17, 18] using the above-mentioned three-equation approach for the mixture of liquid and vapor. The experimental data from the NUPEC PWR Subchannel and Bundle Tests (PSBT) [16] are used to show the quality of results which can be reached by a simple approach based on validated correlations.

The present paper starts with an overview about the main features of the new code. Subsequently, the validation procedure for PWR conditions is described starting from single-channel simulations. The mixing parameter to be used in the PWR-related investigations is deduced from steady-state single-phase flow bundle experiments. The boiling models are validated against void fraction distribution measurements for bundle steady-state and transient flow tests. Finally, the critical heat flux phenomenon is investigated and the SUBCHANFLOW predictions are compared to measurement data. All results are obtained with the current version SUBCHANFLOW 2.1.

2. Main Features of SUBCHANFLOW

2.1. Programming Features

SUBCHANFLOW is a fast running and flexible simulation tool, which is easy to maintain by keeping the code structures as simple as possible. The bases for the source code are the legacy subchannel programs COBRA-IV-I [18] and COBRA-EN [5]. The old methods of data management like Fortran EQUIVALENCE and swapping to hard disk are removed. The Fortran COMMON structure is replaced by a global data structure centralized in one single Fortran module complemented by a description of each global variable name. The thermophysical properties of the coolants and solid materials are summarized in separate modules. All arrays are dynamically allocated depending on the problem-specific input data. The code is prepared to be used as a library called by another simulation tool. There is an error management control that will shut down the code in case of difficulties closing all files and deallocating the memory. The portability of the code is assured by avoiding functions that depend on operating systems. Consequently, SUBCHANFLOW can be compiled under WINDOWS, LINUX, or other UNIX systems using a standard Fortran 95 compiler. The input deck is designed as a text-based “User Interface” with comprehensive keywords and simple tables. Long tables can be fed in by external files named in the input deck. A manifold output is created to be used with different postprocessing tools for example to generate simple curves or more extensive three-dimensional diagrams.

2.2. Modeling Features

SUBCHANFLOW can handle both rectangular and hexagonal fuel bundles and core geometries built from these. As boundary conditions, the total flow rate or a channel-dependent flow rate can be selected. It is possible to distribute the flow automatically to the parallel channels depending on the friction at the bundle inlet. In addition, a pure top-bottom pressure difference boundary can be applied for steady-state calculations. Fluid temperature at the inlet and pressure at the outlet always have to be prescribed as boundary conditions [19].

In opposite to the majority of subchannel codes, SUBCHANFLOW uses rigorously SI units internally in all modules [20]. Modern coolant properties and state functions are implemented for water using the IAPWS-97 formulation (The International Association for the Properties of Water and Steam). In addition, property functions for liquid metals (sodium and lead) and gases (helium and air) are available too. An iterative steady-state numerical procedure is available to determine the power at which critical heat flux conditions appear during the simulation.

In SUBCHANFLOW, profit is taken from the many valuable empirical correlations for pressure drop, heat transfer coefficients, void generation, etc., collected over the last decades. Consequently, it does not follow the general trend to describe two-phase flow by simulating the processes on a microscale basis (e.g., separate conservation equations for liquid droplets, films or vapor bubbles). In SUBCHANFLOW, a three-equation two-phase flow model that is a mixture equation for mass, momentum, and energy balance is implemented. The constitutive relations are expressed as mixture equations for wall friction and wall heat flux as well as a slip velocity relation. In addition, user defined empirical correlations can be implemented. In the present paper, only the correlations used for the actual PSBT validation campaign are mentioned.

2.3. Basic Conservation Equations

In subchannel codes, a channel consists of a finite fraction of the total cross-sectional area of the nuclear reactor core region. The smallest possible channel would be the size of a subchannel surrounded by the fuel rods (see Figures in Table 3). For the numerical simulation, a subchannel is divided into several axial mesh volumes. Transport of mass, momentum, and energy is possible along the axial direction and between the neighboring channels through the gap formed by the fuel rods (lateral direction, cross-flow). The basic transport equations are based on the Euler approach including friction at solid surfaces. In the lateral momentum equation, the convective transport of lateral momentum is neglected because the friction term dominates the cross-flow. In the energy equation, the heat flux from the rod surfaces is the main source term. For transient conditions, a slip between vapor and liquid is taken into account in the enthalpy time derivative. Turbulent transport of momentum and energy between neighboring channels is described by a simple empirical mixing model. For an axial mesh volume (𝑗) in channel (𝑖) surrounded by volumes of neighboring channels (𝑛) through a gap (𝑘), the basic conservation equations in finite difference form are the following.

Mass conservation: 𝐴𝑖,𝑗Δ𝑋𝑗𝜌Δ𝑡𝑖,𝑗𝜌old𝑖,𝑗+𝑚𝑖,𝑗𝑚𝑖,𝑗1+Δ𝑋𝑗𝑘𝑤𝑘,𝑗=0.(1)

Energy conservation: 𝐴𝑖,𝑗Δ𝑡𝜌𝑖,𝑗𝑖,𝑗old𝑖,𝑗+𝑖,𝑗𝜌𝑖,𝑗𝜌old𝑖,𝑗+1Δ𝑋𝑗𝑚𝑖,𝑗𝑖,𝑗𝑚𝑖,𝑗1𝑖,𝑗1+𝑘𝑤𝑘,𝑗𝑘,𝑗=𝑄𝑖,𝑗𝑘𝑤𝑘,𝑗𝑖,𝑗𝑛(𝑘),𝑗.(2)

Axial momentum: Δ𝑋𝑗𝑚Δ𝑡𝑖,𝑗𝑚old𝑖,𝑗+𝑚𝑖,𝑗𝑈𝑖,𝑗+Δ𝑋𝑗𝑘𝑤𝑘,𝑗𝑈𝑘,𝑗=𝐴𝑖,𝑗𝑝𝑖,𝑗𝑝𝑖,𝑗1𝑔𝐴𝑖,𝑗Δ𝑋𝑗𝜌𝑖,𝑗12Δ𝑋𝑓Φ2𝐷𝜌liq+𝐾𝑣𝑖,𝑗||𝑚𝑖,𝑗||𝑚𝑖,𝑗𝐴𝑖,𝑗Δ𝑋𝑗𝑘𝑤𝑘,𝑗𝑈𝑖,𝑗𝑈𝑛(𝑘),𝑗.(3)

Lateral momentum: Δ𝑋𝑗𝑤Δ𝑡𝑘,𝑗𝑤old𝑘,𝑗+𝑈𝑘,𝑗𝑤𝑘,𝑗𝑈𝑘,𝑗1𝑤𝑘,𝑗1=𝑠𝑘𝑙𝑘Δ𝑋𝑗Δ𝑝𝑘,𝑗1𝐾𝐺Δ𝑋𝑣𝑘𝑠𝑘𝑙𝑘𝑗||𝑤𝑘,𝑗||𝑤𝑘,𝑗.(4)

Slip correction for enthalpy derivative based on the two-phase flow basic energy equation (Tong’s function) [21]: 𝜌=𝜌oldfg𝜕𝜓𝜕𝜓=𝜌liq𝑥(1𝛼)𝜌vap𝛼(1𝑥).(5)

Effective specific volume 𝑥𝑣=2𝛼𝜌vap+(1𝑥)2(1𝛼)𝜌liq𝑈=𝑚𝐴𝑣.(6)

2.4. Numerical Solution Procedure

The conservation equations along with the constitutive equations, for example, to calculate the void fraction, represent the system of equations of the mixture two-phase flow model. The basic flow variables such as the axial and lateral flow rates, the pressure, the enthalpies, and the void fractions are calculated in each time step axially layer by layer. For each axial layer, the coolant enthalpies are calculated first from the energy conservation equations. The axial pressure gradients are calculated from the combined axial and transverse momentum equation. The mass flow rates in the lateral directions are calculated from the transverse momentum equation, knowing the axial pressure gradients. The axial mass flow rates are deduced from the mass continuity equation. From the enthalpy of each computational cell, the steam quality and then, through the quality/void correlation, the steam volume fraction and hence the coolant density are computed. This procedure is repeated several times during each time step resulting in a fully implicit scheme. For steady-state calculations, the time step is set to a very large value. The sketched solution algorithm is limited to cases with axial flow rates which always keep positive (upflow). The linear equation system built up by the energy equations in each layer is solved by the SOR (successive overrelaxation) method. The equation system for the pressure gradients can be solved by a direct scheme or by the SOR method which needs much less computer storage. The fuel rod or heater rod temperatures are calculated in each iteration step depending on power release and cladding to coolant heat transfer. For each axial layer, the rod is divided into a number of radial rings to solve the heat conduction equation in radial direction by a finite volume method. Axial heat conduction can be accounted for in transient simulations, if necessary.

3. Validation Using PSBT Benchmark Data

The validation of SUBCHANFLOW started with earlier versions [19, 20, 22] and it was continuously repeated with each new major release. Here, the validation work performed with the latest version will be presented and discussed. A detailed description of the PSBT benchmark as well as of the experimental data is given in [16].

3.1. Flow in a Single Subchannel

A typical NUPEC test assembly consists of four different geometrical types of subchannels: the center channels surrounded by four rods, the center channels surrounded by 3 fuel rods and a guide tube, the side channels surrounded by two rods and a part of the assembly wall and the corner channel surrounded by one rod and two wall parts (see Table 3). A basic test for a subchannel code is to simulate boiling in these kinds of heated single channels. Then, the predicted void fraction can be compared to the one measured in the PSBT tests. The single subchannel test section is uniformly heated along 1555 mm and the void measurement was made at 1400 mm from the bottom of the heated section. Tests are done with different pressure, coolant flow rate, power, and inlet temperature as given in Table 1.

Table 1: Boundary conditions for single channel tests.

The closure correlations used in SUBCHANFLOW are summarized in Table 2 while in Table 3 the geometrical configuration and the global result expressed by the standard deviation (see Equation (7)) between experimental values and simulation results are shown: 𝜎=1𝑁1𝑁𝑖=1𝑇2exp,𝑖𝑇2sim,𝑖.(7)

Table 2: Closure correlations.
Table 3: Global result for different subchannel types.

In Figures 1, 2, 3, and 4 the comparison of the measured and calculated void fraction in the different types of single subchannels is exhibited. The ±0.05 absolute void fraction envelope of the exact agreement is given by the pink straight lines. The standard deviation for the difference of measured and calculated values is given in Table 3. The estimated accuracy for the void measurements is 0.03 (absolute void fraction) [16].

Figure 1: Comparison between calculated and measured void fraction for subchannel type S1.
Figure 2: Comparison between calculated and measured void fraction for subchannel type S2.
Figure 3: Comparison between calculated and measured void fraction for subchannel type S3.
Figure 4: Comparison between calculated and measured void fraction for subchannel type S4.

The majority of the points lie within the 0.05 void fraction error band. Only for the case S3 SUBCHANFLOW considerably overpredicts the void fraction as long as it is below 0.2. There is also a large under prediction of one single point for a void fraction of around 0.5 (see Figure 3). There is no clear and systematic explanation for this behavior. Outliers in the measurement procedure are supposed.

3.2. Determination of Cross-Flow Mixing Coefficient

Under single-phase flow conditions, the cross-flow between neighboring subchannels is divided into two categories: turbulent mixing and diversion cross-flow. The turbulent mixing is an inter-subchannel mixing due to turbulence of the fluid flow. By this mixing, momentum and energy transfer between subchannels take place, but no net mass transfer. A void drift model leading to turbulent exchange of vapor volume between neighboring channels is implemented in SUBCHANFLOW in addition. In the present benchmark, it does not give better results than the standard energy transfer model and is not used. The diversion cross-flow occurs due to lateral pressure gradients, which may be introduced by differences of subchannel geometry or obstructions such as spacers.

In most subchannel codes, the turbulent cross-flow between channel 𝑖 and channel 𝑗 through gap 𝑘is defined by 𝑤𝑖,𝑗=𝛽𝑖,𝑗𝑠𝑘𝐺0.5𝑖+𝐺𝑗,𝑤𝑖,𝑗=𝑤𝑗,𝑖.(8)

The subchannel analysis of bundles is very sensitive to the mixing coefficient 𝛽. Several correlations are published to model the influence of Reynolds number and channel geometry [23, 24]. In case of assemblies containing mixing vane spacers, a simple method is to use a constant value, which may be determined by single-phase flow experiments measuring the exit subchannel temperature [13]. There is a steady-state fluid temperature experiment in the PSBT test series that uses a special lateral power distribution for a 5 by 5 rod array containing 7 mixing vane spacers (see Table 4). These data are used to find a mixing coefficient valid for the assembly under consideration. Figure 5 shows the comparison between calculated and measured temperatures at the outlet of each subchannel for several typical flow conditions. An optimum mixing coefficient of 0.06 was found with a standard deviation of about 5°C.

Table 4: Steady-state mixing experiment conditions.
Figure 5: Comparison of calculated and measured channel outlet temperatures for a mixing coefficient of 0.06 (envelope: ±10°C).
3.3. Steady-State Bundle Void Fraction Distribution

Table 5 gives an overview about all test bundle geometries used in the present paper for void fraction and DNB (departure from nucleate boiling) simulations.

Table 5: Geometry of test bundles for void and DNB measurements.

The main measurement data for bundle experiments provided by PSBT are X-ray densitometer measurements of void fraction at three axial elevations (2.2 m, 2.7 m, and 3.2 m related to a heated length of 3.7 m). The resulting void is an average over the four central subchannels of the bundle. Four steady-state bundle tests are investigated (referenced here as cases B5, B6, B7, B8). The basic rod configuration is the same as for the steady-state single-phase outlet temperature test. B5 has a uniform axial power profile. B6 and B7 have a cosine profile. In the case of B7, the central rod is replaced by a guide tube with a diameter of 12.24 mm. The radial power profile is described by 9 central rods having a relative power of 1.0, whereas the boundary rods are heated with a relative power of 0.85. B8 is a repetition case of B5. In Table 6, the range of boundary conditions is given for all tests performed to measure steady-state bundle void fraction distributions.

Table 6: Boundary conditions for steady state bundle void fraction tests.

In addition to the correlations used for the simulation of the single-channel experiments, SUBCHANFLOW uses the closure laws summarized in Table 7.

Table 7: Empirical correlations used for the bundle tests.

In Figures 6, 7, and 8 the comparison of the calculations with the measured void fractions at the three axial levels is shown for case B5. Again the 0.05 void envelope is indicated. Generally, the majority of the SUBCHANFLOW predictions are in acceptable agreement with the data. But there is a tendency of SUBCHANFLOW to overpredict the void fraction at the lower and middle parts of the bundle for void fractions below 0.4, while SUBCHANFLOW tends to underpredict the void fraction at the upper level of the bundle when the void fraction is between 0.3 and 0.55.

Figure 6: Comparison of calculated and measured void fraction for case B5 at the lower axial level.
Figure 7: Comparison of calculated and measured void fraction for case B5 at the middle axial level.
Figure 8: Comparison of calculated and measured void fraction for case B5 at the upper axial level.

Table 8 summarizes the results of all 4 test cases. The accuracy of the experimental data is about 0.04 (absolute void fraction) [16]. Depending on case and level, the standard deviation is varying between 0.04 and 0.08.

Table 8: Void fraction standard deviations for all steady state test cases.
3.4. Transient Bundle Void Fraction Distribution

Transient void tests were performed with three different bundle types (B5, B6 and B7) representing PWR-relevant scenarios: power increase (PI), flow reduction (FR), depressurization (DP), and inlet temperature increase (TI). In the present paper, only the test cases for bundle type B7 are presented and discussed. All four scenarios have nearly the same initial conditions as documented in Table 9. The transient boundary conditions can be found in Figures 912. All transients lead to an increase of void fraction at the bundle exit. The void measurement technique was the same as that used for the steady-state tests and the void fraction data were taken again at the three axial elevations (lower, middle, upper). They are again averaged over the four central subchannels of the 5 × 5 test bundle.

Table 9: Initial conditions for B7 transients.
Figure 9: Transient power during the PI test.
Figure 10: Transient flow rate behavior in the FR test.
Figure 11: Pressure decrease in the DP test.
Figure 12: Inlet temperature and pressure during the TI test.

The calculated time-dependent void profiles are very similar to the measured values as can be observed in Figures 13, 14, 15, and 16. Based on the comparison of predictions and data, it can be stated that the void fractions predicted at the lower axial position are overestimated. On the other hand, the comparison of the predicted void fraction at the upper level is slightly underpredicted. In general, most of the predictions are qualitatively following the evolution of the measured data during the transient test phase.

Figure 13: Void fraction comparison for power increase transient.
Figure 14: Void fraction comparison for flow reduction transient.
Figure 15: Void fraction comparison for power depressurization transient.
Figure 16: Void fraction comparison for inlet temperature increase transient.

The considerable overprediction of the void fraction by SUBCHANFLOW at the lower bundle part for the overpower (Figure 13) and mass flow reduction (Figure 14) transients indicates that the subcooled boiling model may need further review. The comparison of the SUBCHANFLOW predictions with the measured void for the temperature transient is quite reasonable for all three levels (see Figure 16).

3.5. Steady-State Critical Heat Flux

There are several hundreds of critical heat flux (CHF) correlations for tubes, annular space, and other simple geometries in the open literature, which can be used to calculate the occurrence of critical boiling conditions. For fuel assemblies, a more accurate method is to use a specific CHF correlation or lookup tables especially fitted for the specific geometry and spacer grid type. The geometry and the grid spacers have a considerable influence on the CHF phenomenon. For the same thermohydraulic local conditions, the CHF may strongly vary. The well-known EPRI correlation [25] is a generalized bundle correlation that is used in COBRA-EN [5]. It covers PWR and BWR normal operating conditions as well as loss of coolant accident boundary conditions. This correlation is implemented in SUBCHANFLOW together with an iterative method to determine steady-state critical heat flux without simulating a transient. It is used in the present subchannel mode simulations instead of simple-geometry-based correlations with correction factors because typical bundle properties are implicitly included. The DNB benchmark data provide the power at which the critical heat flux condition is met for various bundle geometries (Table 5). In addition, the axial power profile and the radial power distribution changed in the different experiments. The radial power distribution has the same scheme for all DNB tests. The rods at the edge of the assembly have a relative power of 0.85, as all inner rods have a relative power of 1.0.

The occurrence of DNB during the tests is detected by a rod temperature rise of 11°C measured by the thermocouples. The corresponding power is defined as critical power. The temperature is measured with an accuracy of 1°C, whereas the power has a measurement error of 1% [16].

The SUBCHANFLOW predictions for many tests are compared to the measured critical power in Figure 17 to Figure 22. An error band of ±10% is indicated in the graphs.

Figure 17: Measured and simulated critical power values for assembly A0.

The predictions of the critical power for the tests of bundle type A0 have a standard deviation of 8 % but with a tendency to an overprediction (Figure 17). On the contrary, the critical power predictions for the bundle type A2 are generally underpredicted and several predictions are outside the error band of ±10%. The only difference between the bundle types A0 and A2 is the number of spacer of the mixing vanes type: 5 instead of 7. This may be a reason for the under prediction.

If one compare the results obtained for A2 and A3 (Figures 18 and 19), they are similar in the sense that there is a tendency for an under prediction. But the predictions for the bundle type A3 (Figure 19) are better than the ones for bundle type A2. The only difference between A2 and A3 is the number of rods: 25 and 36, respectively.

Figure 18: Measured and simulated critical power values for assembly A2.
Figure 19: Measured and simulated critical power values for assembly A3.

The axial power profile for bundle types A0, A2, and A3 is uniform which is not really representative for reactor conditions. The following test series with bundle type-A4, A8, A11, A12, and A13 were performed with cosine shaped axial power profile and hence are more close to real reactor conditions. Bundle types A4, A11, and A13 are repetition tests and they do not contain guide tubes, while A8 and A12 consist of 24 simulator rods and one guide tube.

The comparison of the SUBCHANFLOW predictions with the CHF data for these bundle types is good, especially for the test series A13 which is a repetition of the test series A4 (Figures 20 and 22). But also for the most realistic test configuration (A8) the comparison of predictions and data is quite acceptable (Figure 21).

Figure 20: Measured and simulated critical power values for assembly A4.
Figure 21: Measured and simulated critical power values for assembly A8.
Figure 22: Measured and simulated critical power values for assembly A13.

In summary, the following statement can be made: the best simulation result is achieved for cases A0, A4, A8, and A13 with standard deviations of 8%, 9%, 8%, and 5%. A4 and A13 have the same geometry and a cosine axial power shape. A0 and A8 are different in number of spacers, power profile and A8 contains a guide tube. In case of A2 and A3, the simulations underestimate the critical power on average by 10% and 8%. A3 is the only case with a 6 by 6 rod array. A2 has the same geometry as A4, but uses a uniform power shape instead of a cosine shape. Compared to A0, A2 has a different number of spacers, but the power profile is the same. We do not find a clear correlation of the quality of the results by comparing the geometry and power shapes. Finally, the EPRI correlation gives good results for 4 test cases and underestimates the results for two cases, which is conservative in the sense of reactor safety.

For several test runs, the axial and radial location of the occurrence of DNB is documented in the PSBT benchmark. SUBCHANFLOW detects first occurrence of DNB always at the 3 × 3-center rods. For 27% of the documented cases, DNB was experimentally found at peripheral rods. For the axial location the average simulation result is 0.05 m above the average measured value at 2.64 m of the heated length (3.66 m). There is a strong scattering of the simulation results compared to the measurements expressed by a standard deviation of 0.44 m.

3.6. Transient Tests under Critical Heat Flux Conditions

For the investigation of the DNB phenomena, four PWR relevant transient scenarios were carried out at the NUPEC PSBT facility: power increase (PI), flow reduction (FR), depressurization (DP), and temperature increase (TI). Assemblies A4 and A8 are used with an initial inlet temperature of about 291°C and a system pressure of 15.6 MPa. The mass flux at steady-state is about 3100 kg/(sm²). The steady-state power is 2.5 MW. The details of the transient behavior of important parameters are found in [16]. Table 10 shows the comparison of the transient time for the incidence of a critical heat flux condition and the corresponding critical power of the predictions (sim.) and the experiments (exp.).

Table 10: Occurrence of critical heat flux conditions in transient tests.

The maximum deviation between measured and predicted critical power is about 5%. The transient tests were conducted using bundle geometries which give good accordance between simulation and experiment for steady-state, too. So the results are consistent regarding the bundle type used. There is only one case (A8 PI) in which the simulated critical power is significantly larger than the measured one.

4. Summary and Outlook

The main features and validation effort for SUBCHANFLOW regarding PWR relevant phenomena were presented and discussed. A subchannel code based on the experience and empirical formulations of the last decades has been developed and validated using the PSBT benchmark data for typical bundle configurations used in pressurized light water reactors. In a first step, boiling in single subchannels was investigated to validate the basic empirical correlations used for boiling forced flow conditions. Furthermore, the single-phase flow turbulent mixing coefficient was derived from code predictions in comparison to outlet temperatures of a radial nonuniform heated bundle. The predicted void fractions at three axial levels for different bundle configurations and flow boundary conditions agree well with the steady-state and transient void measurements. DNB data for test bundles with uniform and cosine power shape were evaluated by adopting the EPRI critical heat flux model correlation. The standard deviation from the exact accordance of simulation and experiment at in maximum about 10%. For some cases, an underestimation of critical power is observed, which is conservative regarding reactor safety.

The performed investigations clearly demonstrate the prediction capability of SUBCHANFLOW which is an important pillar of the multiphysical and multiscale developments at KIT [2628]. The validation of SUBCHANFLOW related to BWR and innovative reactor phenomena is underway. The integration of this code in the NURESIM platform and the coupling to the reactor dynamic code COBAYA for both square and hexagonal geometries is also advanced [26].


𝐴Subchannel flow area (m²)
𝐷:Hydraulic diameter (m)
𝑓:Single-phase friction coefficient (empirical correlation)
𝑔:Gravity (m/s²)
𝐺:Mass flux (kg/(m²s))
:Specific mixture enthalpy (J/kg)
fg:Evaporation enthalpy (J/kg)
𝐾:Axial pressure loss coefficient (e.g.,) of spacers
𝐾𝐺:Lateral gap pressure loss coefficient (empirical constant)
𝑙:Distance of neighboring subchannels midpoints (m)
𝑚:Mass flow rate at axial cell boundary (kg/s)
𝑁:Number of measurements
𝑝:Pressure at axial cell boundary (Pa)
Δ𝑝:Pressure difference between neighboring channels (Pa)
𝑄:Linear power released to subchannel (W/m)
𝑠:Gap width between two neighboring rods (m)
𝑇exp:Measured temperature (C)
𝑇sim:Calculated temperature (C)
Δ𝑡:Time step (s)
𝑤:Linear mass flow rate through the gap (kg/(ms))
𝑤:Turbulent cross-flow (kg/(ms))
Δ𝑋:Length of axial cell (m)
𝑥:Steam quality
𝛼:Void fraction (empirical correlation, calculated from steam quality)
𝛽:Mixing coefficient (empirical constant)
𝜌:Density (kg/m³)
𝜎:Standard deviation
Φ2:Two-phase friction multiplier (empirical correlation)
Old:Value at previous time step
𝑖:Channel 𝑖
𝑗:Axial cell 𝑗
𝑘:Gap 𝑘
𝑛(𝑘):Channel neighbor belonging to gap 𝑘.


  1. M. J. Thurgood, et al., “COBRA/TRAC—A Thermal-Hydraulic Code for Transient Analysis of Nuclear Reactor Vessel and Primary Coolant Systems,” NUREG/CR-3046, 1983.
  2. C. Y. Payk, et al., “Analysis of FLECHT SEASET 163-Rod Blocked Bundle Data using COBRA-TF,” NRC/EPRI/Westinghouse-12, 1985.
  3. M. Avramova, Development of an innovative spacer grid model utilizing computational fluid dynamics within a subchannel analysis tool [Ph.D. thesis], Penn State University, 2007.
  4. K. Ivanov and M. Avramova, “Challenges in coupled thermal-hydraulics and neutronics simulations for LWR safety analysis,” Annals of Nuclear Energy, vol. 34, no. 6, pp. 501–513, 2007. View at Publisher · View at Google Scholar · View at Scopus
  5. D. Basile, R. Chierici, M. Beghi, E. Salina, and E. Brega, “COBRA-EN, an Updated Version of the COBRA-3C/MIT Code for Thermal-Hydraulic Transient Analysis of Light Water Reactor Fuel Assemblies and Cores,” Report 1010/1 (revised 1.9.99), ENEL-CRTN Compartimento di Milano.
  6. H. Blasius, “Das Ähnlichkeitsgesetz bei Reibungsvorgängen in Flüssigkeiten,” Forschungs-Arbeit des Ingenieur-Wesens 131, 1913 (in German).
  7. A. A. Armand, The Resistance During the Movement of a Two-Phase System in Horizontal Pipes, vol. 828 of AERE-Lib/Trans, Izvestiya Vsesojuznogo Teplotekhnicheskogo Instituta, 1946.
  8. P. W. Dittus and L. M. K. Boelter, “Heat transfer in automobile radiators of the tubular type,” University of California Publications in Engineering, vol. 2, no. 13, pp. 443–461, 1930, reprinted in International Communications in Heat and Mass Transfer, vol. 12, pp. 3—22, 1985.
  9. C. L. Wheeler, et al., “COBRA-IV-I: an Interim Version of COBRA for Thermal-Hydraulic Analysis of Rod Bundles Nuclear Fuel Elements and Cores,” BNWL-1962, Battelle, Pacific Northwest Laboratory, RichIand, Washington, DC, USA, 1976.
  10. R. W. Bowring, “Physical model based on bubble detachment and calculation of steam voidage in the subcooled region of a heated channel,” Tech. Rep. HPK-10, Institutt for Atomenergi, Halden, Norway, 1962.
  11. H. Dorra, S. C. Lee, and S. G. Bankoff, “A Critical Review of Predictive Models For The Onset of Significant Void in Forced-Convection Subcooled Boiling,” WSRC-TR-93-404, Westinghouse Savannah River Corporation, 1993.
  12. B. Chexal, G. Lellouche, J. Horowitz, and J. Healzer, “A void fraction correlation for generalized applications,” Progress in Nuclear Energy, vol. 27, no. 4, pp. 255–295, 1992. View at Scopus
  13. S.-H. Ahn and G.-D. Jeun, “Effect of spacer grids on CHF at PWR operating conditions,” Journal of the Korean Nuclear Society, vol. 33, no. 3, pp. 283–297, 2001.
  14. S. G. Beus, “A Two-Phase Turbulent Mixing Model for Flow in Rod Bundles,” Bettis Atomic Power Laboratory, WAPD-T-2438, 1970.
  15. M. Glück, “Validation of the sub-channel code F-COBRA-TF. Part II. Recalculation of void measurements,” Nuclear Engineering and Design, vol. 238, no. 9, pp. 2317–2327, 2008. View at Publisher · View at Google Scholar · View at Scopus
  16. A. Rubin, A. Schoedel, M. Avramova, H. Utsuno, S. Bajorek, and A. Velazquez-Lozada, “OECD/NRC Benchmark Based on NUPEC PWR Subchannel and Bundle Tests (PSBT), Volume I: Experimental Database and Final Problem Specifications,” US NRC OECD Nuclear Energy Agency, 2010.
  17. D. S. Rowe, “COBRA IIIC: A Digital Computer Program for Steady-State and Transient Thermal Analysis of Rod Cundle Nuclear Fuel Elements,” BNWL-1695, Pacific Northwest Laboratory, 1973.
  18. C. L. Wheeler, et al., “COBRA-IV-I: An Interim Version of COBRA for Thermal Hydraulic Analysis of Rod Bundle Nuclear Fuel Elements and Cores,” BNWL-1962, Pacific Northwest Laboratory, 1976.
  19. U. Imke, V. Sanchez, and R. Gomez, “SUBCHANFLOW: an empirical knowledge based subchannel code,” in Proceedings of the Annual Meeting on Nuclear Technology, pp. 4–6, Berlin, Germany, May 2010.
  20. V. Sánchez, U. Imke, and R. Gomez, “SUBCHANFLOW: a thermal-hydraulic sub-channel program to analyse fuel rod bundles and reactor cores,” in Proceedings of the 17th Pacific Basin Nuclear Conference, pp. 24–30, Cancún, México, October 2010.
  21. L. S. Tong, Boiling Heat Transfer and Two-Phase Flow, John Wiley & Sons, 1965.
  22. A. Berkhan, V. Sánchez, and U. Imke, “Validation of PWR relevant models of SUBCHANFLOW using the NUPEC PSBT Data,” in Proceedings of the 14th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-14), pp. 25–29, Toronto, Canada, September 2011.
  23. J. T. Rogers and R. G. Rosehart, “Mixing by turbulent interchange in fuel bundles. Correlations and interfaces,” in: ASME, 72-HT-53, 1972.
  24. J. T. Rogers and A. E. E. Tahir, “Turbulent interchange mixing in rod bundles and the role of secondary flows,” in ASME, Heat Transfer Conference, San Francisco, Calif, USA, 1975.
  25. D. G. Reddy and C. F. Fighetti, “Parametric Study of CHF Data Volume 2, A generalized subchannel CHF correlation for PWR and BWR fuel assemblies,” EPRI-NP-2609, 1983.
  26. M. Calleja, V. Sanchez, and U. Imke, “Implementation of SUBCHANFLOW in the SALOME Platform and Coupling with the Reactor Dynamic Code COBAYA3,” Jahrestagung Kerntechnik, Berlin, Germany, 17.19. Mai, 2011.
  27. A. Gomez, V. Sanchez, U. Imke, and R. Macian, “DYNSUB: A Best Estimate Coupled System for the Evaluation of Local Safety Parameter,” Jahrestagung Kerntechnik, Berlin, Germany, 17.19. Mai, 2011.
  28. Ivanov, V. Sanchez, and U. Imke, “Application of a coupled code system NCNP-SUBCHANFLOW to study conventional and innovative PWR fuel assembly types,” Jahrestagung Kerntechnik, Berlin, Germany, 17.19. Mai, 2011.