|
|
TWO-DIMENSIONAL CALCULATION
OF FLUID FLOW AND HEAT TRANSFER IN A ROD BUNDLE WITH GEOMETRICAL DISTURBANCE
Vladimir Kriventsev and
Hisashi Ninokata
Research Laboratory for Nuclear Reactors, Tokyo
Institute of Technology
2-12-1 O-okayama, Meguro-ku, Tokyo 152, JAPAN
Phone: +81-3-5734-3056, Fax: +81-3-5734-2959
e-mail: hninokat@nr.titech.ac.jp
Abstract
The results of numerical simulation of fluid
flow and heat transfer in rod bundles with geometrical disturbance to the
hexagonal rod array configuration are presented. The geometry of the rod bundle
was chosen according to the flow and temperature distributions available from
the experimental data.
Reynolds equation for axial velocity component
were simulated in two dimensions. Turbulent shear stresses were modeled by the
turbulent eddy viscosity with anisotropy defined for radial and azimuthal
components. Secondary flows were not taken into consideration. Also, the
averaged energy conservation equation closed with anisotropic turbulent
conductivity coefficients was simulated.
Reynolds and energy conservation equations were
discretizated by the Efficient Finite-Difference (EFD) scheme based on the
"locally exact" analytical solution. A comparison of the accuracy of
the EFD method and traditional central-difference and upwind schemes were
performed.
The benchmark problem was simulated using
components of the Computational Object-Oriented Library for Fluid Dynamics
(COOLFD) which is a new-generation programming tool aimed to improve the
development of the CFD application to complex calculation areas such as rod
bundles in nuclear reactors.
Comparisons of calculated results and
experimental data are presented for the local shear stress, axial velocity and
the wall temperature distributions in the "geometrically disturbed"
region around the dislocated rod.
Introduction: Problem
Definition
The geometry of the rod bundle was chosen
according to the benchmark problem for the 9th IAHR Working Group
Meeting (April 7-9, 1998, Grenoble, France). For such a case, experimental data
for local velocity and wall shear stress distributions were obtained by Hejna
and Mantlik (1998) at NRI (Czech Republic). Another series of the experiments
which provide data on the wall temperature profiles were done by Jukov et al. (1985)
at IPPE (Russia). Both experiments provide complete sets of data for comparison
with the results of numerical simulation.
Governing Equations
The Reynolds equation for steady-state fully
developed turbulent incompressible flows in the two-dimensional Cartesian coordinate
system may be written for the axial velocity component as follows
|
|
where W is the
axial velocity component; v
is the viscosity; 
and 
are the components of
anisotropic turbulent eddy diffusivity; r is the density; P is the
pressure; x1 and
x2 are the
dimensions of the orthogonal coordinate system; and z is the axial coordinate. Here, W is the
unknown dependent variable. Boundary conditions were applied as a non-slip
condition on the walls, zero shear stresses on the outer symmetry lines and as
a constant mean velocity given by conditions of the experiment. The pressure
gradient in the axial direction (a constant in the case of fully-developed flow)
depends on the wall shear stress distribution and, therefore, should satisfy a
flow distribution with the mean velocity above.
Also, an equation of the energy conservation for steady-state fully developed turbulent incompressible flows in two-dimensions was written as follows:
|
|
where T is the
temperature of the fluid; l is the thermal conductivity; 
and 
are the components of anisotropic turbulent eddy thermal
conductivity which are defined using the turbulent Peclet number PeT
as
|
|
(3) |
As was mentioned, secondary flows were not considered in this analysis. The reason is not only the additional complexity of numerical simulation of Navies-Stocks equations for all three velocity components but also uncertainty with the definition of turbulent eddy diffusivity. These coefficients are derived to fit experimental velocity distributions and for most data available it is not clear whether turbulent eddy diffusivity includes the effect of secondary flows or not. It should be noted here that, for those rod bundle experiments where the only axial velocity component is measured, it is impossible (in principal) to calculate true turbulent viscosity components.
Finite-Difference Discretization
A finite-difference method based on the "locally exact" control volume scheme was applied in this work. The detailed description of the Efficient Finite-Differencing scheme was given by Kriventsev and Ninokata (1997). In this section, let us consider an application of the EFD discretization to the governing equations (1) and (2). The main idea of the EFD scheme is based on the idea of using an exact analytical solution of the simplified one dimensional convection-diffusion equation with the finite-difference estimation of the interface flux of the transported quantity (shear stress in the case of Eq. (1) and heat flux in Eq. (2)). This simplified equation can be written for each direction for every control volume. All other terms of the original equation including transport terms from other directions and a transient term are collected in an extra-source term. In doing so, with EFD, we assume that both the extra source term and diffusion coefficient (turbulent diffusivity in our case) are distributed linearly within the neighboring nodes. For example, for the first direction, this simplified equation can be written as follows:
|
|
where 
; 
.
Taking into account the values of the transported quantity at i and i+1 points as boundary conditions we can solve Eq. (4) analytically. It allows us to derive the expression for Reynolds shear stress at the interface between i and i+1 control volumes (the following is for the simplified case when the source term is assumed to be a constant):
|
|
(5) |
while regular upwind and central-difference schemes give the following disccretization:
|
|
(6) |
Expression
(5) was used in this work for finite-difference discretization, as well as for estimation of shear stresses at the rod walls. The EFD scheme is also applicable in a similar way not only in a Cartezian coordinate system, but also in an arbitrary orthogonal coordinate system. Also, similar discretization can be written for energy conservation equation (2) and numerical estimation of heat fluxes on the wall of the rod bundle.
Mesh System
In calculation of this problem, we used an orthogonal coordinate transformation. This transformation was performed using the grid generation components of the Computational Object-Oriented Library for Fluid Dynamics (COOLFD) as described in Kriventsev and Ninokata (1998). First, elementary mesh systems were generated for each typical rod as shown at the Fig.1. All of these elementary meshes were assembled together to fit the whole rod bundle. In doing so, the visual components of the COOLFD library were used as shown at Fig.2.
Eddy diffusivity coefficients
The coefficient of turbulent viscosity in the radial direction was applied from Nijsing and Eifler (1971) as follows:
|
where |
(7) |
The expression of Neelen (1986) was used for eddy diffusivity in the azimuthal direction:
|
where and |
(8) |
In spite of the fact that relations above were obtained for radial and azimuthal directions relatively to the rod wall, we used them directly in the case of the orthogonal system for v1 and v2 correspondingly. The grid lines of the orthogonal meshing system practically coincide with a polar coordinate system near the wall, but for the region in the center of subchannel, the radial and azimuthal direction cannot be defined clearly. It may result in an additional error in numerical solution.
Numerical Procedure
The computational part of the COOLFD library was used to calculate a numerical solution. The details of COOLFD were described in Kriventsev and Ninokata (1998). The Gauss-Zeidel method was used to calculate an improved estimation in internal iteration loop. Next, wall shear velocities u* were recalculated for wall points and then, all the values of turbulent eddy diffusion and thermal conductivity were recalculated for the entire region. An outer iteration loop was used to adjust the pressure gradient which should produce the flow distribution with mean velocity as given by experimental conditions (<W> = 48.34m/s).
The calculation of the flow distribution in the 17-pin rod bundle divided into the 8808 control volumes took about 3min of CPU time on a personal computer with an Intel Pentium-II 400Mhz processor (including mesh generation). In addition, the calculation of the temperature distribution with the energy conservation equation demands about one minute more. A sample of the visual object that serves to calculate the region is shown at Fig.3.
Calculation Results
The results of numerical simulation of fluid flow were compared with experimental data presented by F. Mantlic et. al. at the 9th IAHR Working Group Meeting (April 7-9, 1998, Grenoble, France). Air was used as a working fluid in these experiments. Calculated wall temperature distribution was compared with the liquid metal experiment performed in IPPE by Jukov et al.(1985).
The local shear stress distributions at the wall of rods No 6 and 7 are shown in
Fig.4 and
Fig.5. Here, one can see that the traditional upwind scheme underestimates the shear stresses significantly while the EFD scheme generates fairly satisfactory results (taking into account the complexity of the problem itself as well as insufficient data on eddy diffusivity for rod bundles with geometrical disturbance). The plots in Figs. 6 and 7 give samples of the velocity profiles for rods 6 and 7.
The comparison of calculated temperature distribution around the central dislocated rod with experimental data using sodium as a working fluid is shown in Figs.8 and 9 for Peclet numbers 218.0 and 58.0 respectively.
Conclusion
The results of numerical calculation show a satisfactory agreement with experimental data both for flow and temperature distribution. It proves the applicability of the COOLFD software and EFD scheme for numerical simulation of fluid flows and heat transfer in complex regions, such as a rod bundle of a nuclear reactor. Some additional work should be done to improve the calculation procedure. First of all, a more accurate model of the turbulence should be implemented. Also, the authors believe that three-dimensional simulation, taking into account secondary flows, could improve the numerical accuracy significantly.
References
1. J. Hejna, F. Mantlik, "Benchmark Problem Description", 9th Meeting of IAHR Working Group on Advanced Nuclear Reactor Thermal Hydraulics, CEA-Grenoble, France, 1998
2. A.V. Jukov, P.L. Kirillov, N.M. Matuhin, A.P. Sorokin et al., Thermal Hydraulics Analysis of Fuel Assemblies of Fast Reactor cooled by liquid metal, Moscow, Energoatomizdat, 1985 (in Russian)
3. V. Kriventsev and H. Ninokata, "An Effective Finite-Difference Scheme For Transient Problems of Heat and Mass Transfer", Proc.of 8th Intl. Topical Meeting on Nuclear Reactor Thermal-Hydraulics (Nureth8), Kyoto, (1997), Vol.1, pp.572-580
4. V. Kriventsev and H. Ninokata, "An Application of Object-Oriented Programming in Numerical Simulation of Fluid Dynamics in Complex Geometry", Proc. of Spring Meeting of the JSME, Tokyo, March-April (1998)
|
|
|
Fig.1 Elementary grid components |
Fig.2 Elementary grid components assembled in the rod bundle |
|
|
|
Fig.3 Sample of Visual Component showing the flow distribution at the
experimental 17-pin rod bundle with dislocated central rod |
|
t[Pa] |
t[Pa] |
|
|
|
|
Fig.4 Wall shear stress t; Rod No 6 |
|
|
|
|
Fig.6 Axial Velocity; Rod No 6 |
Fig.7 Axial Velocity; Rod No 7 |
|
|
|
|
Fig.8 Wall Temperature; Rod No 1; Pe = 218 |
Fig.9 Wall Temperature Rod No 1; Pe = 58 |
,
,
,
,
,
.
,
.
,
