Research Article  Open Access
Comparative Analysis of CTF and Trace ThermalHydraulic Codes Using OECD/NRC PSBT Benchmark Void Distribution Database
Abstract
The international OECD/NRC PSBT benchmark has been established to provide a test bed for assessing the capabilities of thermalhydraulic codes and to encourage advancement in the analysis of fluid flow in rod bundles. The benchmark was based on one of the most valuable databases identified for the thermalhydraulics modeling developed by NUPEC, Japan. The database includes void fraction and departure from nucleate boiling measurements in a representative PWR fuel assembly. On behalf of the benchmark team, PSU in collaboration with US NRC has performed supporting calculations using the PSU inhouse advanced thermalhydraulic subchannel code CTF and the US NRC system code TRACE. CTF is a version of COBRATF whose models have been continuously improved and validated by the RDFMG group at PSU. TRACE is a reactor systems code developed by US NRC to analyze transient and steadystate thermalhydraulic behavior in LWRs and it has been designed to perform bestestimate analyses of LOCA, operational transients, and other accident scenarios in PWRs and BWRs. The paper presents CTF and TRACE models for the PSBT void distribution exercises. Codetocode and codetodata comparisons are provided along with a discussion of the void generation and void distribution models available in the two codes.
1. Introduction
In the past few decades, the need of improved nuclear reactor safety analyses has led to a rapid development of advanced methods for multidimensional thermalhydraulic analyses. These methods have progressively become more complex in order to account for variety of physical phenomena anticipated during steadystate and transient Light Water Reactor (LWR) conditions. The newly developed models must be extensively validated against fullscale highquality experimental data. The previous publically available void distribution measurements, which include the ISPRA 16rod tests [1] and the GE 9rod tests [2], have limited databases. Currently the requirements to the numerical modelling of subchannel void distribution dictate an approach that can be applied to a wide range of geometrical and operating conditions. In the past decade, experimental and computational technologies have tremendously improved the study of the flow structure. In that sense, the new OECD/NRC PWR Subchannel and Bundle Tests (PSBT) benchmark [3] has provided an excellent opportunity for validation of innovative models for void distribution and departure from nucleate boiling (DNB) prediction under Pressurized Water Reactors (PWRs) conditions. From 1980s to 1990s, NUPEC (Nuclear Power Engineering Corporation) performed a series of void measurement tests using fullsize mockup tests for both Boiling Water Reactors (BWRs) and PWRs. Based on stateoftheart computer tomography (CT) technology, the void distribution was visualized at the mesh size smaller than the subchannel under actual plant conditions. NUPEC also performed steady state and transient critical power test series based on the equivalent fullsize mockups. Considering the reliability not only of the measured data, but also other relevant parameters such as the system pressure, inlet subcooling, and rod surface temperature, these test series supply the first substantial database for the development of truly mechanistic and consistent models for void distribution and departure from nucleate boiling. The PSBT benchmark was divided into two separate phases, with each consisting of individual exercises. Phase I is Void Distribution Benchmark while Phase II is Departure from Nucleate Boiling Benchmark. The benchmark specification is designed so that it can systematically assess and compare the participants’ numerical models for the prediction of detailed subchannel void distributions and departure from nucleate boiling to full scale experimental data on a prototypical PWR rod bundle.
CTF is a version of the COBRATF code maintained at the Reactor Dynamics and Fuel Management Group (RDFMG) at Pennsylvania State University (PSU) [4]. The original version of COBRATF was developed at the Pacific Northwest Laboratory as a part of the COBRA/TRAC thermalhydraulic code. 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 as well. The code version used at PSU originates from a code version modified during the FLECHT SEASET program [5]. Besides the code utilization to teach and train students in the area of nuclear reactor thermalhydraulic safety analyses, during the last few years at PSU the theoretical models and numerics of COBRATF were substantially improved [6]. The code was subjected to an extensive verification and validation program and was applied to variety of LWR steadystate and transient simulations. CTF is a transient code based on a separated flow representation of the twophase flow. The twofluid formulation, often used in thermalhydraulic 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 timeaveraged 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 semiimplicit, finitedifference 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 hot wall and normal flow regime maps and it is capable of calculating reverse flow, counter flow, and crossflow situations. The code is developed for use with either threedimensional (3D) Cartesian or subchannel coordinates and, therefore, it features extremely flexible nodding for both the thermalhydraulic and the heattransfer solutions. This flexibility allows a fully 3D treatment in geometries amenable to description in a Cartesian coordinate system.
TRACE is a multicomponent solver consolidation of four US NRC computer codes: TRACP, TRACB, RELAP5, and RAMONA. TRACE has been validated and assessed against more than 500 experimental sets of data from separate and integral effect tests, in which comparisons were found to be reasonable in general [7]. TRACE utilizes a finitevolume technique to discretize typical hydraulic components found in a nuclear power plant and calculates the internal energy and equations of motion in each component for two phases. The energy equation is solved using a semiimplicitly numerical scheme and the equations of fluid motion are solve using the stabilityenhancing twostep (SETS) numerical scheme which allows the material Courant limit to be exceeded. This allows very large time steps to be used in slow transients. This set of equations is solved for one and three dimensions in Cartesian and/or cylindrical coordinates. Errors introduced to the solution due to abrupt area changes are corrected by modifying the equations of motion to force Bernoulli flow. For the analysis and results presented in this paper the authors utilized TRACE versionV5.02p2_win32_opt.
The paper presents the CTF and TRACE models for the exercises of the void distribution phase of the OECD/NRC PSBT benchmark. Codetocode and codetodata comparisons are provided along with a discussion of the void generation and void distribution models available in the two codes. The following two sections discuss the void generation and distribution models available in CTF and TRACE with a subsequent codetocode and codetodata comparisons.
2. CTF Models for Vapor Generation and Distribution
The threefield formulation of the twophase flow used in CTF is a straightforward extension of the general twofluid model. Dividing the liquid phase into a continuous liquid field and an entrained liquid drop field allows both fields to have different velocities. The generalized phasic momentum equation is then given as where is the average phase void fraction; is the average kphase density; is the average kphase velocity vector; is the acceleration of gravity vector; is the average kphase viscous stress tensor; is the average supply of momentum to phase k due to mass transfer to phase k; is the average drag force on phase k by the other phases; is the average supply of momentum to phase k due to turbulent mixing and void drift.
In the generalized phasic momentum equation the terms representing the momentum exchange at the interface (interfacial momentum terms) are expressed aswhere is the average drag force per unit volume by the vapor on the continuous liquid and is the average drag force per unit volume by the vapor on the entrained liquid.
The momentum exchange due to mass transfer between the three fields can be written aswhere the is the average rate of vapor generation per unit volume and is the average net rate of entrainment per unit volume. Since both liquid fields contribute to the vapor generation, then .
If denotes the fraction of the total vapor generation coming from the entrained liquid field, then The momentum exchange due to turbulent mixing and void drift is neglected in the entrained liquid field in the annular flow regime: Also, the viscous stress is partitioned into a wall shear and a fluidfluid shear; the fluidfluid shear is neglected: The model for interfacial mass transfer is obtained from the energy jump condition by neglecting the mechanical terms and averaging The interfacial heat transfer, , for phase is given by where is the average interfacial area per unit volume and is a surface heat transfer coefficient. The vapor generation is divided into four components, two for each phase, depending on whether the phase is superheated or subcooled and the total vapor generation rate is given by the sum of these components. The interfacial area per unit volume, , is based on the flow regime, as are the heat transfer coefficients, . The interfacial drag force per unit volume between any two fields is assumed to be a function of the relative velocity between both fields. The interfacial friction coefficients are flow regime dependent and, therefore, neither void correlation nor twophase pressure drop correlation has to be applied. Interfacial drag forces are modeled between continuous liquid and disperse vapor in the bubbly flows and between continuous liquid film and vapor core and entrained droplets and vapor core in the annular flow. The treatment of the interfacial drag is described in Table 1.

Turbulent mixing and void drift phenomena are modeled in CTF by the Lahey and Moody approach [8], where the net twophase 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 singlephase mixing coefficient can be either specified as an input value or calculated using an empirical correlation derived by Rogers and Rosehart [9]. The Beus’ model for twophase turbulent mixing is utilized [10]. In 1980s, both approaches were representing the stateoftheart in turbulent mixing and void drift modeling. Nowadays they are still used in the most of the subchannel codes. A detailed description of the current CTF turbulent mixing and void drift models is given in Table 2.

3. TRACE Model Description
The fully conservative forms of the energy and momentum equations are modified in TRACE to provide a set of internal energy and motion equations. This modification reduces the numerical manipulation and computational time of the solution. This modification is also transferred to the conservation of mass equation.
It is assumed that the volume average of a product is equal to the product of volume averages. Only contributions from wall heat fluxes and heat fluxes at phase interfaces within the averaging volume are normally included in the volume average of the divergence of heat flux. Also, only contributions from the stress tensor due to shear at metal surfaces or phase interfaces within the averaging volume are considered. The only portions of the work terms that contribute to change in bulk kinetic energy of motion are retained excluding viscous heating from most of the cases unless a pump component is used in which case the viscous heating from the pump to the fluid is incorporated by the term of direct heating in the internal energy equation.
These modifications and assumptions yield a set of 6 equations of mass (9) and (10), motion (11) and (12), and internal energy (13) and (14) for gas and gasliquid mixture. An additional mass equation is added for noncondensable gases but in order to still solving only a single set of motion and energy equations, the nongases are assumed to be in mechanical and thermal equilibrium with the steam.
Mass:
Motion:
Internal energy: In (9) and (10), the term is the void fraction; and are the density of the gas and liquid, respectively; and the velocity vectors of gas and liquid; and is the interfacial masstransfer rate (positive from liquid to gas).
In (11) and (12) the additional term is the fluid or total pressure; is the force per unit volume due to shear at the phase interface; is the wall shear force per unit volume acting on the liquid; is the wall shear force per unit volume acting on the gas; is the flow velocity at the phase interface; and is the gravity vector.
The other terms on (13) and (14) are and which are the internal energies of the gas and liquid, respectively. The terms and are the heat transfer rates per unit volume from the wall to gas and from the wall to liquid. The terms and correspond to the power deposited directly to the gas or liquid (without heatconduction process). The term is the interfacial sensible heat transfer. The term accounts for energy carried with mass transfer at the interface, which is the product of mass transfer rate and appropriate stagnation enthalpy at the interface.
The phasechange rate in the set of equations is calculated using the heat conduction limited model: where The term in (16) and (17) is the interfacial area per unit volume, where and are the heat transfer coefficients at the liquid/gas interface and the saturation temperature corresponding to the partial pressure .
The interfacial drag force incorporated in the motion equations (11) and (12) is evaluated by (18). The interfacial drag force is evaluated for vertical pipes and for horizontal/inclined pipes. For vertical pipes the set of correlations is calculated for Pre and Postcriticalheatflux (Pre and PostCHF) condition and for where is the interfacial drag coefficient and the relative velocity: where is the velocity of the gas phase and is the velocity of the liquid phase.
The velocity of the gas phase is evaluated using the local drift velocity (20), where is the volumetric flux: For flow in vertical pipes under PreCHF conditions the interfacial drag coefficient is calculated with (21) and the profile factor with (22) subsequently: A drift flux model approach is used to evaluate local drift velocity () along with the distribution coefficient (). Table 3 summarizes the actual drift flux models used in TRACE for small and large pipes and Bubbly/Slug and the Annular/Mist flow regimes under PreCHF condition.

For PostCHF conditions three principal inverted flow regimes are modeled in TRACE, which are inverted annular, inverted slug, and dispersed flow. These three regimes are defined in terms of void fraction and gas superficial velocity.
The inverted annular regime is used in TRACE for void fractions below 0.6 and the interfacial drag coefficient is calculated using Inverted slug regime is used in TRACE for a void fraction between 0.6 and 0.9 and the dispersed flow for a void fraction over 0.9. In both regimes the interfacial drag coefficient is calculated with (24), where is the drag coefficient for multiparticle and is the projected area per unit volume. The projected area is calculated differently for each regime:
4. Description of Phase I of the OECD/NRC PSBT Benchmark
The first phase of the PSBT benchmark [3] was intended to provide data for the verification of void distribution models in participants’ codes. This phase was composed of four exercises: a steadystate single subchannel benchmark (Exercise I1), a steadystate rod bundle benchmark (Exercise I2), a transient rod bundle benchmark (Exercise I3), and a pressure drop benchmark (Exercise I4). The results presented in this paper are for the steadystate single subchannel benchmark, the steadystate rod bundle benchmark, and the transient rod bundle benchmark.
The range of operating conditions for the facility is given in Table 4 and the operating conditions for the four transient scenarios are given in Table 5.


The properties of each subchannel assembly (Exercise I1) are given in Table 6.
 
Dark square: subjected subchannel; White circle: heated rod; Dark circle: thimble rod. 
The properties of the bundle assemblies (Exercise I2) to be used are given in Table 7.
 
White circle: heated rod; Dark circle: thimble rod; MV: mixing vane; NMV: no mixing vane. Spacer location is distance from bottom of heated length to spacer bottom face. 
Four transient scenarios (temperature increase, power increase, depressurization, and flow reduction) were used in Exercise I3 for each test series, yielding twelve total test cases. The test conditions are summarized in Table 8.

Electrical heaters are used in the PSBT experiments and the heater rod data is specified in the Benchmark Specifications [3] along with radial and axial power distributions to be used for each test simulation.
5. CTF and TRACE Applications to the Void Distribution Phase of the OECD/NRC PSBT Benchmark
The test cases of Exercise I1 were calculated with CTF for four bundle types—S1, S2, S3, and S4. Only the heated length of the subchannel was modeled in an axial discretization of forty equidistant nodes. Codetodata comparisons are given in Figure 1. It can be seen that the CTF predictions stay within the error bound of 10% void (the experimental uncertainties for the steady state void fraction CT scanner measurements were specified as 3% void). This is in agreement with a previously observed tendency of CTF to overpredict the vapor generation rate [11], which is due the utilized interfacial drag modeling in CTF.
(a)
(b)
(c)
(d)
Eight tests of the steadystate series5 ( bundle B5) were modelled by TRACE and CTF. The heated section of the bundle is modeled in TRACE with a threedimension component discretized in 23 axial nodes, 2 radial nodes, and 1 azimuthally node assuming that the power distribution is axis symmetrical. As it can be observed in Figure 2, TRACE predicted the void fraction at the upper part of the bundle with an average error of 2% and maximum error of 7%. In the middle part of the bundle TRACE predicted the void fraction with an average error of 8% and a maximum of 13%. On the other hand, for the lower part of the bundle TRACE overpredicted the void fraction with an average error of 10% and a maximum 16%. The entire B5 bundle was modelled by CTF in a subchannelbysubchannel basis—no symmetry was used. The heated length was divided axially into seventy equidistant nodes. The pressure losses due to spacer grids were calculated as velocity head losses with a loss coefficient of 1.0. The total crossflow between two adjacent subchannels was simulated as a sum of the diversion crossflow due to lateral pressure gradients and the lateral flow due to turbulent mixing and void drift. The steadystate void fraction predictions by CTF show very similar, but slightly better agreement with the measurements as compared to TRACE. The TRACE results obtained from V5.02p2_win32_opt version for this study are similar to the results obtained from other participants in the benchmark and using RELAP5 and TRACE codes [12].
(a)
(b)
(c)
(d)
(e)
(f)
Power increase, flow reduction, depressurization, and temperature increase transients were simulated by NUPEC and selected as benchmark exercise cases. The spaceaveraged instantaneous axial void fraction profiles during the transients were supplied for codetodata comparisons. The Xray densitometer measurements were taken at three intermediate elevations along the heated lengths: 2216 mm, 2669 mm, and 3177 mm. The four transients were simulated with CTF and TRACE for the bundle type B5. Both codes utilised the same configurations used in the steadystate cases. As previously mentioned in TRACE the heated length was divided in 23 axial nodes, where 17 of those nodes upperfaces are located at the same elevation of the spacer grids, in which pressure drops are incorporated into the model with a k factor of 1 as well. CTF and TRACE results are given in Figure 3. The experimental uncertainties were specified as 5% void. As seen in Figures 3, 4, 5, and 6, both codes are capable of reproducing the transient behavior of the bundle average void fraction for the four transient scenarios. The agreement is better at higher axial elevations.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
The time shift observed in the temperature increased transient for both codes should be attributed to the heat capacitance effect of the downcomer region. A study was performed by another benchmark participant, CEAGrenoble [13], regarding the heat capacitance effect of the downcomer region in the bundle test section and the location of the fluid temperature measurement. The fluid temperature measurement occurred just before the coolant inlet nozzle above the downcomer region, so it can be reasonably assumed that there would be some time shift in the flow characteristics due to the time required for the flow to traverse the downcomer. It should be noted that this effect appears to only be significant for the temperature increased transients. The CEAGrenoble team investigated the effect of the downcomer region and determined that the fluid temperature is reduced and shifted when the downcomer is accounted for. A time shift of 10 seconds for the experimental void fraction data was recommended. In addition, it should be noted that DNB will occur earlier when modeling the downcomer region than the experimental data suggests.
6. Conclusions
On behalf of the OECD/NRC PSBT benchmark team, PSU in collaboration with US NRC has performed supporting calculations of the benchmark exercises using its inhouse advanced thermalhydraulic subchannel code CTF and the US NRC system code TRACE. CTF and TRACE were applied to the steadystate and transient void distribution cases. Both codes were able to reproduce the measured data in a good agreement.
Recently, the need to refine models for best estimate calculations based on goodquality experimental data has arisen for various nuclear applications. One of the most extensive and valuable databases available was developed by the Nuclear Power Engineering Corporation (NUPEC) of Japan, consisting of both void distribution and departure from nucleate boiling (DNB) data for a representative pressurized water reactor (PWR) fuel assembly and these data were assembled in the PSBT benchmark.
The objective of the benchmark was twofold. First, the benchmark aimed to evaluate currently available computational approaches in an effort to understand the strengths and weaknesses of current thermalhydraulic codes. Second, the benchmark was intended to encourage the development of the next generation of approaches that focus more on microscopic processes.
In the view of the above statements the PSBT data will be incorporated in the Verification and Validation matrices of both CTF and TRACE codes. The performed comparative analysis in the future will be complemented with uncertainty and sensitivity analysis, which can be used for further improvement of CTF and TRACE models.
References
 H. Herkenrath, W. Hufschmidt, U. Jung, and F. Weckermann, “Experimental investigation of the enthalpy and mass flow distribution in 16rod clusters with BWRPWR geometries and conditions,” EUR 7575 EN, ISPRA, 1981. View at: Google Scholar
 R. T. Lahey Jr., B. S. Shiralkar, and D. W. Radcliffe, TwoPhase Flow and Heat Transfer in MultiRod Geometries: Subchannel and Pressure Drop Measurements in a NineRod Bundle for Diabatic and Adiabatic Conditions, 1970, GEAP13049, GE.
 A. Rubin et al., “OECD/NRC benchmark based on NUPEC PWR subchannel and bundle tests (PSBT). Volume I: experimental database and final problem specifications,” NEA/NSC/DOC, 2010. View at: Google Scholar
 “CTF—a thermalhydraulic subchannel code for LWRs transient analyses. User’s manual,” Technical Report, RDFMG, The Pennsylvania State University, 2009. View at: Google Scholar
 C. Y. Payk et al., “Analysis of FLECHT SEASET 163rod blocked bundle data using COBRATF,” Tech. Rep. NRC/EPRI/Westinghouse12, 1985. View at: Google Scholar
 M. Avramova, D. Cuervo, K. Ivanov et al., “Improvements and applications of COBRATF for standalone and coupled LWR safety analyses,” in Proceeding of the International Conference on the Physics of Reactors (PHYSOR '06), Vancouver, Canada, September 2006. View at: Google Scholar
 TRACE V5. 0 Theory Manual, “Field Equations, Solution Methods, and Physical Models,” USNRC, Washington DC. View at: Google Scholar
 R. T. Lahey and F. J. Moody, The Thermal Hydraulics of a Boiling Water Nuclear Reactor, American Nuclear Society (ANS), 1993.
 J. T. Rogers and R. G. Rosehart, “Mixing by turbulent interchange in fuel bundles, correlations and inferences,” ASME 72HT53, 1972. View at: Google Scholar
 S. G. Beus, “A twophase turbulent mixing model for flow in rod bundles,” Tech. Rep. WAPDT2438, Bettis Atomic Power Laboratory, 1970. View at: Google Scholar
 M. Avramova, K. Ivanov, and L. E. Hochreiter, “Analysis of steady state and transient void distribution predictions for phase I of the OECD/NRC BFBT Benchmark using CTF/NEM,” in Proceedings of the 12th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH12 '07), Pittsburgh, Pa, USA, October 2007. View at: Google Scholar
 M. Avramova, A. Rubin et al., “OECDNEA/USNRC/NUPEC PWR Subchannel and bundle test (PSBT) benchmark, Volume II: final results of phase I on void distribution,” OECD/NEA Report, 2011. View at: Google Scholar
 M. Valette and C. E. A. Grenoble, Private Correspondence.
Copyright
Copyright © 2013 M. Avramova et al. 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.