Validation of the Subchannel Code SUBCHANFLOWUsing the NUPEC PWR Tests ( PSBT )

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.


Introduction
The requirements for computational resources of highresolution 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 twofluid 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 SUB-CHANFLOW as a modern and modular subchannel code starting from the COBRA-family [5,17,18] using the abovementioned 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 Science and Technology of Nuclear Installations 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 steadystate 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.

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.

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, SUB-CHANFLOW 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 SUBCHAN-FLOW, 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.

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 ( j) in channel (i) surrounded by volumes of neighboring channels (n) through a gap (k), the basic conservation equations in finite difference form are the following.
Mass conservation: Energy conservation: ( Axial momentum: ( Lateral momentum: Slip correction for enthalpy derivative based on the twophase flow basic energy equation (Tong's function) [21]: ( Effective specific volume 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.

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].

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.
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:  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].
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.

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 SUB-CHANFLOW 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 i and channel j through gap kis defined by 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. 5 gives an overview about all test bundle geometries used in the present paper for void fraction and DNB (departure from nucleate boiling) simulations.

Steady-State Bundle Void Fraction Distribution. Table
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.
In addition to the correlations used for the simulation of the single-channel experiments, SUBCHANFLOW uses the closure laws summarized in Table 7.
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.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.

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 9-12.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.
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.
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).

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.
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.
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).
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.

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 2 ).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.).
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 10 Science and Technology of Nuclear Installations  used.There is only one case (A8 PI) in which the simulated critical power is significantly larger than the measured one.

Summary and Outlook
The main features and validation effort for SUBCHAN-FLOW 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 [26][27][28].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].Steam quality α: Void fraction (empirical correlation, calculated from steam quality) β: Mixing coefficient (empirical constant) ρ: Density (kg/m 3 ) σ: Standard deviation Φ 2 : Two-phase friction multiplier (empirical correlation) Old: Value at previous time step Liq: Liquid Vap: Vapor i: Channel i j: Axialcell j k: Gapk n(k): Channel neighbor belonging to gap k.

Figure 1 :
Figure 1: Comparison between calculated and measured void fraction for subchannel type S1.

Figure 2 :
Figure 2: Comparison between calculated and measured void fraction for subchannel type S2.

Figure 3 :
Figure 3: Comparison between calculated and measured void fraction for subchannel type S3.

Figure 4 :
Figure 4: Comparison between calculated and measured void fraction for subchannel type S4.

Figure 5 :
Figure 5: Comparison of calculated and measured channel outlet temperatures for a mixing coefficient of 0.06 (envelope: ±10 • C).

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

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

Figure 11 :Figure 12 :
Figure 11: Pressure decrease in the DP test.

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

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

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

Figure 20 :
Figure 20: Measured and simulated critical power values for assembly A4.

Figure 21 :
Figure 21: Measured and simulated critical power values for assembly A8.

Figure 22 :
Figure 22: Measured and simulated critical power values for assembly A13.

NomenclatureA: 11 K:
Subchannel flow area (m 2 ) D h : Hydraulic diameter (m) f : Single-phase friction coefficient (empirical correlation) g: Gravity (m/s 2 ) G: Mass flux (kg/(m 2 s)) h: Specific mixture enthalpy (J/kg) h fg : Evaporation enthalpy (J/kg) Science and Technology of Nuclear Installations Axial pressure loss coefficient (e.g.,) of spacers K G : Lateral gap pressure loss coefficient (empirical constant) l: Distance of neighboring subchannels midpoints (m) m: Mass flow rate at axial cell boundary (kg/s) N: Numberofmeasurements p: Pressure at axial cell boundary (Pa) Δp: Pressure difference between neighboring channels (Pa) Q: Linear power released to subchannel (W/m) s: Gap width between two neighboring rods (m) T exp : Measured temperature ( • C) T sim : Calculated temperature ( • C) Δt: Time step (s) w: Linear mass flow rate through the gap (kg/(ms)) w : Turbulent cross-flow (kg/(ms)) ΔX: Length of axial cell (m) x:

Table 1 :
Boundary conditions for single channel tests.

Table 3 :
Global result for different subchannel types.

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

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

Table 7 :
Empirical correlations used for the bundle tests.

Table 8 :
Void fraction standard deviations for all steady state test cases.

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