# CFD Analysis of a Void Distribution Benchmark of the NUPEC PSBT Tests: Model Calibration and Influence of Turbulence Modelling

Authors: E. Krepper; R. Rzehak
Journal: Science and Technology of Nuclear Installations
 (2012)

---

## Abstract

The paper presents CFD calculations of the void distribution tests of the PSBT benchmark using ANSYS CFX-12.1. First, relevant aspects of the implemented wall boiling model are reviewed highlighting the uncertainties in several model parameters. It is then shown that the measured cross-sectionally averaged values can be reproduced well with a single set of calibrated model parameters for different test cases. For the reproduction of patterns of void distribution cross-sections, attention has to be focussed on the modelling of turbulence in the narrow channel. Only a turbulence model with the capability to resolve turbulent secondary flows is able to reproduce at least qualitatively the observed void distribution patterns.


---

## Body

## 1. Introduction

Based on NUPEC PWR Subchannel and Bundle Tests (PSBT), an international benchmark has been promoted by OECD and NRC and coordinated by Penn State University (PSU). The benchmark includes void distribution and departure from nucleate boiling exercises. In the first exercises the void fraction distribution was investigated in a steady state subchannel grade benchmark. In the present paper some of the single-channel steady state void fraction measurements are analysed to investigate the capabilities of present CFD modelling of wall boiling.For engineering calculations, currently the most widely used CFD approach to model two-phase flows with significant volume fractions of both phases is the Eulerian two-fluid framework of interpenetrating continua (see, e.g., [1–3]). In this approach, balance equations for mass, momentum, and energy are written for each phase, that is, gas and liquid, separately and weighted by the so-called volume fraction which represents the ensemble averaged probability of occurrence for each phase at a certain point in time and space. Exchange terms between the phases appear as source/sink terms in the balance equations. These exchange terms consist of analytical or empirical correlations, expressing the interfacial forces, as well as the heat and mass fluxes, as functions of the average flow parameters. Since most of these correlations are highly problem-specific, their range of validity has to be carefully considered and the entire model has to be validated against experiments.For the case of boiling flows, where heat is transferred into the fluid from a heated wall at such high rates that vapour is generated, additional source terms describing the physics of these processes at the heated wall have to be included. A CFD wall boiling model implemented in CFX following the lines of Kurul and Podowski [4, 5] was calibrated and validated by several authors, for example, Krepper et al. [6] against experimental results of Bartolomej and Chanturiya [7]. In these tests, subcooled flow boiling of water at high pressure flowing upwards in a vertical pipe heated from the outside was investigated and measurements of the axial development of void fraction, wall temperature and cross sectionally averaged liquid temperature were provided.The aim of this work is to investigate the applicability of the CFX models to the PSBT tests. In Section1 the PSBT tests for the void fraction distribution and the parameters of the selected tests are described briefly. In Section 2 the most important model details used for the simulations are described and critically reviewed. In Section 3 the sensitivity of the different parts of the models is investigated in view of the knowledge discussed in Section 2. The parameters are calibrated to only a single test. In the following Section 4 these calibrated models are applied to additional tests without any further change of model parameters. The cross sectionally averaged values of the selected tests are summarized for comparison to the calculations of other participants in the benchmark. In Section 5 the influence of the turbulence model on the void fraction distribution is investigated.


## 2. The OECD Benchmark Test

The PWR Subchannel and Bundle Tests (PSBT) were conducted by NUPEC (1987–1993) within an extensive experimental campaign aimed at verifying the reliability of fuel assemblies used for commercial nuclear power plants [8, 9]. Void fraction measurements and departure from nucleate boiling (DNB) tests were performed under PWR thermal-hydraulic conditions including steady states and transients such as power increase, flow reduction, depressurization, and temperature increase. The void fraction in each experiment was measured by gamma-ray transmission. These tests form the basis of an OECD benchmark for CFD and subchannel codes [10].The subchannel test section, as shown in Figure1, simulates a single subchannel of a PWR fuel assembly. The effective heated length is 1500 mm where the void measurement section is located near close to the top end at 1400 mm from the bottom of the heated section.Figure 1
Test section for the subchannel void measurement [10].For the analysis different tests in different pressure regions were selected (see Table1).Table 1
Selected tests.
Run
Pressure [MPa]
Mass flux [106 kg m−2 h−1]
Inlet temperature [°C]
Power [kW]
Outlet quality
1.2211
15
11
295.4
90
−0.04
1.2223
15
11
319.6
70
0.04
1.2237
15
11
329.6
60
0.08
1.2422
15
5
284.1
60
0.02
1.2423
15
5
299.3
60
0.10
1.4324
10
5
238.9
60
0.00
1.4325
10
5
253.8
60
0.05
1.4326
10
5
268.8
60
0.11
1.4411
10
2
253.7
20
0.00
1.6222
5
5
204.2
50
0.02


## 3. The Models

### 3.1. Modelling of Boiling at a Heated Wall: The General Model Structure

In boiling, heat is transported from the hot wall to the fluid by several different mechanisms. On parts of the wall, where no bubbles reside, heat flows directly to the subcooled liquid in the same way as in single-phase flow.  On parts of the wall where bubbles grow, heat is consumed by the generation of vapour which occurs at the so-called nucleation sites. Moreover, there is a liquid mixing mechanism due to the bubbles which leave the wall. As a consequence of the recirculation around the detaching bubbles, cold liquid from the bulk of the flow is brought into contact with the hot wall which leads to additional cooling. This mechanism is termed quenching. The total supplied heat flux is accordingly expressed as the sum of three contributions as(1)Qtot=QC+QQ+QE,
where QC, QQ, and QE denote the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively. The individual components in this heat flux partitioning are then modelled as functions of the wall temperature and other local flow parameters. Once this is accomplished, (1) can be solved iteratively for the local wall temperature TW, which satisfies the wall heat flux balance. Denoting the fraction of area influenced by the bubbles as AW, the heat flux components are expressed as discussed in the following.The turbulent convection heat flux is calculated in the CFX model version in much the same way as for a pure liquid flow without boiling, but multiplied by the fraction of area unaffected by the bubbles, that is,(2)QC=(1-AW)hC(TW-TL).
Here, hC is the heat transfer coefficient which is written using the temperature wall function T+(y+) known from Kader [11] as
(3)hC=ρCPuτT+,
where nondimensional variables (indicated by superscript “+”) and the friction velocity uτ are defined as usual. Note that (2) and (3) may be evaluated at any location y, provided it is used consistently whenever a variable depends on position. QQ is represented in terms of the quenching heat transfer coefficient hQ:
(4)QQ=AWhQ(TW-TL).
A grid-independent solution for QQ is obtained by evaluating the non-dimensional temperature profile at a fixed value of y+. The evaporation heat flux QE is obtained via the evaporation mass flux at the wall:
(5)QE=m˙WHLG,
where the generated vapour mass  m˙W  is expressed in terms of the bubble diameter at detachment dW, bubble generation frequency f, and nucleation site density  N  as
(6)m˙W=ρGπ6dW3fN.
In terms of bubble detachment diameter and nucleation site density, the wall area fraction AW influenced by vapour bubbles is given by
(7)AW=π(adW2)2N.
Here, a is the so-called bubble influence factor, for which a value of 2 is commonly used [4, 5]. Since AW=1 corresponds to the case where the whole surface is under the influence of bubbles, AW as calculated by (7) has to be limited to values smaller than one. Moreover, it should be kept in mind that already as AW approaches 1 the assumptions of the model are not really satisfied anymore.Correlations for the yet undetermined quantities used in the CFX wall boiling model are discussed in the following.


### 3.2. Bubble Size at Detachment

The bubble size at detachment depends on the liquid subcooling. Also the liquid properties which depend on the system pressure, the flow rate, and the heat flux have an influence.An investigation of the bubble size at detachment was performed by Tolubinsky and Kostanchuk [12] for water at different pressures and subcoolings. The observed dependence on the liquid subcooling at atmospheric pressure can be fitted to a correlation:(8)dW=drefe-((Tsat-TL)/(ΔTrefd)).
A least squares fit gives parameter values are dref=0.0013 m and ΔTrefd=53 K under these conditions.To match the tests of Bartolomej and Chanturiya [7] which were conducted at much higher pressures relevant for typical nuclear energy applications the values of dref and ΔTrefd had to be adjusted to dref=0.6 mm and Trefd=45 K (e.g., [6]).


### 3.3. Nucleation Site Density

The situation concerning data on nucleation site density is much less clear. Most of the time, correlations are expressed in the form of power laws depending on the wall superheat as(9)N=Nref(TW-TsatΔTrefN)p.
However, a recent compilation [13] shows that vastly different parameter values are required to match different data sets. A likely reason for this fact is that nucleation site density is highly dependent on the microscale topography of the boiling surface, which in turn depends strongly on the processes that were used to finish the surface. These processes are very diverse and in most boiling experiments not specifically controlled or characterized.For the tests of  Bartolomej and Chanturiya [7] Nref=0.8×106m-2 and ΔTrefN=10 K were found to yield the best results in the model framework.


### 3.4. Further Model Aspects

The bubble detachment frequencyf  is given according to Cole [14] as a function of the detachment size dW. The quenching heat transfer coefficient is calculated as suggested by Mikic and Rohsenow [15].In the bulk, vapour is assumed to be at saturation condition. Where the liquid is subcooled, that is,TL<Tsat, vapour is condensing and calculated with the transfer rate:
(10)m˙=hLG(Tsat-TL)AIHLG.
Here, hLG is the interfacial heat transfer coefficient, calculated according to Ranz and Marshall [16] and AI is the interfacial area. To close the phase transition model in the bulk bubbly flow with a mean bubble diameter dB, Kurul and Podowski [4] and also Anglart et al. [17] proposed to calculate the bubble diameter dB locally as a linear function of liquid subcooling Tsub:
(11)dB=dB,1(Tsub-Tsub,2)+dB,2(Tsub-Tsub,1)Tsub,2-Tsub,1.
Reference subcooling conditions for typical nuclear energy applications have been given as dB,1=0.1 mm at Tsub,1=-13.5 K and dB,2=2 mm at  Tsub,2=5 K in [17].In a two-phase turbulent flow gaseous bubbles have an influence on the liquid turbulence. This effect is described by increasing the turbulent viscosity according to Sato et al. [18]. Conversely the liquid turbulence flow structures influence the gas distribution, which is described by a turbulent dispersion force according to Burns et al. [19].For momentum exchange between the phases, finally, the Ishii and Zuber [20] drag law was used. Furthermore, a lift force according to Tomiyama et al. [21], a turbulent dispersion force according to Burns et al. [19], and a wall force according to Antal et al. [22] were included.


### 3.5. Evaluated Quantities

Quantities requested for the benchmark have been evaluated as follows. The mixed densityρMIX was calculated according to
(12)ρMIX=αLρL+αGρG,
and the mixed enthalpy HMIX according to
(13)HMIX=αLρLHL+αGρGHGρMIX
with the temperature-dependent static enthalpies for liquid HL and gas HG. The quality X was then calculated according to
(14)X=HMIX-HLSatHGSat-HLSat
with the enthalpies at saturation for liquid HLSat and gas HGSat.


## 3.1. Modelling of Boiling at a Heated Wall: The General Model Structure

In boiling, heat is transported from the hot wall to the fluid by several different mechanisms. On parts of the wall, where no bubbles reside, heat flows directly to the subcooled liquid in the same way as in single-phase flow.  On parts of the wall where bubbles grow, heat is consumed by the generation of vapour which occurs at the so-called nucleation sites. Moreover, there is a liquid mixing mechanism due to the bubbles which leave the wall. As a consequence of the recirculation around the detaching bubbles, cold liquid from the bulk of the flow is brought into contact with the hot wall which leads to additional cooling. This mechanism is termed quenching. The total supplied heat flux is accordingly expressed as the sum of three contributions as(1)Qtot=QC+QQ+QE,
where QC, QQ, and QE denote the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively. The individual components in this heat flux partitioning are then modelled as functions of the wall temperature and other local flow parameters. Once this is accomplished, (1) can be solved iteratively for the local wall temperature TW, which satisfies the wall heat flux balance. Denoting the fraction of area influenced by the bubbles as AW, the heat flux components are expressed as discussed in the following.The turbulent convection heat flux is calculated in the CFX model version in much the same way as for a pure liquid flow without boiling, but multiplied by the fraction of area unaffected by the bubbles, that is,(2)QC=(1-AW)hC(TW-TL).
Here, hC is the heat transfer coefficient which is written using the temperature wall function T+(y+) known from Kader [11] as
(3)hC=ρCPuτT+,
where nondimensional variables (indicated by superscript “+”) and the friction velocity uτ are defined as usual. Note that (2) and (3) may be evaluated at any location y, provided it is used consistently whenever a variable depends on position. QQ is represented in terms of the quenching heat transfer coefficient hQ:
(4)QQ=AWhQ(TW-TL).
A grid-independent solution for QQ is obtained by evaluating the non-dimensional temperature profile at a fixed value of y+. The evaporation heat flux QE is obtained via the evaporation mass flux at the wall:
(5)QE=m˙WHLG,
where the generated vapour mass  m˙W  is expressed in terms of the bubble diameter at detachment dW, bubble generation frequency f, and nucleation site density  N  as
(6)m˙W=ρGπ6dW3fN.
In terms of bubble detachment diameter and nucleation site density, the wall area fraction AW influenced by vapour bubbles is given by
(7)AW=π(adW2)2N.
Here, a is the so-called bubble influence factor, for which a value of 2 is commonly used [4, 5]. Since AW=1 corresponds to the case where the whole surface is under the influence of bubbles, AW as calculated by (7) has to be limited to values smaller than one. Moreover, it should be kept in mind that already as AW approaches 1 the assumptions of the model are not really satisfied anymore.Correlations for the yet undetermined quantities used in the CFX wall boiling model are discussed in the following.


## 3.2. Bubble Size at Detachment

The bubble size at detachment depends on the liquid subcooling. Also the liquid properties which depend on the system pressure, the flow rate, and the heat flux have an influence.An investigation of the bubble size at detachment was performed by Tolubinsky and Kostanchuk [12] for water at different pressures and subcoolings. The observed dependence on the liquid subcooling at atmospheric pressure can be fitted to a correlation:(8)dW=drefe-((Tsat-TL)/(ΔTrefd)).
A least squares fit gives parameter values are dref=0.0013 m and ΔTrefd=53 K under these conditions.To match the tests of Bartolomej and Chanturiya [7] which were conducted at much higher pressures relevant for typical nuclear energy applications the values of dref and ΔTrefd had to be adjusted to dref=0.6 mm and Trefd=45 K (e.g., [6]).


## 3.3. Nucleation Site Density

The situation concerning data on nucleation site density is much less clear. Most of the time, correlations are expressed in the form of power laws depending on the wall superheat as(9)N=Nref(TW-TsatΔTrefN)p.
However, a recent compilation [13] shows that vastly different parameter values are required to match different data sets. A likely reason for this fact is that nucleation site density is highly dependent on the microscale topography of the boiling surface, which in turn depends strongly on the processes that were used to finish the surface. These processes are very diverse and in most boiling experiments not specifically controlled or characterized.For the tests of  Bartolomej and Chanturiya [7] Nref=0.8×106m-2 and ΔTrefN=10 K were found to yield the best results in the model framework.


## 3.4. Further Model Aspects

The bubble detachment frequencyf  is given according to Cole [14] as a function of the detachment size dW. The quenching heat transfer coefficient is calculated as suggested by Mikic and Rohsenow [15].In the bulk, vapour is assumed to be at saturation condition. Where the liquid is subcooled, that is,TL<Tsat, vapour is condensing and calculated with the transfer rate:
(10)m˙=hLG(Tsat-TL)AIHLG.
Here, hLG is the interfacial heat transfer coefficient, calculated according to Ranz and Marshall [16] and AI is the interfacial area. To close the phase transition model in the bulk bubbly flow with a mean bubble diameter dB, Kurul and Podowski [4] and also Anglart et al. [17] proposed to calculate the bubble diameter dB locally as a linear function of liquid subcooling Tsub:
(11)dB=dB,1(Tsub-Tsub,2)+dB,2(Tsub-Tsub,1)Tsub,2-Tsub,1.
Reference subcooling conditions for typical nuclear energy applications have been given as dB,1=0.1 mm at Tsub,1=-13.5 K and dB,2=2 mm at  Tsub,2=5 K in [17].In a two-phase turbulent flow gaseous bubbles have an influence on the liquid turbulence. This effect is described by increasing the turbulent viscosity according to Sato et al. [18]. Conversely the liquid turbulence flow structures influence the gas distribution, which is described by a turbulent dispersion force according to Burns et al. [19].For momentum exchange between the phases, finally, the Ishii and Zuber [20] drag law was used. Furthermore, a lift force according to Tomiyama et al. [21], a turbulent dispersion force according to Burns et al. [19], and a wall force according to Antal et al. [22] were included.


## 3.5. Evaluated Quantities

Quantities requested for the benchmark have been evaluated as follows. The mixed densityρMIX was calculated according to
(12)ρMIX=αLρL+αGρG,
and the mixed enthalpy HMIX according to
(13)HMIX=αLρLHL+αGρGHGρMIX
with the temperature-dependent static enthalpies for liquid HL and gas HG. The quality X was then calculated according to
(14)X=HMIX-HLSatHGSat-HLSat
with the enthalpies at saturation for liquid HLSat and gas HGSat.


## 4. Model Setup

### 4.1. General Setup

For most of the tests calculations were done for only 1/8 of the channel making use of the geometric symmetries. The mesh consisted of 720000 nodes; 1000 axial subdivisions were considered. To check the validity of this simplification for some tests furthermore 1/4 of the channel (850000 nodes) or the whole (1400000 nodes) were simulated with the same results. An inlet condition was set at the bottom. Upstream of the heated length an unheated part of 0.5 m length was simulated to model flow development. For the liquid and gas properties the water steam tables according to IAPWS IF97 were applied. In Section3 the turbulence of the liquid phase was modelled by a Shear Stress Transport (SST) model [23].  In Section 4 two different Reynolds Stress Turbulence models were applied for comparison. The influence of the turbulence model is investigated in Section 5.


### 4.2. Calibration of Boiling Model Parameters

#### 4.2.1. Bubble Size at Detachment

Earlier investigations [24] have shown that the bubble size at detachment has a sensitive influence on the calculated gas fraction since it determines the part of heat transferred by evaporation. To demonstrate the influence of this parameter different values dref (8) were investigated and the resulting cross sectionally averaged gas fraction for the test case 1.4324 is shown in Figure 2. Nevertheless, because of similar conditions as in the Bartolomej and Chanturiya [7] test cases, dref=0.6 mm was assumed for the further simulations.Figure 2
Cross-sectionally averaged gas fraction for different valuesdref (8).


#### 4.2.2. Nucleation Site Density

The nucleation site density has almost no influence on the liquid temperature, a small influence on the gas void fraction but a strong influence on the wall superheatingTW-Tsat. The common lack of information on the nucleation site density can be compensated by adjusting Nref to match the measured wall temperature. Unfortunately, the wall temperatures are not available in the PSBT data. For the test case 1.4324 different values of the reference nucleation site density Nref (9) were assumed to calculate the circumferentially averaged nucleation site density  N  (see Figure 3), bubble influence area AW (see Figure 4), wall temperature TW (see Figure 5), and the cross sectionally averaged gas volume fraction (see Figure 6). Figure 7 shows furthermore for Nref=8×105m-2  the distribution of the area fraction occupied by gas AW on the heated wall. For axial heights > 1 m (Figure 4) and for different areas (Figure 7) AW>1  is found which does not seem reasonable. Therefore, Nref was set to a smaller value of 8×104 m−2 for all simulated cases. For all simulated cases compliance with the condition AW<1 was checked.Figure 3
Nucleation site densityN  (9).Figure 4
Bubble influence areaAW (7).Figure 5
Wall temperature.Figure 6
Cross sectionally averaged gas fraction.Figure 7
Distribution of the area fraction occupied by gasAW (7) on the heated wall for Nref=8×105m-2 (compare Figure 4). Only 1/8 of the channel is simulated considering the geometric symmetries.


## 4.1. General Setup

For most of the tests calculations were done for only 1/8 of the channel making use of the geometric symmetries. The mesh consisted of 720000 nodes; 1000 axial subdivisions were considered. To check the validity of this simplification for some tests furthermore 1/4 of the channel (850000 nodes) or the whole (1400000 nodes) were simulated with the same results. An inlet condition was set at the bottom. Upstream of the heated length an unheated part of 0.5 m length was simulated to model flow development. For the liquid and gas properties the water steam tables according to IAPWS IF97 were applied. In Section3 the turbulence of the liquid phase was modelled by a Shear Stress Transport (SST) model [23].  In Section 4 two different Reynolds Stress Turbulence models were applied for comparison. The influence of the turbulence model is investigated in Section 5.


## 4.2. Calibration of Boiling Model Parameters

### 4.2.1. Bubble Size at Detachment

Earlier investigations [24] have shown that the bubble size at detachment has a sensitive influence on the calculated gas fraction since it determines the part of heat transferred by evaporation. To demonstrate the influence of this parameter different values dref (8) were investigated and the resulting cross sectionally averaged gas fraction for the test case 1.4324 is shown in Figure 2. Nevertheless, because of similar conditions as in the Bartolomej and Chanturiya [7] test cases, dref=0.6 mm was assumed for the further simulations.Figure 2
Cross-sectionally averaged gas fraction for different valuesdref (8).


### 4.2.2. Nucleation Site Density

The nucleation site density has almost no influence on the liquid temperature, a small influence on the gas void fraction but a strong influence on the wall superheatingTW-Tsat. The common lack of information on the nucleation site density can be compensated by adjusting Nref to match the measured wall temperature. Unfortunately, the wall temperatures are not available in the PSBT data. For the test case 1.4324 different values of the reference nucleation site density Nref (9) were assumed to calculate the circumferentially averaged nucleation site density  N  (see Figure 3), bubble influence area AW (see Figure 4), wall temperature TW (see Figure 5), and the cross sectionally averaged gas volume fraction (see Figure 6). Figure 7 shows furthermore for Nref=8×105m-2  the distribution of the area fraction occupied by gas AW on the heated wall. For axial heights > 1 m (Figure 4) and for different areas (Figure 7) AW>1  is found which does not seem reasonable. Therefore, Nref was set to a smaller value of 8×104 m−2 for all simulated cases. For all simulated cases compliance with the condition AW<1 was checked.Figure 3
Nucleation site densityN  (9).Figure 4
Bubble influence areaAW (7).Figure 5
Wall temperature.Figure 6
Cross sectionally averaged gas fraction.Figure 7
Distribution of the area fraction occupied by gasAW (7) on the heated wall for Nref=8×105m-2 (compare Figure 4). Only 1/8 of the channel is simulated considering the geometric symmetries.


## 4.2.1. Bubble Size at Detachment

Earlier investigations [24] have shown that the bubble size at detachment has a sensitive influence on the calculated gas fraction since it determines the part of heat transferred by evaporation. To demonstrate the influence of this parameter different values dref (8) were investigated and the resulting cross sectionally averaged gas fraction for the test case 1.4324 is shown in Figure 2. Nevertheless, because of similar conditions as in the Bartolomej and Chanturiya [7] test cases, dref=0.6 mm was assumed for the further simulations.Figure 2
Cross-sectionally averaged gas fraction for different valuesdref (8).


## 4.2.2. Nucleation Site Density

The nucleation site density has almost no influence on the liquid temperature, a small influence on the gas void fraction but a strong influence on the wall superheatingTW-Tsat. The common lack of information on the nucleation site density can be compensated by adjusting Nref to match the measured wall temperature. Unfortunately, the wall temperatures are not available in the PSBT data. For the test case 1.4324 different values of the reference nucleation site density Nref (9) were assumed to calculate the circumferentially averaged nucleation site density  N  (see Figure 3), bubble influence area AW (see Figure 4), wall temperature TW (see Figure 5), and the cross sectionally averaged gas volume fraction (see Figure 6). Figure 7 shows furthermore for Nref=8×105m-2  the distribution of the area fraction occupied by gas AW on the heated wall. For axial heights > 1 m (Figure 4) and for different areas (Figure 7) AW>1  is found which does not seem reasonable. Therefore, Nref was set to a smaller value of 8×104 m−2 for all simulated cases. For all simulated cases compliance with the condition AW<1 was checked.Figure 3
Nucleation site densityN  (9).Figure 4
Bubble influence areaAW (7).Figure 5
Wall temperature.Figure 6
Cross sectionally averaged gas fraction.Figure 7
Distribution of the area fraction occupied by gasAW (7) on the heated wall for Nref=8×105m-2 (compare Figure 4). Only 1/8 of the channel is simulated considering the geometric symmetries.


## 5. Results

The calculations for the cross-sectional averages listed in Table2 were performed applying a Shear Stress Transport (SST) model for the liquid phase turbulence. This model blends a k-ε model in the bulk and a k-ω model in the near wall region [23]. Furthermore, the corresponding Reynolds Stress model, the Baseline Reynolds Stress model (BSL-RSM) were investigated. This model has the same model structure as the SST model. In some cases the convergence of the BSL-RSM runs was questionable. The k-ω Reynolds Stress model (O-RSM) proved to be numerically more stable. The details of these approaches to turbulence modeling can be found in the ANSYS solver manual.Table 2
Results.
Run
Pressure [MPa]
Turbulence model
Mixed density [kg m−3]
Void fraction [−]
Quality
Δ
P [Pa m−1]
1.2211
15
Exp
610
0.038
−0.040
BSL-RSM
588.6
0.090
−0.043
17481.9
SST
572.9
0.117
−0.037
17805.2
1.2223
15
Exp
456
0.311
0.040
O-RSM
471.2
0.266
0.073
13685.8
SST
494.5
0.212
0.059
15403.9
1.2237
15
Exp
390
0.440
0.080
O-RSM
431.9
0.327
0.110
12175.5
SST
447.2
0.291
0.101
15725.2
1.2422
15
Exp
522
0.182
0.020
BSL-RSM
520.8
0.187
0.015
8865.2
O-RSM
513.2
0.202
0.017
9363.2
SST
507.3
0.215
0.017
9468.4
1.2423
15
Exp
357
0.508
0.100
O-RSM
424.3
0.346
0.095
8927.7
SST
438
0.325
0.082
9401.9
1.4324
10
Exp
600
0.157
0.000
BSL-RSM
593.8
0.166
−0.002
5751.6
SST
578.9
0.190
0.000
6604.5
1.4325
10
Exp
478
0.335
0.050
BSL-RSM
473
0.336
0.048
9796.8
SST
470
0.343
0.046
10279.4
1.4326
10
Exp
353
0.531
0.110
BSL-RSM
363.3
0.503
0.099
9700
SST
361.9
0.506
0.098
10503.2
1.4411
10
Exp
599
0.152
0.000
BSL-RSM
614.6
0.124
0.000
7506.5
O-RSM
616.1
0.121
0.000
7541.9
SST
617.1
0.122
−0.004
7691
1.6222
5
Exp
549
0.306
0.020
O-RSM
578.9
0.263
0.016
7224.8
Exp: experiment.
BSL-RSM: Baseline Reynolds Stress.
O-RSM: Omega Reynolds Stress.
SST: Shear Stress Transport.Calculating the cross sectionally averaged values, no significant differences between the different turbulence models could be discerned (examples for selected tests see Figure8). In Section 5, however, a strong influence of turbulence modelling on the gas fraction distribution is shown.Figure 8
Cross sectionally averaged values for runs 1.4324, 1.4325, and 1.4326; solid lines: Shear Stress Turbulence model (SST), dotted lines: Reynolds Stress Turbulence model (BSL).The experiments do not provide values on bubble sizes or on temperatures of the heated wall. Exploiting this information gap for calibration of the correlations for the bubble size at detachment and the nucleation site density, at least cross-sectionally averaged values for the mixed density and gas fractions can be calculated with good agreement to experimental data.


## 6. Influence of the Turbulence Model on the Gas Fraction Distribution

Considering the cross-sectional distribution of the gas fraction, remarkable differences between the calculations using different turbulence models were found. Applying a turbulence model with a scalar turbulent viscosity, like the SST-model, the maximum of the gas was found near the heated wall in all cases. Figure9(a) shows the calculated distribution for run 1.2422. In contrast, the measured gas distributions published in the benchmark specifications show significant amounts of gas also in the core of the channel. This is more similar to the results obtained by a Reynolds Stress model (see Figure 9(b)).  The reason for these different patterns is probably that the SST model is not able to capture the secondary flows (see Figure 10).Cross sectionally averaged gas fraction at an axial distance of 1.4 m for run 1.2422.
(a)
Shear Stress Turbulence model (SST)
(b)
Reynolds Stress Turbulence model (O-RST)Velocity stream lines projected in a plane at an axial height of 1.4 m for run 1.2422.
(a)
Shear Stress Turbulence model (SST)
(b)
Reynolds Stress Turbulence model (O-RST)Figure11 shows the comparison of different runs with increasing inlet temperature (compare Table 1). With increasing inlet temperature the maximum gas fraction is shifted from the corner zones towards the core regions of the channel. The figure clearly shows the influence of the liquid flow field on the gas distribution similar to the measurements.Figure 11
Tendencies of gas fraction distribution for the runs 1.4324, 1.4325, and 1.4326.


## 7. Conclusion

Current CFD implementations of a wall boiling model are able to calculate the gas fractions in good agreement with measurements. The model review in Section3 has shown which information is necessary to perform meaningful calculations. With suitably calibrated correlations for the bubble size at detachment and the nucleation site density, at least cross sectionally averaged values for the mixed density and gas fractions can be calculated with good agreement to experimental data. In the corresponding experiments besides the measurement of gas fraction the determination of the temperature at the heated wall is urgent necessarily. Further research to reduce the number of model parameters that have to be determined by model calibration would be highly desirable. The distribution of gas within a cross section has been shown to depend strongly on the modelling of turbulence. Indication has been given that common two-equation models are too simplified for subchannel geometries while anisotropic Reynolds Stress models show a better potential. In order to make a stringent qualification of models, however, much more detailed CFD-grade data are urgently needed. While all of our calculations were performed for steady-state conditions, numerical convergence was a subject requiring special attention which could indicate that the real solution might contain transient elements.



---

*Source: 939561-2012-09-20.xml*