Implementation and validation of CHF-models in the two-phase porous-media code TWOPORFLOW Nuclear Engineering and Design

The local deterioration of the heat transfer at certain elevation of the fuel rods of a water-cooled reactor may occur either in case of dry-out of the liquid film (in case of BWR) or due to Departure from Nucleate Boiling (DNB, in case of PWR). Under such local conditions, a sudden cladding temperature escalation may take place challenging the integrity of the cladding material i.e., of a safety barrier. For a reliable evaluation of the safety features of a nuclear reactor under nominal or accidental conditions, the use of validated numerical tools using relevant experimental data is mandatory. In this paper, the implementation of different correlations (Biasi, Bowring, and Groeneveld Look-Up Table) for the prediction of the Critical Heat Flux (CHF) in TWOPORFLOW as well as its validation is presented and discussed. Based on the performed investigations it can be stated that the implemented correlations predict the CHF-value and their location of appearance in acceptable agreement to the measured data.


Introduction
The local deterioration of the heat transfer at certain elevation of the fuel rods of a water-cooled reactor may occur either in case of dry-out of the liquid film (in case of BWR) or due to Departure from Nucleate Boiling (DNB, in case of PWR). Under such local conditions, a sudden cladding temperature escalation may take place challenging the integrity of the cladding material i.e., of a safety barrier. For a reliable evaluation of the safety features of a nuclear reactor under nominal or accidental conditions, the use of validated numerical tools using relevant experimental data is mandatory.
System thermal hydraulic codes like TRACE (USNRC, 2013), RELAP5 (INL, 2018), ATHLET (Austregesilo et al., 2016), CATHARE (Préa et al., 2020) and thermal hydraulic subchannel codes such as COBRA-TF (Salko and Avramova, 2015), Subchanflow (Sánchez-Espinoza et al., 2010), VIPRE (Y. Sung, P. Schueren and A. Meliksetian, 1999), FLICA (Toumi et al., 2000), are widely used for the analysis of the LWRs core under nominal and accidental conditions either as stand-alone or coupled to neutron physics solvers. Considering that, the development of two-phase CFD-tools is still not mature and that the coarse-mesh system thermal hydraulic codes have some model deficiencies, investigations at KIT are focused on the development of the porous-media two-phase 3D code TWOPORFLOW (Jimenez et al., 2014) to fill the gap. It solves a system of six conservation equations for a 3D Cartesian mesh for a two-phase flow based on the porous-media approach using the implicit continuous Eulerian (ICE) method (Imke, 2004). Similar developments are followed in Finland, see PORFLO (Ilvonen et al., 2010), and South Korea CUPID (Jeong et al., 2010). This kind of codes appear as a promising alternative, which is more accurate than 1D or 3D coarse mesh codes and at the same time, it requires much less computational resources than CFD-tools.
The source code and data structure of TWOPORFLOW is continuously improved. The physical models like turbulent viscosity, turbulent heat transport, and void drift for 3D simulations of reactor cores were added (Jauregui . Models related to the pre-CHF heat transfer are validated using void fraction test data (PSBT NUPEC (Jauregui ) in the frame of a doctoral thesis. This paper is focused on the implementation of three CHFcorrelations and their validation using both tube and bundle test data obtained at stationary and transient conditions. The code is capable of simulating the different pre-CHF heat transfer regimes, i.e., forced convection, sub-cooled boiling and saturated boiling. The correlations of Biasi (Biasi et al., 1967) and Bowring (Bowring, 1972) based on tube tests were implemented in the basic code version and in (Jauregui Chavez et al., 2017) a preliminary validation was performed. Later on, the Groeneveld Look-Up- Table was implemented (Groeneveld et al., 2006), which is widely used also in other codes as for example in (USNRC, 2013;INL, 2018;Austregesilo et al., 2016;Préa et al., 2020), and (Salko and Avramova, 2015).
A detailed description of the conservation equations, constitutive equations and numerics of the TWOPORFLOW code is given in (Jauregui . Hence, it is not provided in this paper. The accurate prediction of the critical heat flux of reactor cores is very important when performing safety evaluations. This paper starts with the phenomenology description of the CHF-phenomena in Section 2. A brief presentation of the implemented correlations and their validity range of applications is provided in Section 3. The prediction capabilities of TWOPORFLOW using the CHF-correlations is discussed in Section 4 (40 Becker tube tests (Becker et al., May 1983), 151 steady state and 2 transient tests from the NUPEC BWR Full-size Fine-mesh Bundle Test (BFBT) (Utsuno et al., 2006). Section 5 summarizes the main conclusions and presents the outlook.

Phenomenology of critical heat flux
The CHF is a local physical phenomenon leading to a substantial deterioration of the heat-transfer coefficient (Todreas and Kazimi, 1993). In PWR, the Departure of Nucleate Boiling (DNB) is characterized by the appearance of a vapor film caused by an increased heat flux resulting in lowering the heat transfer coefficient. In BWR, the Dry-out is characterized by the local evaporation of the liquid annular film, Fig. 1, with similar consequences as the DNB.
The CHF phenomenon depends strongly on geometry (pipes, annular flow, rod bundles etc.) and is a very complicated process (Yang et al., 2021). In the present paper, it is described by a more general empirical approach to cover a wide application.
For the used code, the detection of local phenomena is limited by the relative coarse mesh of the porous media model.
The parameter defining the safety margins in nuclear reactors is the ratio between the critical heat flux and the actual calculated heat flux (Todreas and Kazimi, 1993): The minimum value of CHFR (MCHFR or MDNBR) should not be less than for example 1.3 to assure the integrity of the cladding. The specific value depends on the reactor type and the degree of conservatism in the evaluation.

CHF correlations implemented in TWOPORFLOW
Three different CHF-correlations namely the Biasi, Bowring, and Groeneveld Look-Up Table were implemented. They were chosen due to their broad application range not connected to a specific reactor type. In addition, an iterative procedure for the prediction of the surface temperature at CHF conditions was programmed. The physical parameter application ranges are summarized in Table 1.

The Biasi correlation
The CHF Biasi correlation (Biasi et al., 1967) was developed for round ducts and uniform heating with a root-mean-square (rms) error of 7.26% in 4551 data points. 85.5% of the data points are within ± 10% absolute deviation. The following equations are used:  Table  g Vapor phase  Bowring where H(P) and F(P) are functions of the pressure: n is a coefficient dependent of the hydraulic diameter and defined as:

Bowring correlation
The Bowring correlation (Bowring, 1972) was developed also for round tubes with a uniform axial heat flux. The rms error in this correlation is around 7% for 3800 data points. A wide range of application in terms of pressure and mass flux characterizes it. The correlation is described by the following equations: where: and: (11) F 1 to F 4 as well as n are functions of the reduced pressure: P is the pressure given in MPa. The value of n is calculated as: The F coefficients depend on pressure. Forp r between 0.98 and 1.02 they are set to one. For p r < 0.98MPa we have: and: For p R > 1.02MPawe have: and:

The Groeneveld Look-Up-Table
The 2006 Groeneveld Look-Up- Table (LUT) (Groeneveld et al., 2006) is a normalized database for the CHF-prediction depending on various parameters. It contains around 25; 000 data points taken from correlations derived from different experiments, and extrapolations. The rms error at constant inlet flow conditions is 5.86%. The Groeneveld CHFmodels used in TWOPORFLOW are given hereafter: For values of equilibrium quality below − 0.5, a linear extrapolation is made.
where K 1 is a correction factor for the hydraulic diameter taken from (Groeneveld et al., 2005): p factor is calculated if the pressure is larger than 21 MPa else it is set to one.
where the non-dimensional pressure is defined as:

The Becker CHF tests at steady state conditions
In total 14 exercises of a set of experiments done by Becker (Becker et al., May 1983) at the Royal Institute of Technology (KTH) in Stockholm were used for code validation. The test section is uniformly heated along 7000 mm. The temperature measurements were made with thermocouples positioned every 20 cm starting from the bottom up to 3 m high. Above it, the measurements were made every 10 cm. The tube is characterized by a hydraulic diameter of 0.01 m. Details about the tests can be found in (Becker et al., May 1983). Hereafter, the boundary conditions of the tests are listed. The TWOPORFLOW-model of the test is a one-dimensional pipe. TWOPORFLOW-simulations are carried out using the three different CHF-models. It is worth to mention that an experimental error in the heat flux of the order of 1% is reported. Fig. 2 shows the comparison of the CHF predicted by the code with the different correlations and the measured data as function of the mass flux for the different tests. It can be observed that the CHF predicted with the correlations of Bowring and Groeneveld are close to the experimental data for all mass flow rates. The accuracy of the CHF predicted by the Biasi-correlation is dependent on the mass flux: the over-prediction increases with increasing mass fluxes (lower qualities) in most cases. In summary, Biasi shows a CHF standard deviation from the experimental data of 9.19%, Bowring 1.30% and Groeneveld 1.43%. The higher deviations of the Biasi predictions are in the lower steam quality region.

The BFBT CHF tests at steady state conditions
The BFBT critical power tests series includes around 111 experiments for a test section consisting of a bundle with 8 × 8-pin arrangement and a central water rod. The axial heated length is 3708 mm. Three different assemblies C2A, C3, and CB2 are used for the tests. C2A and C3 have radial power distribution similar to the one of the beginning of the reactor cycle, Fig. 3 A, while C2B has a radial power distribution characteristic of a middle of cycle, Fig. 3 B. In addition, C2A and C2B have a cosine axial power shape (Fig. 4 A) while C3 has an inlet-peak power shape, Fig. 4 B.
All different bundles have seven spacer grids distributed along the height with a pressure loss coefficient of 1.2. More details of the geometry are shown in Table 2.
The critical power tests were performed by step-wise increase of the bundle power. During this heat-up phase, the individual thermocouple signals of the heater rods measuring the local cladding temperature are monitored. The accuracy of the temperature measurements of the thermocouples is 1.5 • C and the accuracy of the power measurements is about 1.5%. The critical power was reached when the peak rod surface temperature became 14 • C higher than the one measured at steady-state test conditions. This criterion was chosen by the BFBT experimental team. Dry-out was observed at the rods with the peak power located at the peripheral row adjacent to the channel box. The boiling transition was always observed just upstream of the spacer grids. Fig. 5 shows the radial and axial thermocouple positions. Each thermocouple position is identified as follows: Rod Number -Axial location -Circumferential angle. Measurements are located in the rods 1, 2, 4, 5, 7, 8, 9, 12, 15, 16, 25, 31, 45, 52, 53, 54, 59, and 60, at different axial locations A (3521 mm), B (3009 mm), C (2497 mm), D (1985 mm), and different rotational angles, being sometimes more than one thermocouple for fuel rod. In TWOPORFLOW, it is not possible to calculate the circumferential CHF position on the rod surface. Due to the porous-media approach in TWOPORFLOW, such details cannot be resolved. For that reason, in this study the measured temperatures are compared with the local average temperatures predicted for the rods in the sub-channel, where CHF appears.
All tests for the assembly types C2A, C2B, and C3 were performed for a pressure of 7.2 MPa and the following boundary conditions. More details can be found in (Utsuno et al., 2006).

TWOPORFLOW model of BFBT test Bundle:
The different fuel BFBT-assemblies are modeled in TWOPORFLOW using a rod centered approach. It results in 8 × 8 X-Y-cell arrangement axially subdivided in 24 equidistant cells of 154.5 mm. The pins and water rods are characterized in the porous-media approach by a porosity of 0.57 and 0.17 in axial direction, respectively, Fig. 6. In radial direction, the porosity amounts 0.63. The heat transfer area-density (140.86 m − 1 ) is calculated from the outer radius of the fuel rods. The hydraulic diameters in axial direction (0.009 m corners, 0.016 m edges, 0.011 m inner sub-channels, 0.007 m water rods) and in lateral direction (0.010 m water rods and 0.021 m rest of sub-channels) are used to calculate local pressure loss and heat transfer. The heat conduction within the fuel rod is described by a radial 1-D heat conduction equation. The fuel pellet is subdivided into eight radial nodes. The cladding is represented by two radial nodes. No gap is presented in the fuel rod simulators.
All steady state BFBT critical power tests were analyzed by TWO-PORFLOW using the models described before.
In Figs. 7-9 the deviation (Simulated − Experimental CHF) of the CHF predicted by TWOPORFLOW from the measured data is plotted in dependence on the mass flux used as boundary condition for the different BFBT-assemblies C2A, C2B, and C3.
In Fig. 7 it can be observed that the Bowring correlation overpredicts the CHF for all tests. For mass fluxes lower than 600 kg/m 2 s the over-prediction is lower than 5.11%. For mass fluxes from 880 to 1930 kg/m 2 s the over-predictions increase up to a range of 4.96 to 11.85 %. CHF predicted by the Biasi and Groeneveld correlations show a similar trend. Both correlations under-predict the measured CHF for  In Fig. 8, it can be observed that the Biasi correlation over-predicts the CHF in all the cases with mass fluxes from 580 to 1930 kg/m 2 s in a range of 1.63 to 11.55% and for mass fluxes ~300 kg/m 2 s the differences are in a range of 0 to an under-prediction of 1.91%. Bowring over-predicts all the cases from 300 to 1320 kg/m 2 s in a range of 4.45 to 17.26%. The cases with mass flux ~1610 kg/m 2 s have differences in a range of ±2% and for higher mass fluxes the under-predictions are in a range of 0.48 to 6.26 %. Groeneveld LUT has a similar behavior as Bowring with over-prediction in all the cases from 580 to 1320 kg/m 2 s. For mass fluxes ~ 300 kg/m 2 s the differences are in a range of 0 to an under-prediction of 1.91%. The cases with mass flux ~1610 kg/m2s have differences in a range of ±0.3% and for higher mass fluxes the under-predictions are in a range of 0.48 to 5.85%. All three correlation show an increase in CHF difference from 290 to 890 kg/m 2 s and a decrease for higher mass fluxes.
In Fig. 9 it can be observed that all the methods over-predict the CHF from mass fluxes from 870 to 1910 kg/m 2 s, Biasi in the range of 1.42 to 0.11.04%, Bowring in the range of 7.65 to 12.23%, and Groeneveld in the range of 4.39 to 14.39%. For lower fluxes Biasi presents a range from 0 to a under-prediction of 7.38%, Bowring from 0 to a under-prediction of 1.24%, and Groeneveld a under-prediction from 3.46 to 7.38%.
Considering all tests, Table 3 gives the standard CHF-deviations of the three correlations (Biasi, Bowring, and Groeneveld) predicted for the three different assembly-types (C2A, C2B, C3).
Based on these results, it can be concluded that the three correlations are comparable, but the Biasi and Groeneveld are superior to the Bowring for the majority of the investigated fuel assembly types.

BFBT transient CHF tests
In (Utsuno et al., 2006) not only steady state but also two transient CHF-test were performed in the frame of the BFBT benchmark Phase II Exercise 3. For this purpose, two transient tests with boundary conditions representing BWR-transients such as the turbine trip without bypass and the recirculation pump trip were performed. These experimental data are used for the validation of the CHF-models implemented in TWOPORFLOW.

Tests boundary conditions
In the turbine trip without by-pass, the rapid closure of the turbine isolation valve leads to the propagation of a pressure wave from the main steam line to the core, and this in turn leads to the collapse of the void, a higher coolant density and, in consequence, the neutron moderation improves leading to a power increase. Key parameters of the test are the following: outlet pressure: 7.1 MPa, Inlet mass flow rate: 42 t/h, Inlet temperature: 276.7 • C for assembly type C2A and 275.5 • C for assembly type C3, and Bundle power: 8.5 MW. In Fig. 10, the change in time of the boundary conditions normalized to one are shown, which are used for the simulations with TWOPORFLOW.
In case of the recirculation pump trip, the mass flow rate is reduced after 22 s, increasing the void fraction and as consequence of it, the moderation decreases as well as the fission power. Afterwards, the mass Fig. 9. Deviation of the predicted CHF using the Biasi, Bowring, and Groeneveld LUT from the measured one depending on the mass flux for Assembly C2B.   Fig. 11. In (Utsuno et al., 2006) is indicated that the accuracy of the temperature and power measurements is similar to the ones of the steady state tests i.e., 1.5 • C and 1.5%, respectively.

BFBT model description in TWOFORFLOW
The assemblies C2A and C3 were discretized by TWOPORFLOW using the rod centered approach. As result, 8x8 cells exist in radial X-Yplane while in axial direction 48 cells of 77.25 mm each are considered. The additional parameters such as porosity, hydraulic diameter, heat transfer area-density, and pin radial nodes for the heat conduction are the same as the ones for simulation of the steady state tests.

Simulation of the BFBT turbine trip without bypass test
The transient CHF-test have been simulated with TWOPORFLOW using three correlations: Biasi, Bowring and Groeneveld. In the test, CHF is reached at 23 s, at the top of the spacers. In Table 4, the comparison of the CHF predicted by TWOPORFLOW with different correlations with the experimental data is shown for the two assembly types (CA2, C3). There, the time for appearance of CHF is also indicated.
The predictions of all correlations are close to the measured CHFvalue for both assembly types C2A and C3. In Fig. 12, the evolution of the cladding temperature in X-Y-cell 53 of assembly C3 (Fig. 5) predicted by TWOPORFLOW with three different CHF-correlations is compared with the test data. According to the test, the CHF appears in cell 53 at around 23 s of transient (see red circle in Fig. 12). There, the rapid increase of the cladding temperature due to the heat transfer deterioration (transition boiling, film boiling) can be observed. The gradient of the temperature escalation predicted by the three correlations is similar to the measured one. Nevertheless, the shape, and the time for the strong reduction of the cladding temperature differ from the measured ones. All the three correlations behave in similar manner. The final value of the cladding temperature is the same one of the test. This CHF-phenomenon are challenging for TWOFORFLOW and also for other codes such as COBRA-TF (Burns and Aumiller, 2007). TWOPORFLOW uses a simple interpolation for the transitions boiling curve between the critical heat flux point and the start of pure film boiling. In addition, there is no specific model for rewetting of the hot surface. The timing of the cladding temperature decrease is not described satisfactorily. In addition, TWOPORFLOW over-estimates the maximal cladding temperature by around 20 K.   Table 5 summarizes the comparison of the CHF measured in the two assemblies C2A and C3 with the CHF predicted by TWOPORFLOW for three CHF-correlations. It can be observed that the timing of appearance of CHF is reasonably well captured by the simulations. TWOPORFLOW always under-predicts the measured CHF using the three correlations. For the C2A assembly, the Biasi and Groeneveld correlations are closer to the data followed by Bowring. In case of the C3 assembly, the Biasi correlation predicts the CHF closer to the data, followed by Groeneveld and Bowring. Fig. 13 represents the comparison of the cladding temperature evolution of the cell 31 predicted by TWOPORFLOW with different correlations and measured one for the assembly type C3. The time of CHF appearance predicted by TWOPORFLOW with all three correlations is close to the measured one. In addition, the escalation of the cladding temperature shows a similar slope compared to the experiment. The peak height of the local temperature escalation is under-predicted by all three correlations. It can be state that the time-evolution predicted by TWOPORFLOW with the Biasi CHF correlation is closer to the data than the ones calculated with the Bowring and Groeneveld correlations.

Simulation of the BFBT recirculation pump trip test
Considering the predicted CHF by other codes (Ferrouk et al., 2008;Kim et al., 2018;Choi et al., 2012;Glück, 2007) and their comparison with the measured data, it can be stated that TWOPORFLOW is capable to perform good predictions of the CHF measured by two-BWR relevant tests.

Conclusions and outlook
The prediction capability of TWOPORFLOW using three CHFcorrelations is validated based on the Becker experiments and the BFBT CHF data provided by the benchmark.
For the Becker CHF steady state experiments (tubes) TWOPORFLOW shows the larger deviations using the Biasi correlation 9.19%. The Bowring correlation, on the other hand, presented the most accurate results with deviations of 1.30%. Groeneveld LUT shows similar results than the Bowring correlation with deviation of 1.43%. In the case of the BFBT CHF steady state tests, for assembly C2A the correlations present standard deviations of less than 2.5%. In the case of assembly C2B and C3 the averaged CHF deviations are less than 5%. In the case of BFBT transient tests, the deviation from the measured data observed for all used correlations is less than 2 % for assembly C2A and less than 9% for assembly C3. The axial appearance is inside the given deviation in all the cases. In the case of access to more actual CHF experimental data, a broader validation process will be performed.
Based on the encouraging validation work, the TWOPORFLOW code will be applied to simulate both stationary and accidental conditions of BWR-core representing each fuel assembly as a computational node in radial direction. As next steps, the coupling of TWOPORFLOW with diffusion and transport neutronics solvers such as PARCS (Downar et al., 2018) and PARAFISH (Van Criekingen et al., 2011) is foreseen for the safety related investigations.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.  Fig. 13. Comparison of evolution of the local (i.e., 2497 mm elevation) cladding temperature of the rod located at cell 31 (Fig. 5) predicted by TWOPORFLOW with different CHF-correlations in the BFBT benchmark for assembly type C3.