- Email: [email protected]
Computer simulation of natural convection heat transfer from an assembly of vertical cylinders of PARR-2
Applied Thermal Engineering 27 (2007) 194–201 www.elsevier.com/locate/apthermeng
Computer simulation of natural convection heat transfer from an assembly of vertical cylinders of PARR-2 M. Abdul Basit a, Muhammad Rafique a
b,*
, Imran R. Chughtai a, Mansoor H. Inayat
a
Department of Chemical and Material Engineering, Pakistan Institute of Engineering and Applied Science Nilore, Islamabad 45650, Pakistan b Department of Physics, University of Azad Jammu and Kashmir Muzaffarabad 13100, Azad Kashmir, Pakistan Received 11 August 2005; accepted 1 May 2006 Available online 2 August 2006
Abstract In this paper a Computer Code COSINAC (Computer Simulation of Natural Convection from Assembly of vertical Cylinders) has been developed to simulate the natural convection heat transfer from an assembly of vertical cylinders of Pakistan Research Reactor-2 (PARR-2), under the steady state reactor operation. The momentum and energy equations in cylindrical co-ordinates, representing the thermal hydraulic behavior of a typical fuel pin in Pakistan Research Reactor-2, have been solved numerically for a two dimensional axisymmetric domain. The temperature and velocity profiles and Nusselt number variations have been studied and results have been presented. The computer code COSINAC has been validated against experimental results carried out in previous studies at different occasions. Average outlet coolant temperature simulated by computer code, at different wall heat fluxes, has been found in good agreement with experimental results. 2006 Elsevier Ltd. All rights reserved. Keywords: Natural convection; Pakistan research reactor-2; Assembly of vertical cylinders; Nusselt number
1. Introduction The nuclear reactors are classified on various concepts and basis. These may be classified on the basis of coolant or moderator used, type, the enrichment employed and arrangement of fuel. Most common classification is on the basis of application for which reactor is being used i.e. whether the reactor is meant for fuel production, power generation or research purposes. Whatever may be the type of reactor, an adequate cooling mechanism is a major requirement for safe operation of the reactor. Depending upon the quantity of heat generated by a reactor, different coolants and cooling mechanisms may be utilized. Power reactors are cooled by forced convection while research reactors, with small amount of heat generation, may
*
Corresponding author. Tel.: +92 03009189301/+92 05881034691; fax: +92 51 9223727. E-mail address: rafi[email protected] (M. Rafique). 1359-4311/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2006.05.001
employ free convection as the cooling mechanism. PARR-2, the reactor considered in the present study, consists of assembly of vertical cylinders cooled by natural convection. Extensive study on the natural convection cooling of thin vertical cylinders in various configurations has been done in the past. Steady state natural convection from a single cylinder with known heat flux on the wall was studied experimentally by Nagendra et al. [1] and numerically by Heckel et al. [2]. Valusami and Garge [3] have considered transient heat transfer from a vertical cylinder for a step change in surface heat flux. Keyhani et al. [4] studied steady state natural convection in a vertical annulus enclosed from both ends while Chughtai and Inayat [5] studied the same phenomena for a vertical annulus open from both ends. Keyhani et al. [6] have reported steady state natural convection from an assembly of vertical cylinders enclosed in a canister. For the reactor under consideration, simplified lumpedparameter analyses have been performed to model its
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
195
Nomenclature cp g k L Nuz
specific heat of fluid at constant pressure acceleration due to gravity thermal conductivity of fluid length of the assembly local Nusselt number for a typical fuel pin in assembly of vertical cylinders NuL local Nusselt number for a typical fuel pin in assembly of vertical cylinders at exit NuzCylinder local Nusselt number for vertical cylinder, of same dimensions as fuel pin, in an infinite medium (calculated within the code) q00 wall heat flux r radial position c qgbz4 q00 Ra* modified Rayleigh number = p k2 m2
thermal hydraulics [7–9]. Khan et al. [9] studied the natural convection cooling of PARR-2 in 1993. They used mathematical models to calculate coolant velocity as function of reactor power. The coolant velocity was also determined experimentally by operating the reactor at different power levels. A comparison between the calculated values and experimental values was also done. In the present study, an attempt has been made to study the thermal hydraulic phenomenon around a typical fuel pin of PARR-2. To achieve this end goal a single channel model has been used. An annular region filled with coolant, around a typical fuel pin in the reactor, has been assumed. Heat generation takes place in the fuel rod and heat is then transferred to coolant flowing parallel to the axis of the fuel pin in the annular region. The coolant enters the annular region from below and leaves from top. The thermal boundary-layer equations have been solved for the fluid flowing through the annular region around a typical fuel pin in the reactor, in cylindrical co-ordinates, using finite volume method. Density of fluid has been allowed to vary according to Boussinesq approximation and all other fluid properties have been assumed constant. Calculations have been made for velocity and temperature profiles and Nusselt numbers. 1.1. Reactor description PARR-2 is a tank-in-pool type reactor with a nominal power rating of 27 kW. It is composed of highly enriched uranium (HEU) fuel and is moderated and cooled by light water through natural convection. The HEU fuel is fabricated as thin vertical cylindrical pins of uranium–aluminum alloy with aluminum cladding. Each fuel pin has a heat generating length of 230 mm and an inert portion of 10 mm on both ends. Overall diameter of each fuel pin, including cladding thickness, is 5.9 mm. The fuel pins are arranged in concentric arrays with constant pitch of 10.95 mm so as to make a cylindrical shaped reactor core.
rm
r0 T T1 u v z a b m q
radius of the circle circumscribed about a fuel rod representing the flow area and fluid assignable to that particular fuel rod radius of fuel pin including cladding temperature bulk temperature of the fluid component of velocity in axial direction component of velocity in radial direction axial position thermal diffusivity of fluid volumetric coefficient of thermal expansion kinematics viscosity of fluid density of fluid
The reactor core is surrounded by radial as well as axial beryllium reflectors. There is a central control rod in the reactor core that is used to control the neutron flux and hence the reactor power. Energy generation rate, in the reactor, may be increased by withdrawing the control rod. However the control rod can be withdrawn only at the fixed rate and hence the reactivity insertion through control rod can only be a ramp function of time. The coolant flow path is restricted by placing two orifice plates on the top and bottom of the core (Fig. 1). The restricted coolant flow adds the self power-limiting feature to this reactor. The coolant enters the reactor core through the gap between the lower orifice plate and the annular reflector, rises up due to buoyancy and leaves through the upper orifice plate while cooling the fuel pins. Two k-type thermocouples have been installed near the lower and upper orifice plates to measure the inlet and outlet temperature of the coolant. A calibrated fission chamber type neutron detector has also been installed to measure the neutron flux and hence the reactor power. Part of the coolant coming out of the core is re-entrained into the in-going stream of the coolant and rest rises up into the tank due to buoyancy. The heat rejected from the core is first transferred to the water in the tank and then to the outer pool by natural convection. A brief summary of the reactor specification is given in Table 1. 2. Mathematical model 2.1. Physical and mathematical modeling The actual problem of heat transfer from PARR-2 has a complex 3-dimensional nature; however, due to axisymmetric geometry and physical phenomena, the problem can be studied efficiently in a two dimensional axisymmetric cylindrical coordinate system. It is the case of natural convection heat transfer, in axisymmetric geometry; hence flow is taking place due to buoyancy force caused
196
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
Fig. 1. Schematic diagram and coolant flow patterns in PARR-2 (not to scale).
Table 1 Summary of design features of PARR-2 Item
Value
Rated steady state power Enrichment in U-235 Maximum excess reactivity Thermal neutron flux at 27 kW Moderator and coolant Reflector
27 kW 90% 4.0 mk 1.0 · 1012 #/cm2 s Light water Beryllium
Fuel: Geometry Meal material Clad material Number
Cylindrical UAl4 – Al 303-1 Al 344
Control rod: Material Number
Cadmium 1
Detectors: Neutron Temperature
Fission chamber k-Type thermocouples
by constant heat flux applied on a fuel pin surface. As surface heat flux is same on all the fuel pins and a triangular pitch arrangement has been considered, so flow conditions may be taken same on all fuel rods and flow will be symmetric. The core of PARR-2 cannot be regarded as triangular or square pitch arrangement as shown in Fig. 2, but considering the stagger in fuel pin arrangement, a triangular pitch arrangement seems to be more appropriate. Following a procedure used by Friedland and Bonilla [10] for a tube bundle with an equilateral triangular pitch, a hexagon can be circumscribed about a tube to represent the flow area, and fluid assignable to that particular tube. It is assumed that this situation can be approximated closely enough by replacing the hexagon with a circle of equal area [11] having a radius rm as illustrated in Fig. 3. Coolant
flows through the annulus between fuel pin and this circular region. The dimensions of this annulus have been calculated by dividing the total cross sectional area on 354 fuel pins in PARR 2. Radius of inner cylindrical solid region is 2.95 mm (Radius of fuel pin). Around this solid region there is 3.22 mm thick annular region for coolant flow. So the total radial dimension is 6.17 mm. Consider an assembly of vertical cylinders (fuel rods) aligned in a quiescent (dormant) ambient fluid at temperature T1. Each cylinder in this assembly has a radius of r0. The z-coordinate is taken along axial direction from the leading edge of cylinder and r coordinate in radial direction from the axis of cylinder. The surface of the cylinder is subjected to a constant heat flux q00 . The acceleration due to gravity g is acting downward. All fluid properties are assumed to be constant except for density variation which induces the buoyancy force. Assumption of incompressible laminar boundary-layer flow and the Boussinesq approximation have been employed. The fluid velocities are quite small; therefore the compressive work and viscous dissipation terms in the energy balance equation can be neglected. The three governing conservation equations in the cylindrical co-ordinates, for natural convection, can be written as Continuity equation o o ðruÞ þ ðrvÞ ¼ 0 oz or
ð1Þ
z-Momentum equation u
ou ou t o ou þv ¼ r þ gbðT T 1 Þ oz or r or or
Energy balance equation oT oT a o oT þv ¼ r u oz or r or or
ð2Þ
ð3Þ
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
197
Fig. 2. Arrangement of fuel pins in the PARR-2 core.
Outflow Boundary Condition
Heated Cylinders Hexagonal Area Around Cylinder 6.17
2.95 11.83
φ5.50 Equivalent Unit Cell with Circular X-section Symmetry Boundary
R6.17
Wall Boundary (Constant Heat Flux) 230 Axis of cylindrical domain
z g
Hexagonal flow domain around a typical fuel pin r
Constant Pressure Boundary (T=293 K) 2-D Axisymmetric model used for modeling. (Not to Scale)
Fig. 3. Schematic diagram of the axisymmetric model for a triangular pitch assembly of vertical cylinders (all dimensions are in millimeters).
The boundary conditions for the solution of Eqs. (1)–(3) can be expressed, at the edge of the boundary-layer, as follows for constant heat flux:
For z = 0 and r P r0T = T1 00 For r = r0 and z > 0 u = 0, v = 0, oT ¼ qk or w For r = rm and z P 0 ou ¼ 0, v = 0, oT ¼0 or or
198
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
At z = L, outflow boundary condition [12] has been employed.
330 325
2.2. Method of solution
40x40 Grid points 20x20 Grid points 10x10 Grid points
320
overall trend observed in the domain. It can be seen that the three mesh sizes have produced almost identical results.
3. Discussion
3.2. Temperature and flow patterns
3.1. Grid independence
Fig. 5 shows velocity and temperature profiles for surface heat flux of 19,800 W/m2. In the case of assembly of vertical cylinders, flow is entrained only in axial direction. On one side, axial velocity is zero at the cylinder wall; on the other hand the outer side of the axisymmetric domain is the middle plane between two consecutive heated cylinders. Hence, resistance to flow of coolant is posed only at one side.
Three mesh sizes of 10 · 10, 20 · 20 and 40 · 40 grid points have been selected to study the grid independence of the code results. A uniform mesh size has been considered. Fig. 4 shows temperature profiles in the coolant for surface heat flux of 18,900 W/m2 at an axial distance of 0.158 m from the lower end. These profiles show a typical
Temperature (K)
Eqs. (1)–(3) constitute a system of coupled, non-linear, partial differential equations in the (r, z) co-ordinates. In the present study, a typical semi-implicit finite volume solution was carried out. The details of this numerical solution method are omitted to save space. Some of its highlights are being given. The first step in solving Eqs. (1)–(3) is to convert partial differential equations (2) and (3) given in non-conservation form into conservation form of partial differential equations. These conservation form of equations are then cast into algebraic equations using finite volume method employing Power Law scheme discussed by Patankar [12]. Axial velocity u has been calculated from discretized z-momentum equation (2) employing pointto-point Gauss–Siedel method. Radial velocity v has been calculated from discretized continuity equation. Finally, temperature distribution has been obtained by solving discretized energy equation using point-to-point Gauss–Siedel method.
290 0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
Radial Distance (mm) Fig. 4. Comparison of temperature profiles for natural convection from the axisymmetric model of an assembly of vertical cylinders using three mesh sizes (heat flux on hot surface = 19,800 W/m2, distance from the lower end = 0.158 m).
320
z/L=1.0
310
z/L=1.0
0.01
300
0.00 0.03
320
0.02 0.01 0.00 0.03 0.02
z/L=0.6
0.01 0.00 0.03 0.02
z/L=0.4
0.01
z/L=0.8
310
z/L=0.8
Temperature (K)
Axial velocity (m/s)
305
295
0.02
0.02
310
300
0.03
0.00 0.03
315
300 320
z/L=0.6
310 300 320 z/L=0.4
310 300
Surface of Cylinder
320 z/L=0.2
0.01
z/L=0.2
310 300
0.00 0.003
0.004
0.005
Radial distance (m)
0.006
0.0025 0.0030 0.0035 0.0040 0.0045 0.0050 0.0055 0.0060 0.0065
Radial distance (m)
Fig. 5. Profiles of axial velocity and temperature at various axial positions for natural convection in an axisymmetric model of an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
4
3
Nu z /Nu z Cylinder
Although, along the axial direction, shape of velocity profile changes according to the local buoyancy effects, yet area under a typical velocity profile almost remains unchanged. As the wall temperature is low in the lower portion, velocity profile has maxima at middle plane between two heated cylinders. But as the fluid moves axially along the cylinder, maxima of the velocity profile shifts inwards gradually due to increased buoyancy caused by increase in the surface temperature.
199
2
3.3. Flow entrainment In the case of assembly of vertical cylinders, fluid enters the domain from lower side only hence the entrainment remains almost unchanged with axial distance along the cylinder. Small slope of the line in Fig. 6 is due to temperature dependence of coefficient of volumetric expansion of the coolant.
1
0.00
0.05
0.10
0.15
0.20
Axial distance from the lower end (m) Fig. 7. The ratio of local Nusselt numbers Nuz/NuzCylinder for an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
3.4. Nusselt numbers 1.20 1.18 1.16 1.14
Nu L /Nu L,Cylinder
Fig. 7 shows a plot of Nuz/NuzCylinder versus axial distance along the surface of fuel rod. A vertical cylinder in infinite medium can entrain flow axially as well as radially, whereas a vertical cylinder in an assembly entrains fluid from lower end only, so flow rate in the case of assembly of vertical cylinders is higher in the lower portion as compared with cylinder in infinite medium. The phenomenon of fluid entrainment in assembly of vertical cylinders can be understood by comparing the situation to what happens in a chimney, where fluid rises as it becomes hot and surrounding cooler fluid rushes to the place to occupy the space. This is so called chimney effect. From Fig. 7, it can be seen that Nuz/NuzCylinder is well above 1.0 in lower region but in upper portion Nuz almost approaches to NuzCylinder.
1.12 1.10 1.08 1.06 1.04 1.02 1.00 1011
1012
1013
Modified rayleigh No. Ra* L Fig. 8. The ratio of overall Nusselt numbers Nuz/NuzCylinder for an assembly of vertical cylinders.
6.0x10-6
Flow entrainment (m3 /s)
5.0x10-6
The ratio of overall Nusselt number attained by a cylinder in an assembly to that of a cylinder in an infinite medium has been plotted in Fig. 8. It can be seen that due to so-called chimney effect, NuL for assembly is higher than that for an isolated cylinder for the whole range of RaL . Fig. 9 shows Nuz vs. Raz relation for a cylinder in assembly. It can be observed that Nuz vs. Raz data falls on a single curve for all values of surface heat flux studied in the present research. It is a well known characteristic of natural convection heat transfer that local Nusselt number is a function of modified Rayleigh number only.
4.0x10-6
3.0x10-6
2.0x10-6
1.0x10-6 0.00
0.05
0.10
0.15
0.20
0.25
Axial distance from the lower end (m) Fig. 6. Volumetric flow entrainment as a function of axial distance for an axisymmetric model of an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
3.5. Surface and maximum fuel temperatures Fig. 10 shows temperature variation along fuel pin surface versus axial distance along the surface of cylinder in
200
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
Fuel maximum temperature (K)
10 2
Nu z
Heat flux=3600W/m2 Heat flux=7200W/m 2
10 1
Heat flux=10800W/m2 Heat flux=14400W/m2 Heat flux=18000W/m2 Heat flux=19800W/m2
10 0 103
104
105
106
107
108
109
1010
1011
1012
1013
380 375 370 365 360 355 350 345 340 335 330 325 320 315 310 305 300 295 0.00
Ra*z
345
Heat flux=3600W/m2 Heat flux=7200W/m 2 Heat flux=10800W/m2
340
Heat flux=14400W/m 2
335
Heat flux=18000W/m 2 Heat flux=19800W/m2
Surface temperature (K)
330
0.10
0.15
0.20
0.25
Fig. 11. Variation of maximum fuel temperature with axial distance for the axisymmetric model of an assembly of vertical cylinders.
portion surface temperature rises rapidly in this portion and gradually in the rest of the domain as heat transfer coefficient drops there. Maximum fuel temperature is the maximum temperature occurring at some axis in the fuel. As heat generation is uniform and flow conditions are symmetric so it might occur at centre line of the fuel. It has been calculated employing conduction considerations. Similar trends as attained at surface of fuel pin can be seen in Fig. 11 for maximum fuel temperature and can be explained on same basis.
360 350
0.05
Axial distance from the lower end (m)
Fig. 9. Variation of local Nusselt number with modified Rayleigh number Raz , for the axisymmetric model, of an assembly of vertical cylinders.
355
Heat flux=3600W/m 2 Heat flux=7200W/m 2 Heat flux=10800W/m 2 Heat flux=14400W/m 2 Heat flux=18000W/m2 Heat flux=19800W/m 2
325 320 315 310 305
4. Code validation
300 295
A comparative study of computer code COSINAC has been carried out against experimental results obtained in previous studies. A comparison of the average outlet coolant temperature measured experimentally and the average outlet coolant temperature calculated by code has been given in Table 2 and Fig. 12. It can be seen that code results fairly match with experimental observations. A small difference may be due to incorporation of Boussinesq approximation and assumption of constant material properties independent of temperature. Moreover, outlet temperature has been measured using a k-type thermocouple installed
290 0.00
0.05
0.10
0.15
0.20
0.25
Axial distance from the lower end (m) Fig. 10. Variation of surface temperature with axial distance, for the axisymmetric model of an assembly of vertical cylinders.
assembly, for different surface heat flux values studied. A rise in the surface temperature of cylinder can be observed along the length of cylinder as a uniform heat flux has been applied. Due to high local heat transfer coefficients in lower
Table 2 Steady state operating parameters of PARR-2 S.No.
1 2 3 4 5
Thermal neutron flux (#/cm2 s)
Reactor power
Average value
Nominal power (kW)
2.00 · 1011 ± 4.01 · 108 4.00 · 1011 ± 1.03 · 109 6.01 · 1011 ± 2.68 · 109 8.00 · 1011 ± 2.00 · 1011 1.00 · 1012 ± 2.00 · 1011
5.4 10.8 16.2 21.6 27.0
Inlet temperature Tin (C)
Experimental outlet temperature Tout (C)
COSINAC code result outlet temperature Tout (C)
% Nominal full power
Average value
Average value
Average value
20 40 60 80 100
17.6 18.5 19.5 21.3 23.2
23.6 29.4 33.9 38.2 41.8
23.73 28.10 32.05 36.24 40.22
Coolant Average Outlet Temperature ( o C)
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
201
end and approaches almost to unity at the outlet. During the steady state operation of the reactor, the maximum temperature of the fuel pin has been found to increase along the length but remains far less than the melting point. The local heat transfer coefficient has been found to decrease along the length of the fuel pin.
50 Code Results Experimental Values 40
30
References 20
10
0 20 %
40 %
60 %
80 %
100 %
Reactor power (% Nominal Full Power) Fig. 12. Comparison of average coolant outlet temperature obtained from code with experimental results.
near the outlet of the reactor, so the experimental value might not be the exact representation of the average outlet coolant temperature. 5. Conclusions A computer code COSINAC has been developed to study the steady state heat transfer behavior of PARR-2. Temperature and velocity profiles have been studied in detail. It has been found that at constant dimensionless distance (z/L) with increasing radial distance from the fuel rod the temperature decreases. At a radial distance of r = 3 mm while keeping the value of dimensionless distance equal to unity the temperature at the surface of the fuel rod is found to be 327 K. With increase in the radial distance temperature drop is very small. At the radial distance of 6.17 mm temperature drops to 293.6 K. Averages of these outlet temperatures at different wall heat fluxes computed by using the computer code COSINAC are very close to experimental values as shown in Table 2. The flow entrainment remains almost constant along the axial distance from the lower end of cylinder in assembly. The normalized Nusselt number (Nuz,/Nuz,Cylinder) decreases along the axial distance from the lower end being very large near the lower
[1] H.R. Nagendra, M.A. Tirunarayan, A. Ramachandram, Free convection from vertical cylinder with uniform heat flux, Journal of Heat Transfer, Transactions of the ASME Series C 92 (1970) 191– 194. [2] J.J. Heckel, T.S. Chen, B.F. Armaly, Natural convection along slender vertical cylinders with variable surface heat flux, Journal of Heat Transfer, Transactions of the ASME Series C 111 (1989) 1108– 1111. [3] K. Valusami, V.K. Garg, Transient natural convection over a heat generating vertical cylinder, International Journal of Heat and Mass transfer 35 (5) (1992) 1293–1306. [4] M. Keyhani, F.A. Kulacki, R.N. Christensen, Free convection in a vertical annulus with constant heat flux on the inner wall, Journal of Heat Transfer, Transactions of the ASME Series C 105 (1983) 454– 459. [5] I.R. Chughtai, M.H. Inayat, Numerical simulation of natural convection heat transfer in an open ended vertical annulus, Paper presented at the 6th National Symposium on Frontiers of Physics, Quaid-I-Azam University, Islamabad, Pakistan, 1997. [6] M. Keyhani, F.A. Kulacki, R.N. Christensen, Experimental investigation of free convection in a vertical rod bundle – a general correlation for Nusselt number, Journal of Heat Transfer, Transactions of the ASME Series C 107 (1985) 611–623. [7] Thermal hydraulic parameter calculations of the miniature neutron source reactor, MNSR-Document, China Institute of Atomic Energy, Beijing, 1988. [8] L.A. Khan, S.M. Mirza, N. Ahmed, M.S. Zafar, Study of self-power limiting characteristics of a tank-in-pool type research reactor for positive reactivity-induced over-power transients, Nuclear Science Journal 32 (4) (1995) 288–297. [9] L.A. Khan, N. Ahmed, M.S. Zafar, A study of the natural convection cooling of PARR-2 core, Paper presented at the 2nd All Pakistan Science Conference, Aitchison College, Lahore, Pakistan, 1993. [10] A.J. Friedland, C.F. Bonilla, Analytical study of heat transfer rates for parallel flow of liquid metals through tube bundles, American Institute of Chemical Engineers Journal 7 (1961) 107–122. [11] L.P. Davis, J.J. Perona, Development of free convection axial flow through a tube bundle, International Journal of Heat and Mass transfer 15 (1973) 1425–1438. [12] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980.
Computer simulation of natural convection heat transfer from an assembly of vertical cylinders of PARR-2 M. Abdul Basit a, Muhammad Rafique a
b,*
, Imran R. Chughtai a, Mansoor H. Inayat
a
Department of Chemical and Material Engineering, Pakistan Institute of Engineering and Applied Science Nilore, Islamabad 45650, Pakistan b Department of Physics, University of Azad Jammu and Kashmir Muzaffarabad 13100, Azad Kashmir, Pakistan Received 11 August 2005; accepted 1 May 2006 Available online 2 August 2006
Abstract In this paper a Computer Code COSINAC (Computer Simulation of Natural Convection from Assembly of vertical Cylinders) has been developed to simulate the natural convection heat transfer from an assembly of vertical cylinders of Pakistan Research Reactor-2 (PARR-2), under the steady state reactor operation. The momentum and energy equations in cylindrical co-ordinates, representing the thermal hydraulic behavior of a typical fuel pin in Pakistan Research Reactor-2, have been solved numerically for a two dimensional axisymmetric domain. The temperature and velocity profiles and Nusselt number variations have been studied and results have been presented. The computer code COSINAC has been validated against experimental results carried out in previous studies at different occasions. Average outlet coolant temperature simulated by computer code, at different wall heat fluxes, has been found in good agreement with experimental results. 2006 Elsevier Ltd. All rights reserved. Keywords: Natural convection; Pakistan research reactor-2; Assembly of vertical cylinders; Nusselt number
1. Introduction The nuclear reactors are classified on various concepts and basis. These may be classified on the basis of coolant or moderator used, type, the enrichment employed and arrangement of fuel. Most common classification is on the basis of application for which reactor is being used i.e. whether the reactor is meant for fuel production, power generation or research purposes. Whatever may be the type of reactor, an adequate cooling mechanism is a major requirement for safe operation of the reactor. Depending upon the quantity of heat generated by a reactor, different coolants and cooling mechanisms may be utilized. Power reactors are cooled by forced convection while research reactors, with small amount of heat generation, may
*
Corresponding author. Tel.: +92 03009189301/+92 05881034691; fax: +92 51 9223727. E-mail address: rafi[email protected] (M. Rafique). 1359-4311/$ - see front matter 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.applthermaleng.2006.05.001
employ free convection as the cooling mechanism. PARR-2, the reactor considered in the present study, consists of assembly of vertical cylinders cooled by natural convection. Extensive study on the natural convection cooling of thin vertical cylinders in various configurations has been done in the past. Steady state natural convection from a single cylinder with known heat flux on the wall was studied experimentally by Nagendra et al. [1] and numerically by Heckel et al. [2]. Valusami and Garge [3] have considered transient heat transfer from a vertical cylinder for a step change in surface heat flux. Keyhani et al. [4] studied steady state natural convection in a vertical annulus enclosed from both ends while Chughtai and Inayat [5] studied the same phenomena for a vertical annulus open from both ends. Keyhani et al. [6] have reported steady state natural convection from an assembly of vertical cylinders enclosed in a canister. For the reactor under consideration, simplified lumpedparameter analyses have been performed to model its
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
195
Nomenclature cp g k L Nuz
specific heat of fluid at constant pressure acceleration due to gravity thermal conductivity of fluid length of the assembly local Nusselt number for a typical fuel pin in assembly of vertical cylinders NuL local Nusselt number for a typical fuel pin in assembly of vertical cylinders at exit NuzCylinder local Nusselt number for vertical cylinder, of same dimensions as fuel pin, in an infinite medium (calculated within the code) q00 wall heat flux r radial position c qgbz4 q00 Ra* modified Rayleigh number = p k2 m2
thermal hydraulics [7–9]. Khan et al. [9] studied the natural convection cooling of PARR-2 in 1993. They used mathematical models to calculate coolant velocity as function of reactor power. The coolant velocity was also determined experimentally by operating the reactor at different power levels. A comparison between the calculated values and experimental values was also done. In the present study, an attempt has been made to study the thermal hydraulic phenomenon around a typical fuel pin of PARR-2. To achieve this end goal a single channel model has been used. An annular region filled with coolant, around a typical fuel pin in the reactor, has been assumed. Heat generation takes place in the fuel rod and heat is then transferred to coolant flowing parallel to the axis of the fuel pin in the annular region. The coolant enters the annular region from below and leaves from top. The thermal boundary-layer equations have been solved for the fluid flowing through the annular region around a typical fuel pin in the reactor, in cylindrical co-ordinates, using finite volume method. Density of fluid has been allowed to vary according to Boussinesq approximation and all other fluid properties have been assumed constant. Calculations have been made for velocity and temperature profiles and Nusselt numbers. 1.1. Reactor description PARR-2 is a tank-in-pool type reactor with a nominal power rating of 27 kW. It is composed of highly enriched uranium (HEU) fuel and is moderated and cooled by light water through natural convection. The HEU fuel is fabricated as thin vertical cylindrical pins of uranium–aluminum alloy with aluminum cladding. Each fuel pin has a heat generating length of 230 mm and an inert portion of 10 mm on both ends. Overall diameter of each fuel pin, including cladding thickness, is 5.9 mm. The fuel pins are arranged in concentric arrays with constant pitch of 10.95 mm so as to make a cylindrical shaped reactor core.
rm
r0 T T1 u v z a b m q
radius of the circle circumscribed about a fuel rod representing the flow area and fluid assignable to that particular fuel rod radius of fuel pin including cladding temperature bulk temperature of the fluid component of velocity in axial direction component of velocity in radial direction axial position thermal diffusivity of fluid volumetric coefficient of thermal expansion kinematics viscosity of fluid density of fluid
The reactor core is surrounded by radial as well as axial beryllium reflectors. There is a central control rod in the reactor core that is used to control the neutron flux and hence the reactor power. Energy generation rate, in the reactor, may be increased by withdrawing the control rod. However the control rod can be withdrawn only at the fixed rate and hence the reactivity insertion through control rod can only be a ramp function of time. The coolant flow path is restricted by placing two orifice plates on the top and bottom of the core (Fig. 1). The restricted coolant flow adds the self power-limiting feature to this reactor. The coolant enters the reactor core through the gap between the lower orifice plate and the annular reflector, rises up due to buoyancy and leaves through the upper orifice plate while cooling the fuel pins. Two k-type thermocouples have been installed near the lower and upper orifice plates to measure the inlet and outlet temperature of the coolant. A calibrated fission chamber type neutron detector has also been installed to measure the neutron flux and hence the reactor power. Part of the coolant coming out of the core is re-entrained into the in-going stream of the coolant and rest rises up into the tank due to buoyancy. The heat rejected from the core is first transferred to the water in the tank and then to the outer pool by natural convection. A brief summary of the reactor specification is given in Table 1. 2. Mathematical model 2.1. Physical and mathematical modeling The actual problem of heat transfer from PARR-2 has a complex 3-dimensional nature; however, due to axisymmetric geometry and physical phenomena, the problem can be studied efficiently in a two dimensional axisymmetric cylindrical coordinate system. It is the case of natural convection heat transfer, in axisymmetric geometry; hence flow is taking place due to buoyancy force caused
196
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
Fig. 1. Schematic diagram and coolant flow patterns in PARR-2 (not to scale).
Table 1 Summary of design features of PARR-2 Item
Value
Rated steady state power Enrichment in U-235 Maximum excess reactivity Thermal neutron flux at 27 kW Moderator and coolant Reflector
27 kW 90% 4.0 mk 1.0 · 1012 #/cm2 s Light water Beryllium
Fuel: Geometry Meal material Clad material Number
Cylindrical UAl4 – Al 303-1 Al 344
Control rod: Material Number
Cadmium 1
Detectors: Neutron Temperature
Fission chamber k-Type thermocouples
by constant heat flux applied on a fuel pin surface. As surface heat flux is same on all the fuel pins and a triangular pitch arrangement has been considered, so flow conditions may be taken same on all fuel rods and flow will be symmetric. The core of PARR-2 cannot be regarded as triangular or square pitch arrangement as shown in Fig. 2, but considering the stagger in fuel pin arrangement, a triangular pitch arrangement seems to be more appropriate. Following a procedure used by Friedland and Bonilla [10] for a tube bundle with an equilateral triangular pitch, a hexagon can be circumscribed about a tube to represent the flow area, and fluid assignable to that particular tube. It is assumed that this situation can be approximated closely enough by replacing the hexagon with a circle of equal area [11] having a radius rm as illustrated in Fig. 3. Coolant
flows through the annulus between fuel pin and this circular region. The dimensions of this annulus have been calculated by dividing the total cross sectional area on 354 fuel pins in PARR 2. Radius of inner cylindrical solid region is 2.95 mm (Radius of fuel pin). Around this solid region there is 3.22 mm thick annular region for coolant flow. So the total radial dimension is 6.17 mm. Consider an assembly of vertical cylinders (fuel rods) aligned in a quiescent (dormant) ambient fluid at temperature T1. Each cylinder in this assembly has a radius of r0. The z-coordinate is taken along axial direction from the leading edge of cylinder and r coordinate in radial direction from the axis of cylinder. The surface of the cylinder is subjected to a constant heat flux q00 . The acceleration due to gravity g is acting downward. All fluid properties are assumed to be constant except for density variation which induces the buoyancy force. Assumption of incompressible laminar boundary-layer flow and the Boussinesq approximation have been employed. The fluid velocities are quite small; therefore the compressive work and viscous dissipation terms in the energy balance equation can be neglected. The three governing conservation equations in the cylindrical co-ordinates, for natural convection, can be written as Continuity equation o o ðruÞ þ ðrvÞ ¼ 0 oz or
ð1Þ
z-Momentum equation u
ou ou t o ou þv ¼ r þ gbðT T 1 Þ oz or r or or
Energy balance equation oT oT a o oT þv ¼ r u oz or r or or
ð2Þ
ð3Þ
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
197
Fig. 2. Arrangement of fuel pins in the PARR-2 core.
Outflow Boundary Condition
Heated Cylinders Hexagonal Area Around Cylinder 6.17
2.95 11.83
φ5.50 Equivalent Unit Cell with Circular X-section Symmetry Boundary
R6.17
Wall Boundary (Constant Heat Flux) 230 Axis of cylindrical domain
z g
Hexagonal flow domain around a typical fuel pin r
Constant Pressure Boundary (T=293 K) 2-D Axisymmetric model used for modeling. (Not to Scale)
Fig. 3. Schematic diagram of the axisymmetric model for a triangular pitch assembly of vertical cylinders (all dimensions are in millimeters).
The boundary conditions for the solution of Eqs. (1)–(3) can be expressed, at the edge of the boundary-layer, as follows for constant heat flux:
For z = 0 and r P r0T = T1 00 For r = r0 and z > 0 u = 0, v = 0, oT ¼ qk or w For r = rm and z P 0 ou ¼ 0, v = 0, oT ¼0 or or
198
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
At z = L, outflow boundary condition [12] has been employed.
330 325
2.2. Method of solution
40x40 Grid points 20x20 Grid points 10x10 Grid points
320
overall trend observed in the domain. It can be seen that the three mesh sizes have produced almost identical results.
3. Discussion
3.2. Temperature and flow patterns
3.1. Grid independence
Fig. 5 shows velocity and temperature profiles for surface heat flux of 19,800 W/m2. In the case of assembly of vertical cylinders, flow is entrained only in axial direction. On one side, axial velocity is zero at the cylinder wall; on the other hand the outer side of the axisymmetric domain is the middle plane between two consecutive heated cylinders. Hence, resistance to flow of coolant is posed only at one side.
Three mesh sizes of 10 · 10, 20 · 20 and 40 · 40 grid points have been selected to study the grid independence of the code results. A uniform mesh size has been considered. Fig. 4 shows temperature profiles in the coolant for surface heat flux of 18,900 W/m2 at an axial distance of 0.158 m from the lower end. These profiles show a typical
Temperature (K)
Eqs. (1)–(3) constitute a system of coupled, non-linear, partial differential equations in the (r, z) co-ordinates. In the present study, a typical semi-implicit finite volume solution was carried out. The details of this numerical solution method are omitted to save space. Some of its highlights are being given. The first step in solving Eqs. (1)–(3) is to convert partial differential equations (2) and (3) given in non-conservation form into conservation form of partial differential equations. These conservation form of equations are then cast into algebraic equations using finite volume method employing Power Law scheme discussed by Patankar [12]. Axial velocity u has been calculated from discretized z-momentum equation (2) employing pointto-point Gauss–Siedel method. Radial velocity v has been calculated from discretized continuity equation. Finally, temperature distribution has been obtained by solving discretized energy equation using point-to-point Gauss–Siedel method.
290 0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
Radial Distance (mm) Fig. 4. Comparison of temperature profiles for natural convection from the axisymmetric model of an assembly of vertical cylinders using three mesh sizes (heat flux on hot surface = 19,800 W/m2, distance from the lower end = 0.158 m).
320
z/L=1.0
310
z/L=1.0
0.01
300
0.00 0.03
320
0.02 0.01 0.00 0.03 0.02
z/L=0.6
0.01 0.00 0.03 0.02
z/L=0.4
0.01
z/L=0.8
310
z/L=0.8
Temperature (K)
Axial velocity (m/s)
305
295
0.02
0.02
310
300
0.03
0.00 0.03
315
300 320
z/L=0.6
310 300 320 z/L=0.4
310 300
Surface of Cylinder
320 z/L=0.2
0.01
z/L=0.2
310 300
0.00 0.003
0.004
0.005
Radial distance (m)
0.006
0.0025 0.0030 0.0035 0.0040 0.0045 0.0050 0.0055 0.0060 0.0065
Radial distance (m)
Fig. 5. Profiles of axial velocity and temperature at various axial positions for natural convection in an axisymmetric model of an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
4
3
Nu z /Nu z Cylinder
Although, along the axial direction, shape of velocity profile changes according to the local buoyancy effects, yet area under a typical velocity profile almost remains unchanged. As the wall temperature is low in the lower portion, velocity profile has maxima at middle plane between two heated cylinders. But as the fluid moves axially along the cylinder, maxima of the velocity profile shifts inwards gradually due to increased buoyancy caused by increase in the surface temperature.
199
2
3.3. Flow entrainment In the case of assembly of vertical cylinders, fluid enters the domain from lower side only hence the entrainment remains almost unchanged with axial distance along the cylinder. Small slope of the line in Fig. 6 is due to temperature dependence of coefficient of volumetric expansion of the coolant.
1
0.00
0.05
0.10
0.15
0.20
Axial distance from the lower end (m) Fig. 7. The ratio of local Nusselt numbers Nuz/NuzCylinder for an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
3.4. Nusselt numbers 1.20 1.18 1.16 1.14
Nu L /Nu L,Cylinder
Fig. 7 shows a plot of Nuz/NuzCylinder versus axial distance along the surface of fuel rod. A vertical cylinder in infinite medium can entrain flow axially as well as radially, whereas a vertical cylinder in an assembly entrains fluid from lower end only, so flow rate in the case of assembly of vertical cylinders is higher in the lower portion as compared with cylinder in infinite medium. The phenomenon of fluid entrainment in assembly of vertical cylinders can be understood by comparing the situation to what happens in a chimney, where fluid rises as it becomes hot and surrounding cooler fluid rushes to the place to occupy the space. This is so called chimney effect. From Fig. 7, it can be seen that Nuz/NuzCylinder is well above 1.0 in lower region but in upper portion Nuz almost approaches to NuzCylinder.
1.12 1.10 1.08 1.06 1.04 1.02 1.00 1011
1012
1013
Modified rayleigh No. Ra* L Fig. 8. The ratio of overall Nusselt numbers Nuz/NuzCylinder for an assembly of vertical cylinders.
6.0x10-6
Flow entrainment (m3 /s)
5.0x10-6
The ratio of overall Nusselt number attained by a cylinder in an assembly to that of a cylinder in an infinite medium has been plotted in Fig. 8. It can be seen that due to so-called chimney effect, NuL for assembly is higher than that for an isolated cylinder for the whole range of RaL . Fig. 9 shows Nuz vs. Raz relation for a cylinder in assembly. It can be observed that Nuz vs. Raz data falls on a single curve for all values of surface heat flux studied in the present research. It is a well known characteristic of natural convection heat transfer that local Nusselt number is a function of modified Rayleigh number only.
4.0x10-6
3.0x10-6
2.0x10-6
1.0x10-6 0.00
0.05
0.10
0.15
0.20
0.25
Axial distance from the lower end (m) Fig. 6. Volumetric flow entrainment as a function of axial distance for an axisymmetric model of an assembly of vertical cylinders (heat flux on hot surface = 19,800 W/m2).
3.5. Surface and maximum fuel temperatures Fig. 10 shows temperature variation along fuel pin surface versus axial distance along the surface of cylinder in
200
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
Fuel maximum temperature (K)
10 2
Nu z
Heat flux=3600W/m2 Heat flux=7200W/m 2
10 1
Heat flux=10800W/m2 Heat flux=14400W/m2 Heat flux=18000W/m2 Heat flux=19800W/m2
10 0 103
104
105
106
107
108
109
1010
1011
1012
1013
380 375 370 365 360 355 350 345 340 335 330 325 320 315 310 305 300 295 0.00
Ra*z
345
Heat flux=3600W/m2 Heat flux=7200W/m 2 Heat flux=10800W/m2
340
Heat flux=14400W/m 2
335
Heat flux=18000W/m 2 Heat flux=19800W/m2
Surface temperature (K)
330
0.10
0.15
0.20
0.25
Fig. 11. Variation of maximum fuel temperature with axial distance for the axisymmetric model of an assembly of vertical cylinders.
portion surface temperature rises rapidly in this portion and gradually in the rest of the domain as heat transfer coefficient drops there. Maximum fuel temperature is the maximum temperature occurring at some axis in the fuel. As heat generation is uniform and flow conditions are symmetric so it might occur at centre line of the fuel. It has been calculated employing conduction considerations. Similar trends as attained at surface of fuel pin can be seen in Fig. 11 for maximum fuel temperature and can be explained on same basis.
360 350
0.05
Axial distance from the lower end (m)
Fig. 9. Variation of local Nusselt number with modified Rayleigh number Raz , for the axisymmetric model, of an assembly of vertical cylinders.
355
Heat flux=3600W/m 2 Heat flux=7200W/m 2 Heat flux=10800W/m 2 Heat flux=14400W/m 2 Heat flux=18000W/m2 Heat flux=19800W/m 2
325 320 315 310 305
4. Code validation
300 295
A comparative study of computer code COSINAC has been carried out against experimental results obtained in previous studies. A comparison of the average outlet coolant temperature measured experimentally and the average outlet coolant temperature calculated by code has been given in Table 2 and Fig. 12. It can be seen that code results fairly match with experimental observations. A small difference may be due to incorporation of Boussinesq approximation and assumption of constant material properties independent of temperature. Moreover, outlet temperature has been measured using a k-type thermocouple installed
290 0.00
0.05
0.10
0.15
0.20
0.25
Axial distance from the lower end (m) Fig. 10. Variation of surface temperature with axial distance, for the axisymmetric model of an assembly of vertical cylinders.
assembly, for different surface heat flux values studied. A rise in the surface temperature of cylinder can be observed along the length of cylinder as a uniform heat flux has been applied. Due to high local heat transfer coefficients in lower
Table 2 Steady state operating parameters of PARR-2 S.No.
1 2 3 4 5
Thermal neutron flux (#/cm2 s)
Reactor power
Average value
Nominal power (kW)
2.00 · 1011 ± 4.01 · 108 4.00 · 1011 ± 1.03 · 109 6.01 · 1011 ± 2.68 · 109 8.00 · 1011 ± 2.00 · 1011 1.00 · 1012 ± 2.00 · 1011
5.4 10.8 16.2 21.6 27.0
Inlet temperature Tin (C)
Experimental outlet temperature Tout (C)
COSINAC code result outlet temperature Tout (C)
% Nominal full power
Average value
Average value
Average value
20 40 60 80 100
17.6 18.5 19.5 21.3 23.2
23.6 29.4 33.9 38.2 41.8
23.73 28.10 32.05 36.24 40.22
Coolant Average Outlet Temperature ( o C)
M.A. Basit et al. / Applied Thermal Engineering 27 (2007) 194–201
201
end and approaches almost to unity at the outlet. During the steady state operation of the reactor, the maximum temperature of the fuel pin has been found to increase along the length but remains far less than the melting point. The local heat transfer coefficient has been found to decrease along the length of the fuel pin.
50 Code Results Experimental Values 40
30
References 20
10
0 20 %
40 %
60 %
80 %
100 %
Reactor power (% Nominal Full Power) Fig. 12. Comparison of average coolant outlet temperature obtained from code with experimental results.
near the outlet of the reactor, so the experimental value might not be the exact representation of the average outlet coolant temperature. 5. Conclusions A computer code COSINAC has been developed to study the steady state heat transfer behavior of PARR-2. Temperature and velocity profiles have been studied in detail. It has been found that at constant dimensionless distance (z/L) with increasing radial distance from the fuel rod the temperature decreases. At a radial distance of r = 3 mm while keeping the value of dimensionless distance equal to unity the temperature at the surface of the fuel rod is found to be 327 K. With increase in the radial distance temperature drop is very small. At the radial distance of 6.17 mm temperature drops to 293.6 K. Averages of these outlet temperatures at different wall heat fluxes computed by using the computer code COSINAC are very close to experimental values as shown in Table 2. The flow entrainment remains almost constant along the axial distance from the lower end of cylinder in assembly. The normalized Nusselt number (Nuz,/Nuz,Cylinder) decreases along the axial distance from the lower end being very large near the lower
[1] H.R. Nagendra, M.A. Tirunarayan, A. Ramachandram, Free convection from vertical cylinder with uniform heat flux, Journal of Heat Transfer, Transactions of the ASME Series C 92 (1970) 191– 194. [2] J.J. Heckel, T.S. Chen, B.F. Armaly, Natural convection along slender vertical cylinders with variable surface heat flux, Journal of Heat Transfer, Transactions of the ASME Series C 111 (1989) 1108– 1111. [3] K. Valusami, V.K. Garg, Transient natural convection over a heat generating vertical cylinder, International Journal of Heat and Mass transfer 35 (5) (1992) 1293–1306. [4] M. Keyhani, F.A. Kulacki, R.N. Christensen, Free convection in a vertical annulus with constant heat flux on the inner wall, Journal of Heat Transfer, Transactions of the ASME Series C 105 (1983) 454– 459. [5] I.R. Chughtai, M.H. Inayat, Numerical simulation of natural convection heat transfer in an open ended vertical annulus, Paper presented at the 6th National Symposium on Frontiers of Physics, Quaid-I-Azam University, Islamabad, Pakistan, 1997. [6] M. Keyhani, F.A. Kulacki, R.N. Christensen, Experimental investigation of free convection in a vertical rod bundle – a general correlation for Nusselt number, Journal of Heat Transfer, Transactions of the ASME Series C 107 (1985) 611–623. [7] Thermal hydraulic parameter calculations of the miniature neutron source reactor, MNSR-Document, China Institute of Atomic Energy, Beijing, 1988. [8] L.A. Khan, S.M. Mirza, N. Ahmed, M.S. Zafar, Study of self-power limiting characteristics of a tank-in-pool type research reactor for positive reactivity-induced over-power transients, Nuclear Science Journal 32 (4) (1995) 288–297. [9] L.A. Khan, N. Ahmed, M.S. Zafar, A study of the natural convection cooling of PARR-2 core, Paper presented at the 2nd All Pakistan Science Conference, Aitchison College, Lahore, Pakistan, 1993. [10] A.J. Friedland, C.F. Bonilla, Analytical study of heat transfer rates for parallel flow of liquid metals through tube bundles, American Institute of Chemical Engineers Journal 7 (1961) 107–122. [11] L.P. Davis, J.J. Perona, Development of free convection axial flow through a tube bundle, International Journal of Heat and Mass transfer 15 (1973) 1425–1438. [12] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, New York, 1980.