# CFD Modeling of Boiling Flow in PSBT 5×5 Bundle

Authors: Simon Lo; Joseph Osman
Journal: Science and Technology of Nuclear Installations
 (2012)

---

## Abstract

Three-dimensional computational fluid dynamics (CFD) method was used to model the boiling two-phase flow in one of the PSBT 5-by-5 rod bundle tests. The rod bundle with all the spacers was modeled explicitly using unstructured computational grids. The six-equation, two-fluid model with the wall boiling model was used to model the boiling two-phase flows in the bundle. The computed void fractions compare well with the measured data at the measuring plane. In addition to the averaged void data, the CFD results give a very detailed picture of the flow and void distributions in the bundle and how they are affected by solid structures in the flow paths such as the spacer grids and mixing vanes.


---

## Body

## 1. Introduction

In the NUPEC PWR Subchannel and Bundle Test (PSBT) International Benchmark exercise [1], valuable measured data were made available to test and check the accuracy of numerical simulations of boiling two-phase flows in PWR subchannels and rod bundles. The measured data released by Japan Nuclear Energy Safety (JNES) organization [2] were obtained by the Nuclear Power Engineering Corporation (NUPEC) in Japan who in the period between 1987 and 1995 performed a series of void measurement tests using full-size mock-up tests for both BWRs and PWRs. The PWR tests were considered in PSBT. This paper describes the use of three-dimensional computational fluid dynamics (CFD) to model the boiling two-phase flows in one of the 5-by-5 rod bundle tests.The commercial CFD software STAR-CCM+ v6.06 [3] was used in this study. The rod bundle with all the spacers was modeled explicitly using unstructured computational grids. A brief description of the computational grid is given in Section 2 of the paper. The six-equation, two-fluid model with the wall boiling model was used to model the boiling two-phase flows in the bundle. Full details of the mathematical model are provided in Section 3. The steady-state bundle test B5 Run 5.1121 was analysed, and the results are presented in Section 4. The computed averaged void fraction compares well with the measured data at the upper measuring plane. In addition to the averaged void data, the CFD results give a very detailed picture of the flow and void distributions in the bundle and how they are affected by solid structures in the flow paths such as the spacer grids and mixing vanes.


## 2. Computational Grid

The computational grid was generated by recreating the rod bundle and the spacers using the 3D CAD package in STAR-CCM+. The geometries of the rod bundle and the 3 different spacers were taken from the problem specification report by Rubin et al. [1]. CAD models of the 3 spacers were created separately using the CAD package Autodesk Inventor and imported into STAR-CCM+ via parasolid files. A short section of rods going through each spacer was added. An unstructured polyhedral computational grid was then created for each of the combined rod-spacer sections see Figures 1, 2, and 3. The thickness of the spacer grids is represented by 2 layers of computational grids. The springs and dimples in the spacers were included in the model. The contacts between the springs and dimples with the rods were modeled and shown in Figure 4.Figure 1
CFD grid of spacer with mixing vanes.Figure 2
CFD grid of spacer without mixing vanes.Figure 3
CFD grid of simple spacer.Figure 4
Contacts between rods and spacers.The completed rod bundle assembly was created by connecting together the rod-spacer sections according to the specification given in [1]. The connections between the sections were made by extruding the computational cells at the ends of the sections to create horizontal layers of polyhedral cells. These horizontal cell layers start with a smaller height and expand to larger height to economize on the number of computational cells used. The model includes 7 mixing vane spacers (MV), 2 nonmixing vane spacers (NMV), and 8 simple spacers (SS); see Figure 5. In total the CFD model contains 17,950,126 computational cells.Figure 5
Rods and spacers (red = MV, yellow = NMV, blue = SS).


## 3. Mathematical Model

### 3.1. Two-Fluid Model

The standard six-equation, two-fluid model was used in modeling the boiling two-phase flows considered in this paper. In this model the conservation equations for mass, momentum, and energy are solved for both phases.The conservation of mass for phasek is
(1)∂∂t(αkρk)+∇⋅(αkρkuk)=∑i=1N(m˙ki-m˙ik),
where αk is volume fraction of phase k, ρk is phase density, uk is phase velocity, m˙ki and m˙ik are mass transfer rates to and from the phase, and N is the total number of phases. The sum of the volume fractions is clearly equal to unity:
(2)∑kαk=1.
The conservation of momentum for phase k is
(3)∂∂t(αkρkuk)+∇⋅(αkρkukuk)-∇⋅(αk(τk+τkt))=-αk∇p+αkρkg+M,
where τk and τkt are laminar and turbulence shear stresses, p is pressure, and M is the sum of the interfacial forces that include drag and turbulent dispersion forces in this analysis.The conservation of energy for phasek is:
(4)∂∂t(αkρkhk)+∇⋅(αkρkukhk)-∇⋅(αk(λk∇Tk+μtσh∇hk))=Q,
where hk is phase enthalpy, λk is thermal conductivity, T is temperature, μt is turbulent viscosity, σh is turbulent Prandtl number, and Q is the interfacial heat transfer and other heat sources.The standardk-ɛ turbulence model is solved for both phases to represent the flow turbulence and to provide the turbulent viscosities required in the model. Turbulence modeling in multiphase flows is clearly a highly complex area and not well developed or fully understood yet. Interactions between turbulent eddies and steam bubbles need to be taken into consideration. Turbulent eddies will disperse the bubbles. This effect is modeled by the turbulent dispersion force described below. Bubble motions could induce turbulence as well as dissipating turbulence. Bubble-induced turbulence models are available in the literature but not yet considered in this study. It is believed that the turbulence in the flow considered is generated mainly by the geometry and in particular around the spacer regions. The geometry-generated turbulence effects should be captured reasonably well by the grid and the k-ɛ model. The additional effects of bubble-induced turbulence will be investigated in the planned follow-on study in which the more advanced multiphase turbulence models will be considered.


### 3.2. Interfacial Forces

The drag force between the two phases includes a mean and a fluctuation component. The mean drag force is given by(5)FD=34αdρcCDdαcn|ur|ur=ADur,
where CD is the drag coefficient, ur=(uc-ud) is the relative velocity between the two phases and d is the bubble diameter. Subscript c stands for continuous phase and d for dispersed phase. The exponent n=4 is used to model the effects of high bubble concentration, sometimes called the bubble swarm effects.The fluctuating component of the drag force accounts for the additional drag due to interaction between the dispersed phase and the surrounding turbulent eddies. This force is the turbulent dispersion force or the turbulent drag force:(6)FT=-ADνctαdαcσα∇αd,
where νct is the continuous phase turbulent kinematic viscosity and σα is the turbulent Prandtl number; value of 1 is used.The drag coefficient in (5) is computed according to Tomiyama [4] for a contaminated fluid system:
(7)CD=max⁡[24Re(1+0.15Re0.687),8Eo3(Eo+4)].
The Reynolds and Eotvos numbers in (7) are defined as:
(8)Re=ρcurdμc,Eo=|ρc-ρd|gd2σ,
where g is gravitational acceleration and σ is the surface tension coefficient.The interfacial forces would generally include the lift and wall lubrication forces also. The effects of lift and wall lubrication forces are to move the steam bubbles radially away or towards the rod surfaces. Since in this exercise the computed void distributions are to be averaged across the channel, any information on radial void distribution will be lost in the comparison exercise; hence the inclusion of lift and wall lubrication forces is not important and was left out for simplicity.


### 3.3. Wall Boiling Model

At the heated wall, boiling occurs when the wall temperature exceeds the saturation temperature. The steam generation rate is determined by the wall heat partitioning model as follows,(9)qw=ql+qQ+qe,
where, qw is the total heat flux from the wall, ql is the single phase convection heat flux that takes place outside the influence area of nucleation bubbles, qQ is the quenching heat flux within the bubble influence area, and qe  is the evaporation heat flux.The bubble influence areaAe is defined by
(10)Ae=FAπdd24N′′,
where, FA,dd,N′′ are model constant, bubble departure size, and active nucleation site density, respectively. FA=2 is used in this study.The evaporation heat flux can be expressed as(11)qe=πdd36ρghfgfN′′,
where, ρg is the steam density, hfg is the latent heat and f is the bubble departure frequency.The nucleation site density is obtained from Lemmert and Chawla [5]:
(12)N′′=(mΔTsup⁡)p,
where ΔTsup=Twall-Tsat is the wall superheat and m=185 and p=1.805. Twall is the wall temperature and Tsat is the saturation temperature.The bubble departure diameter is obtained from Tolubinsky and Kostanchuk [6]:
(13)dd=d0exp⁡[-ΔTsubΔT0],
where d0=0.6mm and  ΔT0=45 are model constants. ΔTsub=Tsat-Tl  is the liquid subcooling.The bubble departure frequency is obtained from Cole [7]:
(14)f=43g(ρl-ρg)ddρl.
As the bubble detaches from the wall, the space it occupied is filled by cooler water. Part of the wall heat flux is used in heating the replacement water. This heat flux is the quenching heat flux. Del Valle and Kenning [8] modeled this heat transfer as transient heat conduction in a semi-infinite slab:
(15)qQ=hqAe(Twall-Tl),hq=2ftwρlCplλlπ,
where Cpl is the liquid heat capacity and tw is the waiting time between the bubble departure and the activation of the next bubble:
(16)tw=0.8f.
Convective heating of the liquid occurs over the area not covered by nucleation sites. The convective heat flux can be modeled as
(17)ql=hcAc(Twall-Tl),
where the wall heat transfer coefficient hc is obtained from the CFD wall function model and the area factor Ac is obtained from
(18)Ac=1-Ae.
It is important to recognize that the nucleation site density, (12), and the bubble departure diameter, (13), are used together to form the effective boiling area. Correlations for these the two parameters were tuned as one model against available measured data by many investigators for a wide range of conditions. The authors have recently investigated more advanced models for these two parameters that include measured data from low pressure (1 bar) to high pressure (150 bar); see Lo et al [9]. In the Bartolomei test cases studied by Lo et al. [9] the pressure was 147 bar, close to the present study. It was found that calculating the effective boiling area using the more advanced models and the correlations given in (12) and (13) gave practically the same results. Hence the authors believe the models described are suitable for this study.


### 3.4. Bubble Size Distribution

A large range of bubble diameters can be expected in the flow. Since bubbles are generated by boiling and removed by condensation, the bubble diameter is expected to be function of the liquid temperature. Kurul and Podowski [10] defined the local bubble diameter using linear function between measured bubble diameters at two specified values of liquid subcooling:
(19)d=d1(ΔTsub-ΔT0)+d0(ΔT1-ΔTsub)ΔT1-ΔT0,ΔT0≥ΔTsub≥ΔT1,d=d0,ΔTsub≥ΔT0,d=d1,ΔTsub≥ΔT1.
The following values were used in the model:
(20)d0=1.5×10-4atΔT0=13.5,d1=2×10-3atΔT1=-5.


## 3.1. Two-Fluid Model

The standard six-equation, two-fluid model was used in modeling the boiling two-phase flows considered in this paper. In this model the conservation equations for mass, momentum, and energy are solved for both phases.The conservation of mass for phasek is
(1)∂∂t(αkρk)+∇⋅(αkρkuk)=∑i=1N(m˙ki-m˙ik),
where αk is volume fraction of phase k, ρk is phase density, uk is phase velocity, m˙ki and m˙ik are mass transfer rates to and from the phase, and N is the total number of phases. The sum of the volume fractions is clearly equal to unity:
(2)∑kαk=1.
The conservation of momentum for phase k is
(3)∂∂t(αkρkuk)+∇⋅(αkρkukuk)-∇⋅(αk(τk+τkt))=-αk∇p+αkρkg+M,
where τk and τkt are laminar and turbulence shear stresses, p is pressure, and M is the sum of the interfacial forces that include drag and turbulent dispersion forces in this analysis.The conservation of energy for phasek is:
(4)∂∂t(αkρkhk)+∇⋅(αkρkukhk)-∇⋅(αk(λk∇Tk+μtσh∇hk))=Q,
where hk is phase enthalpy, λk is thermal conductivity, T is temperature, μt is turbulent viscosity, σh is turbulent Prandtl number, and Q is the interfacial heat transfer and other heat sources.The standardk-ɛ turbulence model is solved for both phases to represent the flow turbulence and to provide the turbulent viscosities required in the model. Turbulence modeling in multiphase flows is clearly a highly complex area and not well developed or fully understood yet. Interactions between turbulent eddies and steam bubbles need to be taken into consideration. Turbulent eddies will disperse the bubbles. This effect is modeled by the turbulent dispersion force described below. Bubble motions could induce turbulence as well as dissipating turbulence. Bubble-induced turbulence models are available in the literature but not yet considered in this study. It is believed that the turbulence in the flow considered is generated mainly by the geometry and in particular around the spacer regions. The geometry-generated turbulence effects should be captured reasonably well by the grid and the k-ɛ model. The additional effects of bubble-induced turbulence will be investigated in the planned follow-on study in which the more advanced multiphase turbulence models will be considered.


## 3.2. Interfacial Forces

The drag force between the two phases includes a mean and a fluctuation component. The mean drag force is given by(5)FD=34αdρcCDdαcn|ur|ur=ADur,
where CD is the drag coefficient, ur=(uc-ud) is the relative velocity between the two phases and d is the bubble diameter. Subscript c stands for continuous phase and d for dispersed phase. The exponent n=4 is used to model the effects of high bubble concentration, sometimes called the bubble swarm effects.The fluctuating component of the drag force accounts for the additional drag due to interaction between the dispersed phase and the surrounding turbulent eddies. This force is the turbulent dispersion force or the turbulent drag force:(6)FT=-ADνctαdαcσα∇αd,
where νct is the continuous phase turbulent kinematic viscosity and σα is the turbulent Prandtl number; value of 1 is used.The drag coefficient in (5) is computed according to Tomiyama [4] for a contaminated fluid system:
(7)CD=max⁡[24Re(1+0.15Re0.687),8Eo3(Eo+4)].
The Reynolds and Eotvos numbers in (7) are defined as:
(8)Re=ρcurdμc,Eo=|ρc-ρd|gd2σ,
where g is gravitational acceleration and σ is the surface tension coefficient.The interfacial forces would generally include the lift and wall lubrication forces also. The effects of lift and wall lubrication forces are to move the steam bubbles radially away or towards the rod surfaces. Since in this exercise the computed void distributions are to be averaged across the channel, any information on radial void distribution will be lost in the comparison exercise; hence the inclusion of lift and wall lubrication forces is not important and was left out for simplicity.


## 3.3. Wall Boiling Model

At the heated wall, boiling occurs when the wall temperature exceeds the saturation temperature. The steam generation rate is determined by the wall heat partitioning model as follows,(9)qw=ql+qQ+qe,
where, qw is the total heat flux from the wall, ql is the single phase convection heat flux that takes place outside the influence area of nucleation bubbles, qQ is the quenching heat flux within the bubble influence area, and qe  is the evaporation heat flux.The bubble influence areaAe is defined by
(10)Ae=FAπdd24N′′,
where, FA,dd,N′′ are model constant, bubble departure size, and active nucleation site density, respectively. FA=2 is used in this study.The evaporation heat flux can be expressed as(11)qe=πdd36ρghfgfN′′,
where, ρg is the steam density, hfg is the latent heat and f is the bubble departure frequency.The nucleation site density is obtained from Lemmert and Chawla [5]:
(12)N′′=(mΔTsup⁡)p,
where ΔTsup=Twall-Tsat is the wall superheat and m=185 and p=1.805. Twall is the wall temperature and Tsat is the saturation temperature.The bubble departure diameter is obtained from Tolubinsky and Kostanchuk [6]:
(13)dd=d0exp⁡[-ΔTsubΔT0],
where d0=0.6mm and  ΔT0=45 are model constants. ΔTsub=Tsat-Tl  is the liquid subcooling.The bubble departure frequency is obtained from Cole [7]:
(14)f=43g(ρl-ρg)ddρl.
As the bubble detaches from the wall, the space it occupied is filled by cooler water. Part of the wall heat flux is used in heating the replacement water. This heat flux is the quenching heat flux. Del Valle and Kenning [8] modeled this heat transfer as transient heat conduction in a semi-infinite slab:
(15)qQ=hqAe(Twall-Tl),hq=2ftwρlCplλlπ,
where Cpl is the liquid heat capacity and tw is the waiting time between the bubble departure and the activation of the next bubble:
(16)tw=0.8f.
Convective heating of the liquid occurs over the area not covered by nucleation sites. The convective heat flux can be modeled as
(17)ql=hcAc(Twall-Tl),
where the wall heat transfer coefficient hc is obtained from the CFD wall function model and the area factor Ac is obtained from
(18)Ac=1-Ae.
It is important to recognize that the nucleation site density, (12), and the bubble departure diameter, (13), are used together to form the effective boiling area. Correlations for these the two parameters were tuned as one model against available measured data by many investigators for a wide range of conditions. The authors have recently investigated more advanced models for these two parameters that include measured data from low pressure (1 bar) to high pressure (150 bar); see Lo et al [9]. In the Bartolomei test cases studied by Lo et al. [9] the pressure was 147 bar, close to the present study. It was found that calculating the effective boiling area using the more advanced models and the correlations given in (12) and (13) gave practically the same results. Hence the authors believe the models described are suitable for this study.


## 3.4. Bubble Size Distribution

A large range of bubble diameters can be expected in the flow. Since bubbles are generated by boiling and removed by condensation, the bubble diameter is expected to be function of the liquid temperature. Kurul and Podowski [10] defined the local bubble diameter using linear function between measured bubble diameters at two specified values of liquid subcooling:
(19)d=d1(ΔTsub-ΔT0)+d0(ΔT1-ΔTsub)ΔT1-ΔT0,ΔT0≥ΔTsub≥ΔT1,d=d0,ΔTsub≥ΔT0,d=d1,ΔTsub≥ΔT1.
The following values were used in the model:
(20)d0=1.5×10-4atΔT0=13.5,d1=2×10-3atΔT1=-5.


## 4. Results

The steady-state bundle test B5 Run 5.1121 was analysed. The flow conditions were given as pressureP=167.39 (kg/cm2a), mass flux G=14.96 (106 kg/m2hr), power Q=2990 (kW), and inlet temperature Tin=316.9 (°C). The axial power distribution was uniform, and the radial distribution was Type A in which the central pins have the full power and the pins adjacent to the can walls were at 85% power. The measured void fraction averaged over the central 4 subchannels in the upper measuring plane at z=3.177m was 0.1791.The computed results from the CFD model are shown in Figures6 to 9. Figures 6, 7, and 8 show the nonuniformity of void distribution and flow vectors created by the mixing vanes. Figure 9 shows the computed void distributions at the 3 measuring planes. The computed averaged void over the central 4 subchannels in the upper measuring plane is 0.1576 as compared with the measured value of 0.1791. The computed void is therefore 12% lower than the measured value.Figure 6
Void fraction.Figure 7
Void fraction downstream of spacer with mixing vanes.Figure 8
Effects of vanes on flow.Figure 9
Void distributions in the 3 measuring planes.A major advantage of the 3-dimensional CFD method is the level of details it can provide about the flow pattern, void and temperature distributions, and so forth across the whole bundle as shown in Figures6 to 9 or around individual rod as in Figures 10 and 11. Figure 10 shows the effects of the mixing vanes on void distribution. A higher void region is found in the convex side of the vane and a lower void region in the concave side. Detailed flow distribution around the mixing vane can be studied from vector plots like Figure 11.Figure 10
Void distribution above mixing vane spacer around one rod.Figure 11
Velocity distribution above mixing vane spacer around one rod.


## 5. Conclusions

A 3-dimensional CFD model of the PSBT 5-by-5 rod bundle was constructed using the STAR-CCM+ software. Unstructured polyhedral computational cells were used to model the rod bundle and all the spacer grids explicitly. The six-equation two-fluid model together with the wall boiling model was used to model the boiling two-phase flows in the bundle. The steady-state test B5 Run 5.1121 was studied. The computed void fraction averaged over the 4 central subchannels at the upper measuring plane was 0.1576 which is 12% lower than the measured value of 0.1791. This level of agreement between the results is encouraging given the complexity of the geometry and the boiling two-phase flow physics.A major advantage of 3-dimensional CFD is the level of details it can provide about the flow making it possible to perform detailed design analyses for spacer grid and investigating the effect of mixing vanes. However, before this modelling technology will be accepted for design analysis much more rigorous verification and validation of the models are required. In addition to the comparison of channel-averaged results, comparison against detailed spatial distributions of void, velocity, temperature, and so forth, across the whole bundle is required. Hopefully such detailed measured data will become available soon.



---

*Source: 795935-2012-07-03.xml*