Abstract
In the framework of the OECD/NRC PSBT benchmark, the subchannel grade void distribution data and DNB data were assessed by a subchannel code, MATRA. The prediction accuracy and uncertainty of the zone-averaged void fraction at the central region of the 5 × 5 test bundle were evaluated for the steady-state and transient benchmark data. Optimum values of the turbulent mixing parameter were evaluated for the subchannel exit temperature distribution benchmark. The influence of the mixing vanes on the subchannel flow distribution was investigated through a CFD analysis. In addition, a regionwise turbulent mixing model was examined to account for the nonhomogeneous mixing characteristics caused by the vane effect. The steady-state DNB benchmark data with uniform and nonuniform axial power shapes were evaluated by employing various DNB prediction models: EPRI bundle CHF correlation, AECL-IPPE 1995 CHF lookup table, and representative mechanistic DNB models such as a sublayer dryout model and a bubble crowding model. The DNBR prediction uncertainties for various DNB models were evaluated from a Monte-Carlo simulation for a selected steady-state condition.
1. Introduction
The critical heat flux (CHF) is a parameter of great importance, which constrains the thermal power capability of a light water nuclear reactor (LWR). It is usually predicted by a local parameter CHF correlation accompanied with an appropriate thermal-hydraulic field analysis code to obtain the local subchannel grade conditions in the fuel assembly. For this purpose, the subchannel approach has been widely adopted in the design calculation of an LWR core since it provides reasonably accurate results on the flow and enthalpy distributions in rod bundles with a pertinent computing time.
The OECD/NRC PWR Subchannel and Bundle Tests (PSBT) benchmark was organized on the basis of the NUPEC database. The purposes of the benchmark are the encouragement to develop a theoretically based microscopic approach as well as a comparison of currently available computational approaches. The benchmark consists of two separate phases: a void distribution benchmark and DNB benchmark. Subchannel-grade void distribution data was employed for validation of a subchannel analysis code under steady-state and transient conditions. The DNB benchmark provided subchannel fluid temperature data, which can be used to determine the turbulent mixing parameter for a subchannel code. The steady-state and transient DNB data can be used to evaluate and improve the currently available DNB prediction models in PWR bundles. The NUPEC PWR test facility consists of a high-pressure and high-temperature recirculation loop, a cooling loop, and a data recording system [1]. The void fraction was measured by two different methods: A gamma-ray beam CT scanner system was used to determine the distribution of density/void fraction over the subchannel at a steady-state flow and to define the subchannel averaged void fraction with an accuracy of ±3%. A multibeam system was also used to measure chordal averaged subchannel void fraction in a rod bundle with accuracies of ±4% and ±5% for steady and transient states, respectively.
The purposes of this study are to evaluate the accuracy and uncertainty of the MATRA code for the PSBT benchmark exercises of the subchannel void distribution under steady-state and transient conditions, as well as to provide the analysis results for the subchannel temperature distribution and DNB benchmark problems. The prediction accuracy, which is the closeness of agreement between the predicted value and the true value, was evaluated by comparing the code predictions with the measured data. On the other hand, the prediction uncertainty was evaluated by considering the propagation of uncertainties for the boundary operating conditions and code modelling parameters.
2. MATRA Code Models for PSBT Benchmark Analysis
2.1. Subchannel Analysis Code MATRA
MATRA is a subchannel analysis code, which adopts mixture transport equations for the two-phase flow conditions [2]. The governing equations for the subchannel geometry were derived from the integral balances on an arbitrary fixed control volume. It was assumed that the flow is a transient, single component, and a two-phase mixture of liquid and vapor in a thermodynamic equilibrium state. The continuity, energy, and axial/lateral momentum equations for a vertical subchannel are expressed as follows.
Continuity:
Energy:
Axial momentum:
Lateral momentum: The source terms in the above equations are expressed as
The major unknowns of the governing equations are the coolant density (), axial flow rate (), cross flow (), pressure (), and enthalpy (). Subscript “” means the two-phase mixture property, and superscript “*” means the donor-channel property. The last term of the left-hand side of conservation equations (1), (2), and (3) represents the net exchange of mass, energy, and axial momentum due to a turbulent mixing between subchannel and its surrounding subchannels. By introducing a turbulent mixing parameter,, which is defined as a ratio of the lateral fluctuating mass flux to the axial mass flux of the fluid in the subchannel, the turbulent mixing flow rate per unit length from subchannel to is expressed as
The turbulent mixing parameter is normally determined from a thermal mixing test under single-phase conditions. Two different models for the interchannel turbulent mixing phenomenon are available in the MATRA code: the equal-mass-exchange (EM) model and the equal-volume-exchange and void-drift (EVVD) model. The net fluctuating mass velocity from channel to channel for the EVVD model is expressed as where is the turbulent mixing flow rate per unit length under single-phase conditions and is a two-phase multiplier for the turbulent mixing rate. A detailed description of the correlation parameters is given in [3]. The net fluctuating mass velocity () becomes zero when the EM model is employed for the analysis. Important models of the MATRA code for an analysis of the PSBT benchmark are summarized in Table 1.
2.2. PSBT Benchmark Test Bundles
Geometry information for the PSBT benchmark test bundles is summarized in Table 2. The benchmark data provided the steady-state void fraction averaged over the four central subchannels (CNTR) as shown in Figure 1. The region averaged void fraction was evaluated from the void fractions of subchannels 15, 16, 21, and 22 calculated by the MATRA code.
3. Analysis of Void Distribution Benchmark Problems
3.1. Single-Channel Void Distribution Benchmark
A single subchannel void distribution benchmark provides cross-sectional averaged void fraction at the exit of four different subchannel types found in a PWR assembly: central typical, central thimble, side, and corner subchannel types, which correspond to channel types I, II, III, and IV as shown in Figure 1. The homogeneous equilibrium bulk boiling model tends to overpredict the void fraction in comparison with the modified Armand model, while the low-void fraction was not adequately predicted if a subcooled boiling model was not employed. As the result, this model revealed relatively large scatterings against the void fraction measurement error () [1] as shown in Figure 2(a). On the other hand, the effects of subcooled boiling models in the MATRA code were examined for the single-channel void distribution benchmark exercises. The subcooled and bulk boiling models described in Table 1 revealed a slightly overprediction of the channel exit void fraction for central typical and side subchannels, as shown in Figure 2(b).
(a) Homogeneous equilibrium model
(b) Subcooled boiling model
3.2. Steady-State Bundle Void Distribution Benchmark
The experimental data include a chordal averaged void fraction at three axial elevations, 2216 mm, 2669 mm, and 3177 mm, from the bottom of the heated length. The void prediction error, which represent the difference between the predicted and measured void fractions, are plotted in Figure 3 at the lower and upper axial levels. As the axial elevation increases, the error of mean (P-M) decreases from 5.1% to −3.4%, while the standard deviation remains within 6.0%~6.9% for all axial levels. The maximum error of mean (P-M) was calculated as 11% at the lower elevation of B7, which has a central unheated rod and cosine axial power shape. The statistics of the void prediction error are summarized in Table 3.
(a)
(b)
An uncertainty of the predicted void fraction may be caused by various parameters involved in the calculations, such as the hydraulic and thermal input variables, dimensions of the bundle components, and modelling parameters used in the code. A prediction uncertainty of the void fraction was evaluated for a case selected from the benchmark exercises. The pressure, inlet temperature, mass flux, and bundle power for the selected case (Run no. 7.6322) were 4.87 MPa, 168.6°C, 2239 kg/m2s, and 3536 kW, respectively. The selected parameters for this uncertainty analysis are described in Table 4. The accuracy and probability distribution function for the operating parameters were determined from the benchmark specifications [1]. In this table, the accuracy ranges represented one standard deviation for the normal distribution, or the boundary values for the flat distribution. The uncertainty data for the modelling parameters, which includes the bundle pressure loss factors and the turbulent mixing parameter, were approximately determined from the design information for a typical PWR [4]. The uncertainties of the selected parameters were propagated to the uncertainty of the void fraction through the physical models of two-phase flows employed in the MATRA code.
The uncertainty of a void fraction was evaluated from a Monte-Carlo simulation with the MATRA code by employing a simple random sampling technique for the selected parameters. The frequency distributions of the void fractions calculated by the MATRA code at lower and upper elevations are plotted in Figure 4 from 2,000 samplings combining each parameter. As shown in the figure, it was revealed that the uncertainty of the calculated void fraction tends to reduce as the axial elevation increases.
(a)
(b)
The two-sided tolerance limits with 95% probability and a 95% confidence level were evaluated for estimating the magnitude of uncertainty for this selected case, where and s are the sample mean and sample standard deviation of the predicted void fraction, respectively. A nonparametric method [5] was applied to determine the tolerance limits, since it was revealed that the data was not drawn from a normally distributed population at a 5% significance level. As a result, the interval of tolerance limits for the region averaged void fraction at the axial levels of 2216 mm, 2669 mm, and 3177 mm were calculated as 4.6%, 2.3%, and 1.5%, respectively.
3.3. Transient Bundle Void Distribution Benchmark
A transient bundle void distribution benchmark provided subchannel averaged void fraction at 3 different axial levels of the CNTR region under four transient conditions: power increase (PI), flow reduction (FR), depressurization (DP), and temperature rise (TI). These data are important for a benchmark of the subchannel analysis codes in terms of predicting a CHF for a reactor transient or accident conditions. The implicit scheme employed in the MATRA code allowed no restrictions in the time step size for the transient calculations. The analysis results revealed the maximum deviation between the predicted and measured values of the averaged void fraction for the four different transients of PI, FR, DP, and TI as 0.247, 0.216, 0.209, and 0.285, respectively. For the FR transient, the maximum/minimum values of (P-M) of the void fraction at CNTR region were calculated as 0.076/−0.126 for B5 and 0.216/−0.077 for B7, respectively. Similar to the steady-state bundle void distribution benchmark, the maximum deviation between the benchmark data and the MATRA prediction was found at the lower elevation of the B7 bundle for this transient problem.
An uncertainty analysis for a transient bundle void distribution benchmark was conducted for the parameters described in Table 4. The FR transients for test bundles B5 and B7 were selected for the uncertainty analysis. From a Monte-Carlo simulation with the MATRA code, the range of the predicted void fraction at each time step was plotted as the grey region in Figure 5. The maximum deviations of the void fraction between the maximum and minimum predicted values (i.e., the width of grey region) at each time step were calculated as 11% for B5 and 12%~13% for B7 at different axial levels during the FR transients.
(a)
(b)
(c)
(d)
4. Analysis of Mixing and DNB Benchmark Problems
4.1. Subchannel Exit Temperature Distribution Benchmark
The accuracy of a subchannel analysis code is fairly dependent on the modeling of interchannel exchanges between adjacent subchannels, such as a diversion cross flow and turbulent mixing [6]. The diversion cross flow caused by a lateral pressure difference between subchannels is calculated from the lateral momentum equation. The turbulent mixing term is considered in the subchannel equations by (6), which implies a turbulent mixing parameter defined as where , , and are the mean turbulent eddy diffusivity of the heat, the mixing distance, and the fluctuating velocity between subchannels and , respectively. Traditionally, a constant value of has been applied to the whole subchannel analysis region.
In the PSBT subchannel exit temperature distribution benchmark, the subchannel exit temperatures were measured at the exit of a 5 × 5 mixing test bundle (A1), which has a large gradient of radial power distribution. The average power in the hot region is 4-times higher than that in the cold region. The optimum value of is determined at the condition when the rms error of the subchannel exit temperatures between the predicted and measured values becomes minimum for various values of . An increase of in the MATRA code will reduce the predicted difference of the temperature rise between the hot and the cold regions. Figure 6 shows a comparison of the subchannel exit temperature distributions between the MATRA results and PSBT data. If there is a temperature gradient in each region due to additional mixing effect, it may reduce the exit temperature differences between the two regions, and finally increase the optimum in comparison with the case of no temperature gradient in each region. As shown in the Figure 6, a temperature gradient between HH-HC, as well as CH-CC, was observed in the experimental data which did not appear in the MATRA results. This discrepancy results in a significantly higher value of the optimum , as shown in Table 5. For various operating conditions, the mean value of was calculated as 0.08 as summarized in Table 5.
It was thought that the temperature gradient appearing in the experimental data was attributed to the thermal mixing in the diagonal direction of the test bundle, which may be caused by the alignment of mixing vanes mounted in the spacer grids. A numerical analysis using a CFD code (ANSYS, version 12.1) was conducted for a selected case (Test number 01-5237) to investigate the influence of spacer grids with mixing vanes on the subchannel temperature gradient in test bundle A1. A simplified grid model was adopted for the analysis which focused on the effect of the mixing vanes. The computational meshes for the upper and lower parts of the mixing vanes were generated by applying a hexagonal sweep method, while tetrahedron meshes were applied to the mixing vane region. The standard k-ε turbulence model and the convergence criterion of the 10−6 RMS residual were used for the analysis of the flow characteristics. The geometry of the mixing vanes remarkably affects the flow distribution, as shown in Figure 7.
Figure 7 shows the streamline from the inlet. To investigate the mixing between each region, the start surface and expressed colour of the streamline were divided into two surfaces and colours. The blue and red colours indicate the streamline starting at the HC and HH regions, respectively. The cross-sectional view is configured from the end of the first mixing vane. It was noted that the flow direction is drastically changed after going through the mixing vane, and the cross flow is formed as shown in this Figure 7. Moreover, the cross flow generated by the mixing vane induces vigorous mixing in the diagonal direction between the HC and CH regions.
The diagonal mixing effects caused by the mixing vanes cannot be reflected adequately by a single optimum in the subchannel analysis code. From the experimental result, it was inferred that the thermal mixing inside the HC region may occur more vigorously than the region HH as shown in Figure 6. A regionwise model was introduced into the MATRA code to examine the local mixing effect. For the selected case of 01-5237, the regionwise values were empirically fitted from the subchannel-wise exit temperature distribution data, which results 0.015, 0.4, 0.2, and 0.03 for the HH, CH, HC, and CC regions, respectively. As shown in Figure 8, the regionwise mixing model may improve the accuracy of the subchannel exit temperatures, specifically for the rod bundles with nonhomogeneous mixing characteristics due to mixing vanes. However, it is necessary to develop an appropriate algorithm for determining nonhomogeneous turbulent mixing parameters from the experimental data and/or CFD analysis results.
(a)
(b)
4.2. Steady-State DNB Benchmark
The steady-state DNB benchmark data provided the power at which DNB occurred and the corresponding location in the bundle. The prediction accuracy of the MATRA code was evaluated for uniform and nonuniform axial power shapes by employing various DNB prediction models: a generalized bundle CHF correlation, a generalized tube CHF correlation, and two kinds of mechanistic DNB models.
The EPRI correlation [7] is a generalized bundle CHF correlation that has been developed on the basis of the local fluid conditions obtained with a subchannel code. The applicable range of the correlation covers normal operating conditions of PWR and BWR as well as hypothetical LOCA conditions.
A CHF lookup table method [8] was adopted as a generalized tube CHF correlation. It provides CHF values for water-cooled tubes at discrete values of pressure, mass velocity, and critical quality. Linear interpolation between table values gives the CHF for a specific condition, and several correction factors were introduced to extend the CHF table to various shapes of the boiling channel. The applicability of the CHF lookup table method has been assessed for rod bundle CHF data with a subchannel analysis code [9].
Two mechanistic models for DNB were selected: a sublayer dryout model and a bubble crowding model. A sublayer dryout model was developed on the basis of dryout mechanisms of a thin liquid sublayer underneath a vapor blanket (microlayer) resulting from the Helmholtz instability at the microlayer-vapor interface. Lin and his coworkers [10] suggested an improved model for low quality flows. Key parameters associated with this model are sublayer thickness, blanket velocity, and blanket length. The sublayer thickness was determined by a force balance on the vapor blanket in the radial direction between the evaporation thrust and the lift force. The velocity of vapor blanket was determined by a balance between the buoyancy force and drag force. The mean length of a vapor blanket was assumed to be equal to the critical Helmholtz wavelength at the vapor-liquid interface. They examined the model for selected rod bundle data with the COBRA-IIIC/MIT-1 code [11].
The near-wall bubble crowding model postulated that CHF at low qualities occurs when the bubbles near the heated wall coalesce into a vapor film. The bubble layer near the wall was assumed to become so thick that it inhibits enthalpy transport between the fluid in the core region and the liquid near the wall. Weisman and his coworkers suggested that CHF occurs when the void fraction in the bubbly layer just exceeds the critical value of 0.82. During the modelling of a turbulent interchange between the bubbly layer and the core regions, two empirically determined parameters have been adopted. Weisman and Ying [12] suggested an extended version of the bubble crowding model to high quality and low velocity conditions. They evaluated the model for rod bundle data under PWR conditions using the COBRA-IV-I code [13].
It is known that two different methodologies can be used for the prediction of CHF by a local parameter CHF correlation: the direct substitution method (DSM) and the heat balance method (HBM) [14]. For a calculation of CHF in rod bundles by HBM, it was assumed that (1) the local mass velocity distribution within a test bundle remained unchanged at different power levels, and (2) the enthalpy rise in the subchannel was proportional to the bundle power. The EPRI correlation was applied with DSM, while the other three CHF prediction models were applied with HBM. At a given experimental condition, the CHF values were calculated in all subchannels with the local thermal hydraulic parameters calculated by the MATRA code to find out the minimum DNBR location in the test bundle.
For the test bundles with nonuniform axial power shapes, the predicted CHF value was determined at the location where the critical bundle power becomes minimum. An iterative calculation was required for the HBM, until the predicted critical bundle power coincided with the bundle power, which was used for a determination of the local thermal hydraulic conditions. During this iterative calculation, it was assumed that the ratio of the hot-subchannel enthalpy rise to the bundle-averaged enthalpy rise remains constant for various bundle powers. Pertinent correction factors for a nonuniform axial power shape were employed for the EPRI correlation and CHF lookup table method, while an additional correction factor was not required for the mechanistic DNB models.
The statistics of predicted-to-measured CHF ratio (P/M) for the steady-state DNB benchmark exercises are summarized in Table 6. The mean relative bias error () is defined as follows:
The parametric behaviours of P/M with respect to the bundle average mass flux are shown in Figure 9 for different DNB models. As a result, it was found that the EPRI correlation and CHF lookup table revealed similar prediction accuracy. The sublayer dryout model underpredicted the CHF by 4.9%, while the bubble crowding model overpredicted it by 16.4%. The standard deviation for all of the benchmark data was calculated to be within 7.6%~8.9% for the four different DNB prediction models.
(a)
(b)
(c)
(d)
The uncertainty of predicted DNBR was evaluated by considering the uncertainty parameters described in Table 4. The frequency distributions of DNBR at a selected condition (test bundle A0, Run 5350) from a Monte-Carlo simulation with 2,000 samplings are shown in Figure 10 for various DNB models. As a result, it was appeared that the DNBR frequency distribution for a bubble crowding model was somewhat different from those for the empirical DNB models as shown in the Figure 10. It may be attributed to the uncertainty of void fraction which propagated to the uncertainty of mechanistic DNB model through the modelling parameters, such as the effective quality or the turbulent transport velocity between core and bubbly layer. The relative standard deviations of DNBR for the four different DNB prediction models lay in 1.0%~1.2%, while the relative standard deviation of the local void fraction was calculated as 5.5%.
(a)
(b)
(c)
(d)
5. Conclusions
The PSBT benchmark data for the void fraction and DNB exercises were evaluated by a subchannel analysis code, MATRA. The prediction accuracy of the MATRA code for the steady-state bundle void distribution exercises revealed a mean (P-M) error of 2.9%~5.1%. The MATRA code tended to over-predict the void fraction at lower elevations under the steady and transient conditions.
The uncertainty of the MATRA code under steady-state and transient conditions was evaluated by considering the uncertainties of selected operating and modelling parameters with a Monte-Carlo simulation method. For Run no. 7.6322 of the steady-state bundle void distribution benchmark, the interval of tolerance limits for the region averaged void fraction at the axial levels of 2216 mm, 2669 mm, and 3177 mm were calculated as 4.6%, 2.3%, and 1.5%, respectively. For the flow reduction transient, the maximum/minimum values of (P-M) of the void fraction at CNTR region were calculated as 0.076/−0.126 for B5 and 0.216/−0.077 for B7, respectively.
For the fluid temperature benchmark data, the best estimated value of the turbulent mixing parameter () was calculated as 0.08 with a standard deviation of 0.011. As a result of the CFD analysis, the temperature gradients inside the hot and cold regions were possibly explained by the mixing vane effect. A regionwise model was examined to reflect the diagonal mixing effects caused by mixing vanes. It was revealed that the temperature gradients in the hot and cold regions could be predicted approximately by using the regionwise model.
DNB data for test bundles with uniform and cosine axial power shape was evaluated by employing various CHF prediction models. The mean error and standard deviation of P/M for the two different mechanistic DNB models, that is, a sublayer dryout model and a near wall bubble crowding model, were calculated as −4.9%~+16.4% and 8.3%~8.9%, respectively. The AECL-IPPE 1995 CHF lookup table with HBM and EPRI CHF correlation revealed the mean error to be −3.0% and −2.7%, respectively. As a result of DNBR uncertainty analysis for a selected case (test bundle A0, Run 5350), the relative standard deviations were calculated as 1.0%~1.2% with various DNB models.
Further research works for improving the subcooled boiling models in rod bundles, or for an accurate modeling of the spacer grid effects on the thermal-hydraulic interactions between subchannels in connection with a CFD analysis would contribute to enhance the accuracy of the subchannel codes.
Nomenclature
| : | Channel flow area (m) |
| : | Channel hydraulic diameter (m) |
| : | Single-phase friction factor |
| : | Turbulent momentum factor |
| : | Gravity acceleration (m/sec) |
| : | Mass flux (kg/m-sec) |
| : | Flow enthalpy (kJ/kg) |
| : | Mixture enthalpy (kJ/kg) |
| : | Thermal conductivity of fluid (kW/m-) |
| : | Two-sided tolerance factor |
| : | Form loss factor |
| : | Crossflow resistance factor |
| : | Void drift coefficient |
| : | Centroid distance between adjacent channels (m) |
| : | Mass flow rate (kg/sec) |
| : | Pressure (kg/m-sec) |
| : | Surface heat flux (kW/m) |
| : | Sample standard deviation |
| : | Gap between channel and (m) |
| : | Fluid temperature () |
| : | Axial flow velocity (m/sec) |
| : | Lateral flow velocity (m/sec) |
| : | Effective specific volume (m/kg) |
| : | Crossflow from channel to (kg/m-sec) |
| : | Turbulent mixing flow rate per unit axial length (kg/m-sec) |
| : | Sample mean |
| : | Void fraction |
| : | Turbulent mixing parameter |
| : | Mean relative bias error |
| : | Two-phase friction multiplier |
| : | Two-phase mixture density (kg/m) |
| Standard deviation | |
| : | Heated perimeter (m). |
| avg: | Averaged for channel and |
| : | Channel index |
| : | Rod index. |
| : | Donor channel property. |
Acknowledgments
The authors thank Dr. Utsuno at NUPEC, Professor Avramova at Pennslvania State University and Dr. Choi at OECD for their efforts pertaining to the benchmark program. This work has been carried out under the R&D program supported by the Ministry of Education, Science and Technology of the Republic of Korea.