Comparative Analysis of CTF and Trace Thermal-Hydraulic Codes Using OECD / NRC PSBT Benchmark Void Distribution Database

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 thermal-hydraulics 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 in-house advanced thermalhydraulic subchannel code CTF and the US NRC system code TRACE. CTF is a version of COBRA-TF 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 steady-state thermal-hydraulic behavior in LWRs and it has been designed to perform best-estimate 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. Code-to-code and code-to-data comparisons are provided along with a discussion of the void generation and void distribution models available in the two codes.


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 thermal-hydraulic analyses.These methods have progressively become more complex in order to account for variety of physical phenomena anticipated during steady-state and transient Light Water Reactor (LWR) conditions.The newly developed models must be extensively validated against full-scale high-quality experimental data.The previous publically available void distribution measurements, which include the ISPRA 16-rod tests [1] and the GE 9-rod 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 full-size mock-up tests for both Boiling Water Reactors (BWRs) and PWRs.Based on stateof-the-art 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 full-size 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 Science and Technology of Nuclear Installations 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 COBRA-TF code maintained at the Reactor Dynamics and Fuel Management Group (RDFMG) at Pennsylvania State University (PSU) [4].The original version of COBRA-TF 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 thermal-hydraulic safety analyses, during the last few years at PSU the theoretical models and numerics of COBRA-TF were substantially improved [6].The code was subjected to an extensive verification and validation program and was applied to variety of LWR steady-state and transient simulations.CTF is a transient code based on a separated flow representation of the two-phase flow.The two-fluid formulation, often used in thermal-hydraulic codes, separates the conservation equations of mass, energy, and momentum to vapor and liquid.CTF extends this treatment to three fields: vapor, continuous liquid, and entrained liquid droplets, which results in a set of nine time-averaged conservation equations.The conservation equations for each of the three fields and for heat transfer from and within the solid structure in contact with the fluid are solved using a semi-implicit, finite-difference numerical technique on an Eulerian mesh, where time intervals are assumed to be long enough to smooth out the random fluctuations in the multiphase flow, but short enough to preserve any gross flow unsteadiness.The code is able to handle both hot wall and normal flow regime maps and it is capable of calculating reverse flow, counter flow, and cross-flow situations.The code is developed for use with either three-dimensional (3D) Cartesian or subchannel coordinates and, therefore, it features extremely flexible nodding for both the thermal-hydraulic and the heat-transfer 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: TRAC-P, TRAC-B, 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 finite-volume 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 semi-implicitly numerical scheme and the equations of fluid motion are solve using the stability-enhancing two-step (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 version-V5.02p2win32 opt.
The paper presents the CTF and TRACE models for the exercises of the void distribution phase of the OECD/NRC PSBT benchmark.Code-to-code and code-to-data 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 code-to-code and code-to-data comparisons.

CTF Models for Vapor Generation and Distribution
The three-field formulation of the two-phase flow used in CTF is a straightforward extension of the general two-fluid 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 k-phase void fraction;   is the average k-phase density;   is the average k-phase velocity vector;  is the acceleration of gravity vector;   is the average k-phase 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 as where   ,vap liq is the average drag force per unit volume by the vapor on the continuous liquid and   ,vap ent 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 as where 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 Γ  = Γ  liq + Γ  ent .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 fluid-fluid shear; the fluid-fluid 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, h.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 two-phase 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 two-phase mixing (including void drift) is assumed to be proportional to the nonequilibrium void fraction gradient.The void drift is only assumed to occur in bubbly, slug, and churn flow, where liquid is the continuous phase and vapor is the dispersed phase.The single-phase mixing coefficient can be either specified as an input value or calculated using an empirical correlation derived by Rogers and Rosehart [9].The Beus' model for two-phase turbulent mixing is utilized [10].In 1980s, both approaches were representing the state-of-theart 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.

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 gas-liquid 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 heat-conduction 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 phase-change 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 Postcritical-heat-flux (Pre and Post-CHF) 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 Pre-CHF 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 Pre-CHF condition.
For Post-CHF 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 C , is the drag coefficient for

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 steady-state single subchannel benchmark (Exercise I-1), a steady-state rod bundle benchmark (Exercise I-2), a transient rod bundle benchmark (Exercise I-3), and a pressure drop benchmark (Exercise I-4).The results presented in this paper are for the steady-state single subchannel benchmark, the steady-state 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 I-1) are given in Table 6.
The properties of the bundle assemblies (Exercise I-2) to be used are given in Table 7.
Four transient scenarios (temperature increase, power increase, depressurization, and flow reduction) were used in Exercise I-3 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.

CTF and TRACE Applications to the Void Distribution Phase of the OECD/NRC PSBT Benchmark
The test cases of Exercise I-1 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.Code-to-data 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.
Eight tests of the steady-state series-5 (5 × 5 bundle B5) were modelled by TRACE and CTF.The heated section of the bundle is modeled in TRACE with a three-dimension 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 subchannel-by-subchannel 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 cross-flow between two adjacent subchannels was simulated as a sum of the diversion cross-flow 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].
Power increase, flow reduction, depressurization, and temperature increase transients were simulated by NUPEC and selected as benchmark exercise cases.The space-averaged instantaneous axial void fraction profiles during the transients were supplied for code-to-data comparisons.The X-ray 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 steady-state cases.As previously mentioned in TRACE the heated length was divided in 23 axial nodes, where 17 of those nodes upper-faces 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.
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, CEA-Grenoble  [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 CEA-Grenoble 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.

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 in-house advanced thermal-hydraulic subchannel code CTF and the US NRC system code TRACE.CTF and TRACE were applied to the steady-state 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 good-quality 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 thermal-hydraulic 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.

Figure 1 :
Figure 1: CTF predictions of steady-state void fraction in single-subchannel tests.
Middle level (2.669 m) steady-state test series 5 Middle level (2.669 m) steady-state test series 5

Figure 2 :
Figure 2: TRACE and CTF predictions of steady-state test series 5 of the 5 × 5 bundle B5.

Figure 4 :
Figure 4: CTF (a) and TRACE (b) predictions of void fraction in bundle type B5 during power increase transient.

Figure 5 :
Figure 5: CTF (a) and TRACE (b) predictions of void fraction in bundle type B5 during flow depressurization transient.

Figure 6 :
Figure 6: CTF (a) and TRACE (b) predictions of void fraction in bundle type B5 during temperature increase transient.

Table 4 :
Range of NUPEC PWR test facility operating conditions.

Table 5 :
Transient parameters of NUPEC PWR test facility.

Table 8 :
Test conditions for transient void measurement test series 5T, 6T, and 7T.