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.
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 (Figure 1), 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.
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:where , , and are velocity vector, density, and pressure while the constants and 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, , 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, to , 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.
Standard 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 , 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.
Figure 1 is the benchmark result for a cylinder-bank problem with in-line 24 rows where the diameter of each tube is .8 mm, and the pitch is 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, , where and 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, .
(a)
(b)
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].
(a)
(b)
value is defined for , and the normal length of grid at the first one neighboring the wall is where 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 as in whole computational domain. However, the required number of grids can be much reduced with use of wall functions [7]:where is the turbulence kinetic energy and is a coefficient. Therefore, the grids in the range of are neglected to be unnecessary, and less than 30 is sufficient for the capture of outer layer in a low Reynolds number flow for engineering computation. Therefore, 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 command paraFoam [8].
3.3.1. 3D Grid Generation
There are two kinds of methods for the grid generation in OpenFOAM. With blockMesh 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 command fluentMeshToFoam 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.
(a)
(b)
3.3.2. Computational System
With the simpleFoam solver, the 3D moderator model in Figure 5 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 the polyMesh 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.
(a)
(b)
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. Table 1 is the list of solvers used in the present research.
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 is (pressure) and (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 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.
(a)
(b)
The Reynolds number in Table 2 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 model aswhere 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
Figure 9 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, value in (3) is checked to . 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.
As shown in Teyssedou et al. [10], the 2D computation tends to underestimate the pressure drop, and all codes using 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.
(a)
(b)
(c)
In the 3D computation, Table 4, 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.
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 iswhere 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 , and the velocity magnitude is measured inside the tube bank [2]. The porosity in this problem is defined as [14]
In the comparison of pressure drops, 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, 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.
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 at s with OpenFOAM. The pressure drop between PT1 and PT3 in Figure 8(b) is calculated at s for laminar flow where the predicted value is 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.
(a)
(b)
4.3. Sensitivity to Turbulence Models
Various turbulence models are tested for the same benchmark problems in Section 4.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 , , SST (shear stress transport) [16], and SA (Spalart-Allmaras) models [17]. As shown overall, the 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 models: R, RNG, and so forth, although they are all 2D computations.
In the agreement of pressure drop in the STERN laboratory test, the model seems better than the 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 () [10]. The SST model seems not a good solution for this problem because it is originally based on the 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 the 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.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
This work was supported by (1) the National Research Foundation of Korea (NRF) grant funded by the Korea government (Ministry of Science, ICT, and Future Planning) (no. NRF-2012M2A8A4025964) and (2) the Human Resources Development Program (Grant no. 20144030200590) of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) grants funded by the Korea government (Ministry of Trade, Industry, and Energy).