# The Feasibility of Multidimensional CFD Applied to Calandria System in the Moderator of CANDU-6 PHWR Using Commercial and Open-Source Codes

Authors: Hyoung Tae Kim; Se-Myong Chang; Jong-Hyeon Shin; Yong Gwon Kim
Journal: Science and Technology of Nuclear Installations
 (2016)

---

## Abstract

The moderator system of CANDU, a prototype of PHWR (pressurized heavy-water reactor), has been modeled in multidimension for the computation based on CFD (computational fluid dynamics) technique. Three CFD codes are tested in modeled hydrothermal systems of heavy-water reactors. Commercial codes, COMSOL Multiphysics and ANSYS-CFX with OpenFOAM, an open-source code, are introduced for the various simplified and practical problems. All the implemented computational codes are tested for a benchmark problem of STERN laboratory experiment with a precise modeling of tubes, compared with each other as well as the measured data and a porous model based on the experimental correlation of pressure drop. Also the effect of turbulence model is discussed for these low Reynolds number flows. As a result, they are shown to be successful for the analysis of three-dimensional numerical models related to the calandria system of CANDU reactors.


---

## Body

## 1. Introduction

The CANDU reactors, developed from 1950s to mid-1960s in Canada, have been adopted in Korea since late 1980s, and now we are operating four pressurized heavy-water reactors constructed in Wolsong nuclear power plant [1]. One of the most important safety issues in the CANDU reactors is the moderator subcooling margin [2] during a large loss of coolant accident to ensure the integrity of the fuel channel submerged in the heavy-water moderator. The moderator subcooling is estimated by the three-dimensional (3D) CFD analysis code which predicts the local distributions of moderator flow as well as temperature. Therefore, the validation of the CFD code is required for accurate prediction of the 3D flow behavior in the moderator system of a CANDU reactor.However, the flow in the system of geometrical complexity in a calandria tank of CANDU reactor is transported through the distributed heat tubes, exerting the pressure drop to the internal coolant flow, so the phenomena should be considered as complicated multiphysics problems. Thanks to the rapid development of CFD (Computational Fluid Dynamics) technology, the 1D model codes have been attempted to be substituted to the 3D simulation, which nowadays becomes not so expensive that we can enjoy the benefit of the innovation [3].In the CANDU-6 reactor equipped with 380 calandria tubes inside the tank of moderator system using pressurized heavy water, the hydrothermal problems have been studied in various research literatures [4, 5]. Modern computational methods are still being developed, and various commercial or open-source codes should be used to simulate physics in the field of nuclear engineering.The moderator system consists of a closed circuit (Figure1), which circulates the heavy water through the calandria and heat exchangers to remove the heat generated in the moderator during reactor operation. The hot heavy water is drawn from the bottom of the calandria through two outlet pipes, pumped by either of the two pumps, and discharged into the two heat exchangers. The cooled heavy water from the two heat exchangers is returned to the calandria through eight inlet nozzles directed upwards and located on diametrically opposite sides of the calandria vessel at the horizontal mid-plane. In this study the modeling of the moderator flow is only considered within the calandria. Since a square array of horizontal calandria tubes (outside boundary of the fuel channels) is located inside the calandria, it is necessary to investigate the moderator flow around the rod bundle.Figure 1
Configuration of the moderator circulation test.In this paper, we set up the 3D numerical method of analysis using three kinds of codes: the commercial codes, COMSOL Multiphysics [6], and ANSYS-CFX [7] are used for validation as well as the solution of specific problems, and OpenFOAM [8], an open-source code, is remodeled for the compatibility of the commercial ones. The calandria tank system is first modeled as a simplified multidimensional geometry to see the essential physics and to test the feasibility of using the present numerical models.


## 2. Numerical Method

### 2.1. Governing Equations

The hydraulic governing equations based on the single-phase incompressible flow are written in the vector form:(1)
∇
·
V
=
0
,
(2)
ρ
∂
V
∂
t
+
V
·
∇
V
=
-
∇
p
+
ρ
g
+
μ
∇
2
V
+
f
V
,where V, ρ, and p are velocity vector, density, and pressure while the constants μ and f
V are dynamic viscosity and body force per unit volume. Equation (1) is the continuity equation for incompressible flow, and the Navier-Stokes momentum equation (2) is decoupled from energy equation in the source term, f
V, or buoyancy force from Boussinesq approximation, which is neglected in the present implementation.


### 2.2. Turbulence Modeling

As the flow is so fast where the Reynolds number is more than critical one,1
0
4 to 1
0
5, which is low, the effect of turbulence cannot be neglected. Jet flows through tube bundles are generally unsteady with boundary layer separation and vortex shedding around the tubes.Standardk
-
ε model performs poorly for separated flows due to its overly diffusive nature, but it is robust for mesh types, so preferred in general [9]. Also we used other turbulence models such as k
-
ω, SST (Shear Stress Transport), and Spalart-Allmaras model for the comparison of result [10].


### 2.3. Boundary Conditions

The essential boundary conditions in this problem are listed as follows:(1)
Velocities: no-slip conditions at walls, and the inlet velocity is specified from the volume flow rate of the moderator system.
(2)
Pressure: zero pressure gradient conditions at walls and inlet, which should be valid under the assumption that the thickness of boundary layer is very thin. The outlet pressure is fixed by the moderator system.


### 2.4. Numerical Methods

For the commercial code COMSOL, FEM (finite element method) based numerical schemes are used: a segregated algorithm for the pressure term with FGMRES linearization and Petrov-Galerkin least-square artificial dissipation techniques [6]. For the commercial code ANSYS-CFX, element-based FVM (finite volume method) schemes are applied for the vertex-centered complicated geometry using shape function concept [11]. In the computation using OpenFOAM, SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm, a kind of FVM is applied for the iteration until the steady state for (1)~(2). In this method, the pressure gradient term in (2) is isolated, and subiterations should be performed between predictor and corrector [12].


## 2.1. Governing Equations

The hydraulic governing equations based on the single-phase incompressible flow are written in the vector form:(1)
∇
·
V
=
0
,
(2)
ρ
∂
V
∂
t
+
V
·
∇
V
=
-
∇
p
+
ρ
g
+
μ
∇
2
V
+
f
V
,where V, ρ, and p are velocity vector, density, and pressure while the constants μ and f
V are dynamic viscosity and body force per unit volume. Equation (1) is the continuity equation for incompressible flow, and the Navier-Stokes momentum equation (2) is decoupled from energy equation in the source term, f
V, or buoyancy force from Boussinesq approximation, which is neglected in the present implementation.


## 2.2. Turbulence Modeling

As the flow is so fast where the Reynolds number is more than critical one,1
0
4 to 1
0
5, which is low, the effect of turbulence cannot be neglected. Jet flows through tube bundles are generally unsteady with boundary layer separation and vortex shedding around the tubes.Standardk
-
ε model performs poorly for separated flows due to its overly diffusive nature, but it is robust for mesh types, so preferred in general [9]. Also we used other turbulence models such as k
-
ω, SST (Shear Stress Transport), and Spalart-Allmaras model for the comparison of result [10].


## 2.3. Boundary Conditions

The essential boundary conditions in this problem are listed as follows:(1)
Velocities: no-slip conditions at walls, and the inlet velocity is specified from the volume flow rate of the moderator system.
(2)
Pressure: zero pressure gradient conditions at walls and inlet, which should be valid under the assumption that the thickness of boundary layer is very thin. The outlet pressure is fixed by the moderator system.


## 2.4. Numerical Methods

For the commercial code COMSOL, FEM (finite element method) based numerical schemes are used: a segregated algorithm for the pressure term with FGMRES linearization and Petrov-Galerkin least-square artificial dissipation techniques [6]. For the commercial code ANSYS-CFX, element-based FVM (finite volume method) schemes are applied for the vertex-centered complicated geometry using shape function concept [11]. In the computation using OpenFOAM, SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm, a kind of FVM is applied for the iteration until the steady state for (1)~(2). In this method, the pressure gradient term in (2) is isolated, and subiterations should be performed between predictor and corrector [12].


## 3. Multidimensional Implementation

In this study, three codes are used for the implementation of models for the calandria moderator system: COMSOL Multiphysics (version 4.4), ANSYS-CFX (version 14.5), and OpenFOAM (version 2.3.1). The first two are commercial codes while the last one is an open-source code.


### 3.1. COMSOL Multiphysics

Despite the deficiency for memory control, COMSOL has been widely used for multiphysical analyses thanks to its accuracy and compatibility for variable kinds of physics. From version 4, the CFD module becomes independent to be powerful for convergence and stability. Two-dimensional computation is performed for the test problem and extendable to 3D for simple problems if the computer memory is sufficient.Figure1 is the benchmark result for a cylinder-bank problem with in-line 24 rows where the diameter of each tube is d
=
4.8 mm, and the pitch is l
=
10.37 mm where the driver medium is air. Tubes are modeled precisely to get a high resolution of computational result with total 640,740 prism-type meshes for the quarter symmetric domain of computation: see Figure 2(a). The velocity field is visualized in Figure 2(b), and the pressure normalized as a dimensionless parameter, C
p
≡
p
-
p
r
e
f
/
ρ
V
m
2
/
2, where p
r
e
f and V
m are the reference pressure at the first-row station and the incident speed normal to the cylinders, respectively. The normalized pressure drop is plotted with the counterpart experimental data in Figure 3 where the horizontal distance is normalized with the distance between centers of two adjacent cylinders, x
/
S
x.Figure 2
COMSOL computation: (a) grids system and (b) velocity field.
(a)
(b)Figure 3
The pressure drop along with the tapping line. Computation of normalized pressure drop along with the lower edge line of Figure2(b) with COMSOL, which is compared with the experimental result, courtesy to Professor Tongbeum Kim, University of the Witwatersrand, South Africa.


### 3.2. ANSYS-CFX

The strong point of ANSYS-CFX is its high adaption to the complex geometry and quick convergence using multigrids technique and parallel processing. The use of shape function used in FEM for each cell or the element-based FVM makes it possible to construct a vertex-centered scheme, which is contrast to ANSYS-FLUENT using a cell-centered FVM [7]. Figures 4(a) and 4(b) show an example of full-scale computation for the 2D geometry model with a buoyancy term in momentum equation, (2) where the total number of meshes is 672,912 [13].Figure 4
An example of numerical result for the computation of a 2D moderator model of CANDU-6 reactor with ANSYS-CFX: a cross-sectional view, (a) velocity vectors (partial), and (b) process for grid construction [13].
(a)
(b)y
+ value is defined for Δ
y, and the normal length of grid at the first one neighboring the wall is (3)
y
+
=
Δ
y
ρ
μ
Δ
u
t
Δ
y
y
=
y
w
,where u
t is the tangentially transformed velocity component along the wall [7]. To get the proper result of computation without wall functions, the dimensionless wall distance should be guaranteed y
+
<
1 as in whole computational domain. However, the required number of grids can be much reduced with use of wall functions [7]:(4)
y
+
=
m
a
x
ρ
μ
Δ
y
C
μ
1
/
4
k
,
11.06
,where k is the turbulence kinetic energy and C
μ
=
0.09 is a coefficient. Therefore, the grids in the range of y
+
<
11.06 are neglected to be unnecessary, and y
+ less than 30 is sufficient for the capture of outer layer in a low Reynolds number flow for engineering computation. Therefore, y
+
<
30 is guaranteed for all the computational domains in Figure 4(b) with safety margin.


### 3.3. OpenFOAM

OpenFOAM (Open Field Operation and Manipulation) has been developed by Henry Weller and Hrvoje Jasak in Imperial College. The source code has been opened to the public since 2004. This code is operated on the Linux based O/S such as Ubuntu, so the copyright is absolutely free for every CFD program developer. This code is originated from the OOP (object oriented programming) concept based on C++ program language. Solvers and libraries are defined as C++ classes. With the postprocessor ParaView, the graphical visualization becomes possible with a commandparaFoam [8].


#### 3.3.1. 3D Grid Generation

There are two kinds of methods for the grid generation in OpenFOAM. WithblockMesh command, the OpenFOAM can produce a grid system as an ASCII code file under the specific syntax. However, for the complicated grids, the transformation from the specific format supported by a commercial code to the OpenFOAM mesh is much more efficient. To produce the FLUENT mesh, the ANSYS ICEM CFD software compatible with ANSYS-CFX is used from which we produce a ∗.cfx format file, and then this file is converted to ∗.cas (FLUENT mesh format, ASCII). With a commandfluentMeshToFoam supported by OpenFOAM, the grid data also generates boundaries and cell connectivity: see Figure 5(a). Figure 5(b) is an example view of 3D mesh.Figure 5
3D grids system in OpenFOAM: (a) marked with boundaries and (b) example view of full grids system, approximately 200 k hexahedral meshes.
(a)
(b)


#### 3.3.2. Computational System

With the simpleFoam solver, the 3D moderator model in Figure5 is preliminarily computed in Figures 6(a) and 6(b), velocity and pressure per unit density fields for the given grids presenting an asymmetric distribution. Figure 7 shows the folder structure where the grids file for mesh information is located in thepolyMesh directory. Other numerical constraints are restored at text-mode files in the corresponding folders such as initial or boundary conditions, turbulence model, material properties, solver-control numerical coefficients, computational scheme, and tolerance for steady state.Figure 6
3D computation result using OpenFOAM fields of (a) velocity,V
m
a
x
=
0.5073 m/s, and (b) pressure per unit density, p
/
ρ
m
a
x
=
0.2686 (m/s)2 gauge.
(a)
(b)Figure 7
Schematic of the folder structure for OpenFOAM computation.The user can select a proper solver that is supported by OpenFOAM classes. Even he/she can compose a new solver based on the existing solvers as the C++ open-source code is edited in the environment of Linux operating system. Table1 is the list of solvers used in the present research.Table 1
OpenFOAM solvers used in the present computation.
Types
Name of solver
Related physics
Incompressible(isothermal)
(porousSimpleFoam)simpleFoam
(Porous media)Steady/unsteady-stateIncompressibleTurbulent (model)
Heat transfer
buoyantBoussinesqSimpleFoam
Steady/unsteady-stateBuoyant: free convectionIncompressibleTurbulent (model)In all computation with OpenFOAM, the numerical algorithm is based on SIMPLE scheme, and the relaxation factors are controlled from 0.3 (pressure) to 0.7 (others), under-relaxation for the control of computational stability, and tolerance for convergence criteria for SIMPLE algorithm is1
×
1
0
-
3 (pressure) and 1
×
1
0
-
4 (others) in the relative values for the total residual during time marching and also global steady state.


## 3.1. COMSOL Multiphysics

Despite the deficiency for memory control, COMSOL has been widely used for multiphysical analyses thanks to its accuracy and compatibility for variable kinds of physics. From version 4, the CFD module becomes independent to be powerful for convergence and stability. Two-dimensional computation is performed for the test problem and extendable to 3D for simple problems if the computer memory is sufficient.Figure1 is the benchmark result for a cylinder-bank problem with in-line 24 rows where the diameter of each tube is d
=
4.8 mm, and the pitch is l
=
10.37 mm where the driver medium is air. Tubes are modeled precisely to get a high resolution of computational result with total 640,740 prism-type meshes for the quarter symmetric domain of computation: see Figure 2(a). The velocity field is visualized in Figure 2(b), and the pressure normalized as a dimensionless parameter, C
p
≡
p
-
p
r
e
f
/
ρ
V
m
2
/
2, where p
r
e
f and V
m are the reference pressure at the first-row station and the incident speed normal to the cylinders, respectively. The normalized pressure drop is plotted with the counterpart experimental data in Figure 3 where the horizontal distance is normalized with the distance between centers of two adjacent cylinders, x
/
S
x.Figure 2
COMSOL computation: (a) grids system and (b) velocity field.
(a)
(b)Figure 3
The pressure drop along with the tapping line. Computation of normalized pressure drop along with the lower edge line of Figure2(b) with COMSOL, which is compared with the experimental result, courtesy to Professor Tongbeum Kim, University of the Witwatersrand, South Africa.


## 3.2. ANSYS-CFX

The strong point of ANSYS-CFX is its high adaption to the complex geometry and quick convergence using multigrids technique and parallel processing. The use of shape function used in FEM for each cell or the element-based FVM makes it possible to construct a vertex-centered scheme, which is contrast to ANSYS-FLUENT using a cell-centered FVM [7]. Figures 4(a) and 4(b) show an example of full-scale computation for the 2D geometry model with a buoyancy term in momentum equation, (2) where the total number of meshes is 672,912 [13].Figure 4
An example of numerical result for the computation of a 2D moderator model of CANDU-6 reactor with ANSYS-CFX: a cross-sectional view, (a) velocity vectors (partial), and (b) process for grid construction [13].
(a)
(b)y
+ value is defined for Δ
y, and the normal length of grid at the first one neighboring the wall is (3)
y
+
=
Δ
y
ρ
μ
Δ
u
t
Δ
y
y
=
y
w
,where u
t is the tangentially transformed velocity component along the wall [7]. To get the proper result of computation without wall functions, the dimensionless wall distance should be guaranteed y
+
<
1 as in whole computational domain. However, the required number of grids can be much reduced with use of wall functions [7]:(4)
y
+
=
m
a
x
ρ
μ
Δ
y
C
μ
1
/
4
k
,
11.06
,where k is the turbulence kinetic energy and C
μ
=
0.09 is a coefficient. Therefore, the grids in the range of y
+
<
11.06 are neglected to be unnecessary, and y
+ less than 30 is sufficient for the capture of outer layer in a low Reynolds number flow for engineering computation. Therefore, y
+
<
30 is guaranteed for all the computational domains in Figure 4(b) with safety margin.


## 3.3. OpenFOAM

OpenFOAM (Open Field Operation and Manipulation) has been developed by Henry Weller and Hrvoje Jasak in Imperial College. The source code has been opened to the public since 2004. This code is operated on the Linux based O/S such as Ubuntu, so the copyright is absolutely free for every CFD program developer. This code is originated from the OOP (object oriented programming) concept based on C++ program language. Solvers and libraries are defined as C++ classes. With the postprocessor ParaView, the graphical visualization becomes possible with a commandparaFoam [8].


### 3.3.1. 3D Grid Generation

There are two kinds of methods for the grid generation in OpenFOAM. WithblockMesh command, the OpenFOAM can produce a grid system as an ASCII code file under the specific syntax. However, for the complicated grids, the transformation from the specific format supported by a commercial code to the OpenFOAM mesh is much more efficient. To produce the FLUENT mesh, the ANSYS ICEM CFD software compatible with ANSYS-CFX is used from which we produce a ∗.cfx format file, and then this file is converted to ∗.cas (FLUENT mesh format, ASCII). With a commandfluentMeshToFoam supported by OpenFOAM, the grid data also generates boundaries and cell connectivity: see Figure 5(a). Figure 5(b) is an example view of 3D mesh.Figure 5
3D grids system in OpenFOAM: (a) marked with boundaries and (b) example view of full grids system, approximately 200 k hexahedral meshes.
(a)
(b)


### 3.3.2. Computational System

With the simpleFoam solver, the 3D moderator model in Figure5 is preliminarily computed in Figures 6(a) and 6(b), velocity and pressure per unit density fields for the given grids presenting an asymmetric distribution. Figure 7 shows the folder structure where the grids file for mesh information is located in thepolyMesh directory. Other numerical constraints are restored at text-mode files in the corresponding folders such as initial or boundary conditions, turbulence model, material properties, solver-control numerical coefficients, computational scheme, and tolerance for steady state.Figure 6
3D computation result using OpenFOAM fields of (a) velocity,V
m
a
x
=
0.5073 m/s, and (b) pressure per unit density, p
/
ρ
m
a
x
=
0.2686 (m/s)2 gauge.
(a)
(b)Figure 7
Schematic of the folder structure for OpenFOAM computation.The user can select a proper solver that is supported by OpenFOAM classes. Even he/she can compose a new solver based on the existing solvers as the C++ open-source code is edited in the environment of Linux operating system. Table1 is the list of solvers used in the present research.Table 1
OpenFOAM solvers used in the present computation.
Types
Name of solver
Related physics
Incompressible(isothermal)
(porousSimpleFoam)simpleFoam
(Porous media)Steady/unsteady-stateIncompressibleTurbulent (model)
Heat transfer
buoyantBoussinesqSimpleFoam
Steady/unsteady-stateBuoyant: free convectionIncompressibleTurbulent (model)In all computation with OpenFOAM, the numerical algorithm is based on SIMPLE scheme, and the relaxation factors are controlled from 0.3 (pressure) to 0.7 (others), under-relaxation for the control of computational stability, and tolerance for convergence criteria for SIMPLE algorithm is1
×
1
0
-
3 (pressure) and 1
×
1
0
-
4 (others) in the relative values for the total residual during time marching and also global steady state.


## 3.3.1. 3D Grid Generation

There are two kinds of methods for the grid generation in OpenFOAM. WithblockMesh command, the OpenFOAM can produce a grid system as an ASCII code file under the specific syntax. However, for the complicated grids, the transformation from the specific format supported by a commercial code to the OpenFOAM mesh is much more efficient. To produce the FLUENT mesh, the ANSYS ICEM CFD software compatible with ANSYS-CFX is used from which we produce a ∗.cfx format file, and then this file is converted to ∗.cas (FLUENT mesh format, ASCII). With a commandfluentMeshToFoam supported by OpenFOAM, the grid data also generates boundaries and cell connectivity: see Figure 5(a). Figure 5(b) is an example view of 3D mesh.Figure 5
3D grids system in OpenFOAM: (a) marked with boundaries and (b) example view of full grids system, approximately 200 k hexahedral meshes.
(a)
(b)


## 3.3.2. Computational System

With the simpleFoam solver, the 3D moderator model in Figure5 is preliminarily computed in Figures 6(a) and 6(b), velocity and pressure per unit density fields for the given grids presenting an asymmetric distribution. Figure 7 shows the folder structure where the grids file for mesh information is located in thepolyMesh directory. Other numerical constraints are restored at text-mode files in the corresponding folders such as initial or boundary conditions, turbulence model, material properties, solver-control numerical coefficients, computational scheme, and tolerance for steady state.Figure 6
3D computation result using OpenFOAM fields of (a) velocity,V
m
a
x
=
0.5073 m/s, and (b) pressure per unit density, p
/
ρ
m
a
x
=
0.2686 (m/s)2 gauge.
(a)
(b)Figure 7
Schematic of the folder structure for OpenFOAM computation.The user can select a proper solver that is supported by OpenFOAM classes. Even he/she can compose a new solver based on the existing solvers as the C++ open-source code is edited in the environment of Linux operating system. Table1 is the list of solvers used in the present research.Table 1
OpenFOAM solvers used in the present computation.
Types
Name of solver
Related physics
Incompressible(isothermal)
(porousSimpleFoam)simpleFoam
(Porous media)Steady/unsteady-stateIncompressibleTurbulent (model)
Heat transfer
buoyantBoussinesqSimpleFoam
Steady/unsteady-stateBuoyant: free convectionIncompressibleTurbulent (model)In all computation with OpenFOAM, the numerical algorithm is based on SIMPLE scheme, and the relaxation factors are controlled from 0.3 (pressure) to 0.7 (others), under-relaxation for the control of computational stability, and tolerance for convergence criteria for SIMPLE algorithm is1
×
1
0
-
3 (pressure) and 1
×
1
0
-
4 (others) in the relative values for the total residual during time marching and also global steady state.


## 4. Benchmark Validation

### 4.1. Definition and Problem

A benchmark test for the performance of each code is proposed for the well-known STERN laboratory experiment [8]. This problem is often used for the comparison with CFD results [10, 13, 14] and uses a reduced-scale CANDU-6 moderator model: see Figure 8(a). The central part of tube bundles is isolated with flat plates, and the pressure drop is precisely measured from pressure taps located in the three stations, from PT1 to PT3 in Figure 8(b) where V
m is the mean flow velocity incident to the rectangular channel. The diameter of each tube is 33.02 mm; the spacing is 71.4 mm: 4 × 24 arrays of cylinders, and the dimension of channel is 2,000 (length) × 285.6 (width) × 200 (depth) in mm unit. Two columns of square blocks do not contain the cylinder at inlet and outlet in Figure 8(b). The experimental conditions are listed in Table 2.Table 2
Three cases in STERN laboratory experiment,d
=
33.02
m
m.
Case
Physical properties of heavy water
V
m [m/s]
ρ [kg/m3]
μ [Ns/m2]
R
e
d
=
ρ
V
m
d
/
μ
#1 (T.P. 326)
0.054
992.25
0.000653
2,709
#2 (T.P. 299)
0.070
981.00
0.000440
5,153
#3 (T.P. 306)
0.103
971.60
0.000355
9,308Figure 8
STERN laboratory experiment: (a) experimental configuration and (b) the measurement of pressure [13, 15].
(a)
(b)The Reynolds number in Table2 shows that the flow regime lies in the flow instability region, not the fully turbulence. However, the incident flow usually contains some inherent turbulence even if the boundary condition of an experiment is very well controlled. The intake turbulence is defined with the variables of k
-
ε model as(5)
k
=
1
2
V
′
2
,
(6)
ε
=
C
μ
0.75
k
1.5
d
,where V
′ is the mean time-perturbation amount of flow velocity [14]. The inlet turbulence intensity in (5) is assumed to be 5%, which is a general value in many experiments, and the reference distance level is the diameter of a cylinder in (6).


### 4.2. Comparison of Results

Figure9 shows the mesh generation procedure for the pressure drop test. The mesh generation tool used in the present work is the ANSYS ICEM CFD [7]. Only two rows of cylinders are considered as the symmetric boundary condition is imposed on the top surface. To capture the local flow in the boundary layer around the rod bundles, fine mesh layers are generated as a hexagonal surface mesh, and this 2D mesh in the region of the tube pitch is extruded to obtain the 3D mesh. At the wall, y
+ value in (3) is checked to y
+
<
30. Half of the test section is modeled by applying the symmetric boundary condition on the upper plane of computational domain to enhance the efficiency of computation with a reduced number of mesh elements. A series of convergence test for grid has been done to make sure that the number of grids should be sufficient for the precise solution. Finally, the total number of hexagonal elements is estimated as 424,320. The grid model in Figure 9 is used commonly to 3D computations using ANSYS-CFX and OpenFOAM. In the reduced 2D computation, the grid in the depth direction is merged as a cell for OpenFOAM.Figure 9
Grid system for the benchmark validation.As shown in Teyssedou et al. [10], the 2D computation tends to underestimate the pressure drop, and all codes using k
-
ε turbulence model presented the same trend and similar order of qualitative values: see Table 3. The hydraulic diameters in the 2D model are far different from the original 3D experimental setup, so flow resistance in the depth direction cannot be neglected, generating additional pressure drop. The velocity fields visualized in Figure 10 show that the 2D flow in Figure 10(a) is different from the 3D flow in Figures 10(b) and 10(c), characterized as the irregularity of velocity vectors in the depth direction.Table 3
Comparison of pressure drop for the STERN laboratory experiment, 2D.
Case
Δ
p (Pa)
Experiment(STERN)
COMSOL(2D)
ANSYS-FLUENT (2D)∗
OpenFOAM(2D)
#1
28.2
15.0
16.0
14.0
#2
41.3
22.5
20.3
21.4
#3
78.2
45.2
36.4
38.0
∗Quoted for the data with k
-
ε turbulence model in FLUENT version 15 [10].Figure 10
Velocity fields from the computation with various codes, Case  3: (a) COMSOL (2D), (b) ANSYS-CFX (3D), and (c) OpenFOAM (3D).
(a)
(b)
(c)In the 3D computation, Table4, the error between computational and experimental value is far reduced, but there is still the discord among the experimental result and computations. The cause of error might be originated from the limitation of turbulence model because the flow regime lies in the low Reynolds number region from 2,709 to 9,308. The Kármán instability of the stack of cylinders makes the numerical solution difficult to converge for the case of 2D computation where often the pressure drop is underestimated [10] especially in the low Reynolds number region, but this can be far more relaxed with the multidimensional turbulence effect in 3D results because of the increase in degrees of freedom. Case  3 flow field is compared for ANSYS-CFX and OpenFOAM in Figures 10(b) and 10(c), equivalently.Table 4
Comparison of pressure drop for the STERN laboratory experiment, 3D.
Case
Δ
p (Pa)
Experiment(STERN)
Equation (7)∗
COMSOL(3-D)
ANSYS-CFX (3-D)
OpenFOAM(3-D)
#1
28.2
27.0
22.3
23.6
23.6
#2
41.3
40.1
36.3
37.6
38.3
#3
78.2
77.7
76.9
77.5
79.9
∗The porous model is used with the experimental correlation [15].Hadaller et al. [15] developed the experimental pressure drop correlations of this problem for aligned and staggered arrangement cylinders. In the aligned or in-line correlation, the pressure drop per unit channel length is(7)
Δ
p
Δ
L
=
4.54
R
e
d
-
0.172
l
cos
⁡
θ
1
2
ρ
V
m
2
;
R
e
d
=
ρ
V
m
d
μ
,where l
=
71.4 mm is the pitch of cylinders and θ is the incident angle of free stream, zero for this problem. The free-stream flow velocity is defined as V
m
=
ε
p
V, and the velocity magnitude V
=
∑
i
V
i
2 is measured inside the tube bank [2]. The porosity in this problem is defined as [14](8)
ε
p
=
1
-
π
/
4
d
2
l
2
.In the comparison of pressure drops,Δ
p between PT1 and PT3 in Figure 8(b) from all the codes is given in Table 4. ANSYS-CFX and OpenFOAM data are from the direct simulation of (1)~(2) while (7) is used for the corrected body force per unit volume in the momentum equation (2): that is, f
V
,
x
=
-
Δ
p
/
Δ
L in the source term of the COMSOL code [14]. The correlation in (7) shows a good agreement with the experimental data, but the numerical data also result in comparable agreement. Figure 11 is the relative errors from the experimental data and shows the full performance of codes for the present benchmark problem. All the computational results except for the OpenFOAM of Case  3 predict underestimating of the experimental data. Some artificial factors like surface roughness or incident turbulence could have effect on the experimental date for the increase of pressure drop.Figure 11
Comparison of relative errors for various 3D codes from the experimental data.However, in spite of the uncertainty of experiment, Case  1 shows the worst prediction of numerical codes (20.9% for COMOSOL; 16.3% for both ANSYS-CFX and OpenFOAM), which is strongly suspected to the defect of turbulence model at low Reynolds number as 2,709. The laminar computation with no turbulence model is never converged to the steady solution but oscillates within 100 times of the tolerance of error. For Case  1, the turbulent solution is converged att
=
1,451 s with OpenFOAM. The pressure drop between PT1 and PT3 in Figure 8(b) is calculated at t
=
10,000 s for laminar flow where the predicted value is Δ
p
=
24.6 Pa, and the error from experimental data is mitigated to 12.8%. The wake flow fields computed from laminar and turbulent model are compared for Case  1 streamlines in Figures 12(a) and 12(b). In the laminar model, as the shear layer bounding on the separated flow from cylinders blocks the flow along the centerline of symmetry, that is, the upper plane, the flow resistance makes the horizontal pressure gradient slightly increase more than that of the turbulent model.Figure 12
Comparison of wake flow field in the streamlines for Case  1, OpenFOAM with (a) laminar, or no turbulence model, and (b)k
-
ε turbulence model.
(a)
(b)


### 4.3. Sensitivity to Turbulence Models

Various turbulence models are tested for the same benchmark problems in Section4.2. 2D models with COMSOL and 3D models with ANSYS-CFX are presented in Tables 5 and 6, respectively. Four models are tested such as standard k
-
ε, k
-
ω, SST (shear stress transport) [16], and SA (Spalart-Allmaras) models [17]. As shown overall, the k
-
ε model shows the most reliable agreement with the experimental data: see Tables 5 and 6. Reference [10] also shows that the error of standard model is less than those of various k
-
ε models: R, RNG, and so forth, although they are all 2D computations.Table 5
Sensitivity to the turbulence models, 2D (COMSOL).
Case
Δ
p (Pa)
k
-
ε
k
-
ω
SST
Spalart-Allmaras
#1
15.0
12.3
9.5
11.4
#2
22.5
17.3
15.1
—∗
#3
45.2
33.0
31.0
—∗
∗They do not converge under the present grids set.Table 6
Sensitivity to the turbulence models, 3D (ANSYS-CFX).
Case
Δ
p (Pa)
k
-
ε
k
-
ω
SST
Spalart-Allmaras
#1
23.6
21.8
20.6
26.7
#2
37.6
34.7
31.8
42.1
#3
77.5
74.8
73.9
115.2In the agreement of pressure drop in the STERN laboratory test, thek
-
ε model seems better than the k
-
ω model for both 2D and 3D, which is thought to the fact that they are using different wall functions for “the enhanced wall treatment” inside the inner boundary layer (y
+
≈
1) [10]. The SST model seems not a good solution for this problem because it is originally based on the k
-
ω model. In the 3D low Reynolds number case (Case  1), the error from the experimental data is decreased to 5.3% with the one-equation Spalat-Allmaras model in Table 6. However, it is not accurate for the case of the highest Reynolds number (Case  3) as shown in Table 6 and even it does not converge to a steady solution (Case  2 and 3) for 2D in Table 5, where the present grids system is not adequate or too coarse for the one-equation turbulence model.


## 4.1. Definition and Problem

A benchmark test for the performance of each code is proposed for the well-known STERN laboratory experiment [8]. This problem is often used for the comparison with CFD results [10, 13, 14] and uses a reduced-scale CANDU-6 moderator model: see Figure 8(a). The central part of tube bundles is isolated with flat plates, and the pressure drop is precisely measured from pressure taps located in the three stations, from PT1 to PT3 in Figure 8(b) where V
m is the mean flow velocity incident to the rectangular channel. The diameter of each tube is 33.02 mm; the spacing is 71.4 mm: 4 × 24 arrays of cylinders, and the dimension of channel is 2,000 (length) × 285.6 (width) × 200 (depth) in mm unit. Two columns of square blocks do not contain the cylinder at inlet and outlet in Figure 8(b). The experimental conditions are listed in Table 2.Table 2
Three cases in STERN laboratory experiment,d
=
33.02
m
m.
Case
Physical properties of heavy water
V
m [m/s]
ρ [kg/m3]
μ [Ns/m2]
R
e
d
=
ρ
V
m
d
/
μ
#1 (T.P. 326)
0.054
992.25
0.000653
2,709
#2 (T.P. 299)
0.070
981.00
0.000440
5,153
#3 (T.P. 306)
0.103
971.60
0.000355
9,308Figure 8
STERN laboratory experiment: (a) experimental configuration and (b) the measurement of pressure [13, 15].
(a)
(b)The Reynolds number in Table2 shows that the flow regime lies in the flow instability region, not the fully turbulence. However, the incident flow usually contains some inherent turbulence even if the boundary condition of an experiment is very well controlled. The intake turbulence is defined with the variables of k
-
ε model as(5)
k
=
1
2
V
′
2
,
(6)
ε
=
C
μ
0.75
k
1.5
d
,where V
′ is the mean time-perturbation amount of flow velocity [14]. The inlet turbulence intensity in (5) is assumed to be 5%, which is a general value in many experiments, and the reference distance level is the diameter of a cylinder in (6).


## 4.2. Comparison of Results

Figure9 shows the mesh generation procedure for the pressure drop test. The mesh generation tool used in the present work is the ANSYS ICEM CFD [7]. Only two rows of cylinders are considered as the symmetric boundary condition is imposed on the top surface. To capture the local flow in the boundary layer around the rod bundles, fine mesh layers are generated as a hexagonal surface mesh, and this 2D mesh in the region of the tube pitch is extruded to obtain the 3D mesh. At the wall, y
+ value in (3) is checked to y
+
<
30. Half of the test section is modeled by applying the symmetric boundary condition on the upper plane of computational domain to enhance the efficiency of computation with a reduced number of mesh elements. A series of convergence test for grid has been done to make sure that the number of grids should be sufficient for the precise solution. Finally, the total number of hexagonal elements is estimated as 424,320. The grid model in Figure 9 is used commonly to 3D computations using ANSYS-CFX and OpenFOAM. In the reduced 2D computation, the grid in the depth direction is merged as a cell for OpenFOAM.Figure 9
Grid system for the benchmark validation.As shown in Teyssedou et al. [10], the 2D computation tends to underestimate the pressure drop, and all codes using k
-
ε turbulence model presented the same trend and similar order of qualitative values: see Table 3. The hydraulic diameters in the 2D model are far different from the original 3D experimental setup, so flow resistance in the depth direction cannot be neglected, generating additional pressure drop. The velocity fields visualized in Figure 10 show that the 2D flow in Figure 10(a) is different from the 3D flow in Figures 10(b) and 10(c), characterized as the irregularity of velocity vectors in the depth direction.Table 3
Comparison of pressure drop for the STERN laboratory experiment, 2D.
Case
Δ
p (Pa)
Experiment(STERN)
COMSOL(2D)
ANSYS-FLUENT (2D)∗
OpenFOAM(2D)
#1
28.2
15.0
16.0
14.0
#2
41.3
22.5
20.3
21.4
#3
78.2
45.2
36.4
38.0
∗Quoted for the data with k
-
ε turbulence model in FLUENT version 15 [10].Figure 10
Velocity fields from the computation with various codes, Case  3: (a) COMSOL (2D), (b) ANSYS-CFX (3D), and (c) OpenFOAM (3D).
(a)
(b)
(c)In the 3D computation, Table4, the error between computational and experimental value is far reduced, but there is still the discord among the experimental result and computations. The cause of error might be originated from the limitation of turbulence model because the flow regime lies in the low Reynolds number region from 2,709 to 9,308. The Kármán instability of the stack of cylinders makes the numerical solution difficult to converge for the case of 2D computation where often the pressure drop is underestimated [10] especially in the low Reynolds number region, but this can be far more relaxed with the multidimensional turbulence effect in 3D results because of the increase in degrees of freedom. Case  3 flow field is compared for ANSYS-CFX and OpenFOAM in Figures 10(b) and 10(c), equivalently.Table 4
Comparison of pressure drop for the STERN laboratory experiment, 3D.
Case
Δ
p (Pa)
Experiment(STERN)
Equation (7)∗
COMSOL(3-D)
ANSYS-CFX (3-D)
OpenFOAM(3-D)
#1
28.2
27.0
22.3
23.6
23.6
#2
41.3
40.1
36.3
37.6
38.3
#3
78.2
77.7
76.9
77.5
79.9
∗The porous model is used with the experimental correlation [15].Hadaller et al. [15] developed the experimental pressure drop correlations of this problem for aligned and staggered arrangement cylinders. In the aligned or in-line correlation, the pressure drop per unit channel length is(7)
Δ
p
Δ
L
=
4.54
R
e
d
-
0.172
l
cos
⁡
θ
1
2
ρ
V
m
2
;
R
e
d
=
ρ
V
m
d
μ
,where l
=
71.4 mm is the pitch of cylinders and θ is the incident angle of free stream, zero for this problem. The free-stream flow velocity is defined as V
m
=
ε
p
V, and the velocity magnitude V
=
∑
i
V
i
2 is measured inside the tube bank [2]. The porosity in this problem is defined as [14](8)
ε
p
=
1
-
π
/
4
d
2
l
2
.In the comparison of pressure drops,Δ
p between PT1 and PT3 in Figure 8(b) from all the codes is given in Table 4. ANSYS-CFX and OpenFOAM data are from the direct simulation of (1)~(2) while (7) is used for the corrected body force per unit volume in the momentum equation (2): that is, f
V
,
x
=
-
Δ
p
/
Δ
L in the source term of the COMSOL code [14]. The correlation in (7) shows a good agreement with the experimental data, but the numerical data also result in comparable agreement. Figure 11 is the relative errors from the experimental data and shows the full performance of codes for the present benchmark problem. All the computational results except for the OpenFOAM of Case  3 predict underestimating of the experimental data. Some artificial factors like surface roughness or incident turbulence could have effect on the experimental date for the increase of pressure drop.Figure 11
Comparison of relative errors for various 3D codes from the experimental data.However, in spite of the uncertainty of experiment, Case  1 shows the worst prediction of numerical codes (20.9% for COMOSOL; 16.3% for both ANSYS-CFX and OpenFOAM), which is strongly suspected to the defect of turbulence model at low Reynolds number as 2,709. The laminar computation with no turbulence model is never converged to the steady solution but oscillates within 100 times of the tolerance of error. For Case  1, the turbulent solution is converged att
=
1,451 s with OpenFOAM. The pressure drop between PT1 and PT3 in Figure 8(b) is calculated at t
=
10,000 s for laminar flow where the predicted value is Δ
p
=
24.6 Pa, and the error from experimental data is mitigated to 12.8%. The wake flow fields computed from laminar and turbulent model are compared for Case  1 streamlines in Figures 12(a) and 12(b). In the laminar model, as the shear layer bounding on the separated flow from cylinders blocks the flow along the centerline of symmetry, that is, the upper plane, the flow resistance makes the horizontal pressure gradient slightly increase more than that of the turbulent model.Figure 12
Comparison of wake flow field in the streamlines for Case  1, OpenFOAM with (a) laminar, or no turbulence model, and (b)k
-
ε turbulence model.
(a)
(b)


## 4.3. Sensitivity to Turbulence Models

Various turbulence models are tested for the same benchmark problems in Section4.2. 2D models with COMSOL and 3D models with ANSYS-CFX are presented in Tables 5 and 6, respectively. Four models are tested such as standard k
-
ε, k
-
ω, SST (shear stress transport) [16], and SA (Spalart-Allmaras) models [17]. As shown overall, the k
-
ε model shows the most reliable agreement with the experimental data: see Tables 5 and 6. Reference [10] also shows that the error of standard model is less than those of various k
-
ε models: R, RNG, and so forth, although they are all 2D computations.Table 5
Sensitivity to the turbulence models, 2D (COMSOL).
Case
Δ
p (Pa)
k
-
ε
k
-
ω
SST
Spalart-Allmaras
#1
15.0
12.3
9.5
11.4
#2
22.5
17.3
15.1
—∗
#3
45.2
33.0
31.0
—∗
∗They do not converge under the present grids set.Table 6
Sensitivity to the turbulence models, 3D (ANSYS-CFX).
Case
Δ
p (Pa)
k
-
ε
k
-
ω
SST
Spalart-Allmaras
#1
23.6
21.8
20.6
26.7
#2
37.6
34.7
31.8
42.1
#3
77.5
74.8
73.9
115.2In the agreement of pressure drop in the STERN laboratory test, thek
-
ε model seems better than the k
-
ω model for both 2D and 3D, which is thought to the fact that they are using different wall functions for “the enhanced wall treatment” inside the inner boundary layer (y
+
≈
1) [10]. The SST model seems not a good solution for this problem because it is originally based on the k
-
ω model. In the 3D low Reynolds number case (Case  1), the error from the experimental data is decreased to 5.3% with the one-equation Spalat-Allmaras model in Table 6. However, it is not accurate for the case of the highest Reynolds number (Case  3) as shown in Table 6 and even it does not converge to a steady solution (Case  2 and 3) for 2D in Table 5, where the present grids system is not adequate or too coarse for the one-equation turbulence model.


## 5. Concluding Remark

With various CFD codes including an open-source code, which are COMSOL Multiphysics, ANSYS-CFX, and OpenFOAM, this research has shown the feasibility of application to the simulation of thermal hydraulics in the moderator system of a PHWR. The governing equations are Navier-Stokes equations with turbulence model where, in certain cases, some correction is required for the physics such as porous media model for the calandria tube stacks, heat source distribution, and Boussinesq approximation for thermal flow with temperature field. The computational environment for each code has been set up independently.2D and 3D application using the same level of grid system is tested for various codes on the STERN laboratory experiment as a benchmark problem. The 2D computation cannot predict the exact pressure drop in the experiment. The 3D results from the present codes are quantitatively compared, including the porous model with Hadaller’s correlation, with each other for the frictional loss of internal flow in the moderator tank of CANDU-6 reactor. Case  1, at the lowest Reynolds number, results in 16.3 to 20.9% error from the experimental data, but the flow regime lies in the instability region where the turbulent model should be defected. The laminar computation for Case  1 shows that the error can be decreased to 12.8%. The sensitivities for the change of turbulence models are traced for 2D and 3D computations, and thek
-
ε model gives the best result of agreement with experimental date for this given grids system. For Case  1 (the lowest Reynolds number), the Spalart-Allmaras model improves the result to 5.3% error, but the one-equation turbulence model should be very carefully implemented because of its computational instability. The results of commercial and open-source codes lie within an acceptable agreement for the benchmark problem, revealing the feasibility of 3D CFD for the solution of hydrothermal problems applying to PHWR, under the conditions for not only normal operation but also serious accident situations in potential.



---

*Source: 3194839-2016-09-05.xml*