Applied Thermal Engineering 36 (2012) 288e295

Contents lists available at SciVerse ScienceDirect

Applied Thermal Engineering

journal homepage: [Link]/locate/apthermeng

CFD analysis of a CO2 based natural circulation loop with end heat exchangers
Ajay Kumar Yadav, M. Ram Gopal*, Souvik Bhattacharyya
Department of Mechanical Engineering, Indian Institute of Technology Kharagpur, Kharagpur 721302, India

a r t i c l e i n f o a b s t r a c t

Article history: 3-D steady flow simulation of a CO2 based, single-phase (subcritical and supercritical) rectangular natural
Received 4 August 2011 circulation loop (NCL) with end heat exchangers is presented. Viscous dissipation in fluid and axial
Accepted 18 October 2011 conduction in fluid as well as in solid wall are considered. Results are obtained for various inlet
Available online 3 November 2011
temperatures of water in the hot heat exchanger, for a fixed inlet temperature of cooling water in the cold
heat exchanger. Effect of loop operating pressure (60e100 bar) on system performance is also investi-
Keywords:
gated. Results show that due to the presence of bends and strong buoyancy effects near pseudo-critical
Natural circulation loop
zone, local velocity and temperature vary in all three dimensions. It is noticed that at a given loop
CFD
Turbulent flow
pressure, the variation of heat transfer rate with loop fluid temperature changes sign in the vicinity of
Supercritical carbon dioxide pseudo-critical region. Steady-state Reynolds number and heat transfer coefficients obtained through the
Heat transfer CFD analysis are compared with those calculated using available correlations. Although qualitative trends
match, quantitative differences exist that may be attributed to difference in configuration. Additionally,
new correlations are proposed for Re in terms of modified Gr, and Nu in terms of Re and Pr which are
expected to be useful in the design and analysis of NCLs with end heat exchangers.
Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Zhang et al. [11] and Chen et al. [12] reported studies on the effects
of heat transfer and the instabilities of supercritical CO2 flow in
Natural circulation loop (NCL) based secondary fluid systems a 2-D NCL at a fixed operating pressure of 90 bar operating over
are simple and reliable due to the absence of any moving a large heat source temperature range. However, their analysis
components such as pumps. While water based NCLs are widely considered the case of isothermal heat source and heat sink only.
used in applications such as solar collectors, nuclear reactors, etc, There are many practical applications in which end heat
a wide variety of brines is used for low temperature applications. exchangers are used for heat source and sink (e.g. refrigeration and
In recent years, a growing popularity of carbon dioxide as the air-conditioning systems, solar collectors, solar power plants). In
secondary fluid has been witnessed in both forced as well as in studies on estimation of in-tube cooling heat transfer coefficient
natural circulation loops. This may be attributed to the favourable for turbulent forced convection of supercritical CO2, correlations
thermo-physical properties of CO2 in addition to its environment developed by Pitla et al. [13] and Du et al. [14] are usually
friendliness. Studies show that for low temperature refrigeration employed. However, Du et al. [14] suggested the use of 3-D model
and air-conditioning applications, CO2 based NCLs are very to study the behaviour of supercritical CO2 near the pseudo-critical
compact in comparison to other conventional working fluids [1] point. Studies on heat transfer coefficient for supercritical CO2
and have been proposed for various heat transfer application based NCLs are scarce in literature. In addition, to account for the
such as new generation nuclear reactors [2], in chemical extraction strong buoyancy effects (local) near pseudo-critical zone, and the
[3,4], cryogenic refrigeration [5], heat pump [6], electronic cooling effect of bends in pipe etc., it becomes essential to consider
systems [7], in geothermal applications [8,9], etc. However, a three-dimensional model for greater accuracy. However, inves-
detailed modelling and analyses of CO2 based NCLs are relatively tigations employing 3-D models of CO2 based NCLs with end heat
sparse in the literature. Kiran Kumar and Ram Gopal [10] reported exchangers are not available in the open literature. To fill in that
one-dimensional steady-state analysis of a rectangular NCL with void, this study presents CFD analysis of a 3-D CO2 based NCL with
end heat exchangers for low temperature applications. Recently end heat exchangers. Results are presented on the steady-state
behaviour of the loop at various operating pressures and temper-
atures. Operating parameter range is chosen such that the loop
* Corresponding author. Tel.: þ91 3222 282986. fluid (CO2) exists either as a subcritical single-phase fluid or as
E-mail address: ramg@[Link] (M. Ram Gopal). a supercritical single-phase fluid.

1359-4311/$ e see front matter Ó 2011 Elsevier Ltd. All rights reserved.
doi:10.1016/[Link].2011.10.031
 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295 289

Nomenclature Q heat transfer rate (W)

Re Reynolds number
A area (m2) T temperature of fluid (K)
cp specific heat capacity (J/kgK) V velocity (m/s)
d diameter of inner pipe or loop diameter (m)
D diameter of outer pipe of CHX and HHX (m) Greek letters
g gravitational acceleration (m/s2) b volumetric expansion coefficient (1/K)
Grm modified Grashof number ΔT temperature difference across the CHX/HHX (K)
h heat transfer coefficient (W/m2K) m dynamic viscosity (kg/ms)
H0 total height of vertical pipes (m) r density of fluid (kg/m3)
k thermal conductivity (W/mK)
L length of CHX and HHX (m) Subscripts
L0 total length of a horizontal pipe (m) avg average
L1 adiabatic pipe length on horizontal pipe (m) CHX cold heat exchanger
Lt total length of the loop (m) HHX hot heat exchanger
m mass flow rate (kg/s) m modified
Pr Prandtl number w water

2. Physical model and mathematical formulations The following geometric and material specifications are
considered for the model:
2.1. Physical model Internal diameter of NCL (internal pipe), d ¼ 4 cm
Internal diameter of CHX and HHX (external pipe), D ¼ 6 cm
Fig. 1 shows the schematic of a rectangular NCL which consists Length of CHX and HHX, L ¼ 120 cm
of a cold heat exchanger (CHX), a hot heat exchanger (HHX), a riser Total width and height of NCL (L0  H0) ¼ 154 cm  124.5 cm
and a downcomer. The loop fluid is heated sensibly by extracting Total length of the loop (Lt) ¼ 545 cm
heat from the external fluid in HHX and is cooled sensibly by Insulated pipe length in horizontal pipe (L1) ¼ 17 cm
rejecting heat to the external fluid in the CHX. Circulation of the Tube wall thickness ¼ 3 mm; material: copper.
loop fluid is maintained due to the buoyancy effect caused by
heating at the bottom and cooling at the top.
2.2. Mathematical formulation and solution
The following assumptions are made in the analysis:

The governing mass, momentum and energy equations are

(a) The internal (CO2) and external (water) fluid streams are in
solved using the CFD code, FLUENT. A three-dimensional geometry
single-phase throughout the loop.
was prepared in GAMBIT 2.4.6 and steady-state simulation was
(b) The system is operating at steady-state.
carried out employing Fluent 6.3.26, where the implicit-coupled
(c) The riser, downcomer and external pipes of CHX and HHX are
finite-volume method was used to discretize the governing equa-
perfectly insulated.
tions. The pressure-implicit with splitting of operators (PISO) algo-
(d) The external hot and cold fluid properties are estimated at inlet
rithm was used to solve the coupling model between velocity and
temperature as the temperature variation is kept small.
pressure.
(e) Fouling in both the heat exchangers are neglected.
The momentum and energy terms in the governing equations
are iterated with a second-order upwind scheme that uses the
upstream values and gradients to compute the control volume face
values. Turbulent parameters (k, ε etc.) are also iterated with
a second-order upwind scheme. The PRESTO (Pressure Staggering
Option) scheme is used to discretize the pressure term. No-slip
boundary condition is used on the pipe walls. Turbulent models
for supercritical fluid are less developed and still under intense
study [15]. Therefore, in the present simulation a general Renorm-
alization Group (RNG) k  ε model is selected as the first step to
introduce the expression of turbulent effect. This method has also
been used successfully in previous studies on supercritical CO2
turbulent flow yielding accurate results [16]. For the near wall
treatment, a standard wall function has been assumed in the case of
turbulent flow [17]. Axial conduction and viscous dissipation in
fluid are considered. Axial conduction along the tube wall is also
taken into account. Convergence is obtained when various resid-
uals of the parameters (temperature, velocity, pressure etc.) do not
change with the iterations. Convergence criterion is maintained
below103 for the residuals of continuity and below106 for energy.
Conservation of mass and energy are also checked for all the cases
in the analysis.
Fig. 2 shows the meshing of cross section of the loop which has
Fig. 1. Schematic of NCL with end heat exchangers (CHX and HHX). a maximum grid size of 4 mm in radial direction and coarse
290 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295

ZA
cp T rVdA
0
T ¼ (2)
ZA
cp rVdA
0

Steady-state Reynolds number and modified Grashof number

[18] are defined as,

4m
Re ¼ (3)
pdm

g bd3 r2 QH0
Grm ¼ (4)
Am3 cp
where, Q ¼ mw cp;w DTw
All properties for Re and Grm are calculated at the average
temperature (Tavg) of the loop, which is expressed as,

TCHX þ THHX
Tavg ¼ (5)
2
where, TCHX and THHX are the water inlet temperatures at CHX and
HHX, respectively.
Fig. 2. Meshing of a cross section.
Modified Grashof number expressed in Eq. (4) is somewhat
different from classical Grashof number. Bau and Torrance [19]
prefer to call it a non-dimensional heating rate. However, the
meshing is done for axial direction (5 mm grid size in CHX and term QH/(Amcp) has the same dimensions of temperature difference
HHX, 20 mm for riser and downcomer). Mesh generation, imple- and if we consider it as reference temperature difference, then it is
mented in the GAMBIT platform, yielded a total of 131,904 nodes. appropriate to refer to it as a modified Grashof number [18].
Grid independence test was carried out to ensure optimal choice of Recently most of the authors [11,20,21] are using same definition of
fineness. For computational purposes, the following terms are modified Grashof number in the study of natural circulation loop.
defined.
Mass flow rate at any cross section is defined as,
2.3. Calculation of thermo-physical properties of CO2
ZA
rVdA Since the state of CO2 inside the loop varies from subcritical to
m ¼ (1)
supercritical and the variation in thermo-physical properties near
0
the critical point is extremely large, it is essential to calculate the
Mass weighted average (bulk mean) temperature of the fluid at properties by taking into account the temperature variation prop-
any cross section is defined as, erly. However, as shown in the literature, due to very small variation
in operating pressure throughout the NCL the effect of variation of
pressure on the properties of single-phase carbon dioxide is not
expected to be significant [10,11]. For example, results from the
14
present study show that the maximum variation of CO2 pressure is
90bar
12 35
Specific heat, c p (kJ/kgK)

cp from NIST database

30
10 cp curve fit for present study
60 bar
Specific heat (kJ/kgK)

25 70 bar
8 80 bar
90 bar
20 100 bar

6
15

4 10

5
2
0
304 312 320 328 336 344 352 250 275 300 325 350 375 400
Temperature (K) Temperature (K)

Fig. 3. Curve fit data for specific heat of supercritical CO2. Fig. 4. Variation of specific heat of CO2 with temperature at various pressures.
 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295 291

10.22 kPa at an operating pressure of 90 bar (0.11% of operating with specific heat reaching the maximum at a given pressure. Fig. 4
pressure). Hence, for a given operating pressure, the properties of exhibits the effect of operating pressure on pseudo-critical point
CO2 at any point in the loop are calculated at the fixed operating with bulk mean temperature of fluid. As pressure increases,
pressure and local temperature. The required properties of CO2 pseudo-critical point for specific heat shifts towards higher
including density, specific heat, thermal conductivity and viscosity temperature.
are obtained from the NIST Standard Reference Database REFPROP
[22] and are considered to be temperature dependent in the calcu- 3. Results and discussions
lation. Properties of CO2 for the operating temperature range at
a temperature difference of 1 K are added to the fluid properties list In this study, water is employed as the external fluid in both CHX
in the Fluent, and a piecewise-linear interpolation approach is used and HHX. Inlet temperature of water in CHX is kept constant at
to calculate the properties within 1 K temperature difference. Fig. 3
shows the comparison of specific heat data used in Fluent for this
study with REFPROP data. Maximum variation with temperature is a 324
0.7
observed in case of specific heat (cp), while for the remaining
thermo-physical properties (density, thermal conductivity and 322
0.6
viscosity) the regression correlation coefficient is superior.
320 0.5

Temperature (K)

Velocity (m/s)
2.4. Pseudo-critical region
318 0.4

Pseudo-critical region is characterised by the temperature at

T, sec. 5 V, sec. 5 0.3
which thermo-physical properties of the fluid (specific heat, 316
T, sec. 6 V, sec. 6
density, viscosity and conductivity) exhibit abrupt fluctuations T, sec. 7 V, sec. 7 0.2
314
0.1
a 312
330 0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
CHX RIGHT LEG HHX LEFT LEG
Diameter (cm)
325
b 328
Temperature (K)

313 K; 323 K 0.7

333 K; 343 K
353 K 326 0.6
320

0.5
Temperature (K)

324 T, sec. 11 V, sec. 11

Velocity (m/s)
315 T, sec. 12 V, sec. 12
0.4
T, sec. 13 V, sec. 13
322
0.3
310
320 0.2

0 1 2 3 4 5 0.1
Length (m) 318
0.0
b 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
CHX RIGHT LEG HHX LEFT LEG
2 Diameter (cm)
T o tal P ressu re d ifferen ce (kP a)

c 0.75
0
320.8
0.60
-2
320.4
Temperature (K)

Velocity (m/s)

0.45
-4
313 K
323 K 320.0 T, sec. 1 0.30
-6 V, sec. 1
333 K T, sec. 2 V, sec. 2
343 K T, sec. 3 V, sec. 3
T, sec. 4 V, sec. 4 0.15
-8 353 K 319.6

-10 0.00
319.2
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0
0 1 2 3 4 5
Diameter (cm)
Length (m)
Fig. 6. (a) Temperature and velocity profile at different cross sections in CHX at 333 K.
Fig. 5. (a) Variation of bulk mean temperature throughout the loop length. (b) Vari- (b) Temperature and velocity profile at different cross sections in HHX at 333 K.(c)
ation of total pressure difference throughout the loop length. Temperature and velocity profile at different cross sections in riser at 333 K.
292 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295

Fig. 7. Temperature contour plot for HHX water inlet temperature of 333 K at (a) Left leg (section 3), (b) CHX (section 6), (c) Right leg (section 9) and (d) HHX (section 12).

305 K and that in HHX is varied from 313 K to 353 K in steps of 10 K. region where density of the fluid is maximum leading to higher
Results are obtained for five operating pressures of CO2, i.e., 60, 70, pressure drop.
80, 90 and 100 bar. Operating pressure of the system is defined at
the centre of HHX. Mass flow rate of water in CHX as well as in HHX 3.2. Temperature and velocity profile at different cross sections
is assumed to remain constant at a value of 0.3 kg/s which ensures
turbulent flow on water side. Fig. 6 (a), (b) and (c) show the local temperature and local axial
velocity profiles at different cross sections in the CHX, HHX, and left
leg (riser) at a hot external fluid inlet temperature of 333 K and
3.1. Variation of temperature and pressure throughout the loop at a CO2 pressure of 90 bar. Variation of temperature and velocity
90 bar shown in the figures are drawn across the diameter at different
cross sections of the loop as defined in Fig. 1. Due to large variation
Fig. 5(a) and (b) show the variation of bulk mean temperature in density with temperature at this pressure strong buoyancy is
and total pressure difference of CO2 through the loop for varying produced. Unlike vertical sections, the combined effect of gravity
water inlet temperatures of HHX at an operating pressure of 90 bar. and local buoyancy in horizontal sections (CHX and HHX) leads to
Flow reversal is observed at a hot fluid inlet temperature of 343 K asymmetric variation of velocity and temperature. Thus, tempera-
and thus the functions of riser and downcomer are reversed. Hence, ture of fluid in the horizontal segments (CHX and HHX) increases
the CHX and HHX become parallel flow heat exchangers, since the toward the top at all the segments except at the exit of bends (in the
flow direction of external fluid is kept the same for all conditions. diagram, 0 is the bottom and 4 cm is the top). In the riser, at
One of the possible reasons for flow reversal is that the loop is a section far away from the bend (section 4), there is negligible
thermally and geometrically symmetrical. In this study it is temperature variation across the diameter of pipe. Velocity profiles
observed that external fluid temperature difference in CHX or HHX are also influenced by the centrifugal action caused by bends. Inside
is less than 1 K which shows that the loop is thermally almost the close loop, fluid at the outer part has higher velocity than at the
symmetrical. Chen et al. [12] also observed a similar flow reversal inner part, as is exhibited by all the plots.
phenomenon. Due to the assumption of adiabatic condition in the
riser and downcomer, the variation in temperature is almost 3.3. Temperature and velocity contour plots
negligible in these parts, whereas variations can be seen in CHX and
HHX due to heat transfer with the external fluid. Variation of total Figs. 7 and 8 show the velocity (m/s) and temperature (K)
pressure can be seen in riser and downcomer due to height of the contour plot for HHX water inlet temperature of 333 K at the centre
loop. Maximum variation of pressure in a loop is 10.22 kPa (0.11%) of left leg (section 3), CHX (section 6), right leg (section 9) and HHX
at an operating pressure of 90 bar. As temperature increases, (section 12). The isotherms and velocity vectors clearly show the
pressure drop for the loop decreases which is due to the decrease in non-uniform temperature and velocity variation inside the tube at
density at higher temperature. At a pressure of 90 bar and HHX a given cross section. As mentioned before, temperature at upper
water inlet temperature of 313 K, loop fluid is in the pseudo-critical part of CHX and HHX is higher than that at the lower part due to

Fig. 8. Velocity contour plot for HHX water inlet temperature of 333 K at (a) Left leg (section 3), (b) CHX (section 6), (c) Right leg (section 9) and (d) HHX (section 12).
 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295 293

local buoyancy effect. In the riser (Fig. 7(a)) and downcomer comparison to mass flow rate starts at a lower pressure with HHX
(Fig. 7(c)), temperature is almost uniform over the cross section as water inlet temperature remaining the same. To know the pseudo-
there is no heat transfer from the wall except axial conduction. critical point, bulk mean temperature is needed which is taken as
Velocity contour also shows the same trend as explained in Fig. 6. an approximate average temperature of CHX and HHX water inlet
Since the temperature variation has a bearing on the stability of the temperatures.
loop, it is essential to capture it in order to understand the loop Fig. 9(b) also shows that as water inlet temperature of HHX
behaviour. This evidently justifies the use of a 3-dimensional increases, the heat transfer rate increases due to increase in
model. buoyancy effect. Even though there is decrease in mass flow rate at
higher temperature at a given pressure due to lower density, heat
3.4. Variation of loop fluid mass flow rate and loop heat transfer transfer rate does not decrease. This may be attributed to the fact
rate that at constant pressure, with increase in temperature, density of
the fluid decreases and near the critical point the viscosity and
Fig. 9(a) and (b) show the variation of CO2 mass flow rate and coefficient of thermal expansion decrease with temperature.
loop heat transfer rate with operating pressure at various HHX Decrease in viscosity with increase in temperature enhances the
water inlet temperatures. Results show that as operating temper- Reynolds number that contributes in enhancing heat transfer rate.
ature of loop fluid increases, mass flow rate as well as heat transfer
rate increase up to the pseudo-critical point and starts decreasing 3.5. Variation of heat transfer coefficient
beyond that. Mass flow rate and heat transfer rate are influenced
mainly by density and specific heat, respectively. Pseudo-critical Fig. 10 shows the variation of CO2 heat transfer coefficient with
point (temperature) for density at any pressure is lower than that operating pressure at different HHX water inlet temperatures for
of specific heat. In other words for a constant temperature situa- cooling. Heat transfer coefficient based on wall function is pre-
tion, pseudo-critical pressure for density is higher than that of sented here for cooling (CHX). Results show that at higher pressure
specific heat. Therefore, decrease in heat transfer rate in (90 and 100 bar), as operating temperature of loop fluid increases
heat transfer coefficient also increases initially and reaches
a maximum and then starts decreasing at higher temperature.
a 0.35
Results also show that as operating pressure increases peak of the
graph shifts towards the higher temperature due to shifting of the
313 K
0.30 pseudo-critical point. Heat transfer coefficient for heating is also
323 K
333 K calculated for all the cases and a similar trend is observed.
Mass flow rate (kg/s)

343 K
0.25 353 K 3.6. Comparison of heat transfer coefficients

0.20 Fig. 11 shows the comparison of obtained heat transfer coeffi-

cient employing the well known DittuseBoelter and Gnielinsky
correlations [23,24] at an operating pressure of 90 bar. Correlation
0.15 proposed by Filonenko [25] is used to calculate the friction factor
required for Gnielinsky correlation. Heat transfer coefficient based
0.10 on wall function is obtained from simulation for cooling (CHX) and
heating (HHX). Thermo-physical properties of the loop fluid are
estimated at the average temperature for the correlations. The
0.05 same Reynolds number as was calculated from mass flow rate at
60 70 80 90 100
bulk mean temperature is used in the correlations to obtain the
Pressure (bar)
2800
b 2400
Cooling 60 bar
70 bar
Heat transfer coefficient (W/m K)

313 K 2400 80 bar

2

2000 323 K
333 K 90 bar
100 bar
Heat transfer rate (W)

343 K 2000
1600 353 K

1600
1200

1200
800
800
400
400
0
60 70 80 90 100 310 320 330 340 350
Pressure (bar) Temperature (K)
Fig. 9. (a) Variation of mass flow rate of CO2 with HHX inlet temperature. (b) Variation Fig. 10. Variation of heat transfer coefficient for different operating pressure for
of heat transfer rate of CO2 with HHX inlet temperature. cooling (CHX).
294 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295

heat transfer coefficient. Even though the entire set of graphs

exhibits extremely close agreement qualitatively, quantitatively
there is some difference. The predictions from this study do exhibit
much better comparison with those obtained from Gnielinsky
correlation. It may be noted that these correlations employed here
are broad-spectrum general correlations for fully developed forced
circulation and may not be very suitable for natural circulation
supercritical CO2. It is also expected that near the pseudo-critical
region, the fluid flow cannot be fully developed due to large
property variation. The plot clearly shows that the heat transfer
coefficient obtained in the present study using CFD is higher than
that obtained through the correlations. In addition, it may be noted
that the obtained cooling heat transfer coefficient is slightly higher
than that in case of heating.

3.6.1. Proposed correlation for Nusselt number

Based on the generated data, new correlations have been Fig. 12. Validation of obtained result with experimental data.
developed for cooling and heating. The proposed correlation is
dimensional parameters Reynolds number (Re) and modified Gra-
based on the functional form of DittuseBoelter correlation with
shof number (Grm) which are calculated at the average temperature
modified values of constant and indices.
(Tavg) of the loop. Fig. 12 shows that even though there is a good
Nu ¼ CRem Prn (6) trend match between the results, quantitative difference exists
between the numerical and the experimental results. This may be
where, C, m and n are constants. The proposed correlation is attributed to the dissimilarity in the method of heating the loop
specified by the following: fluid in the two studies. Hence, to predict the performance accu-
In case of cooling: C ¼ 0.02279, m ¼ 0.8231, n ¼ 0.4216. rately, it is essential to analyze the NCL with loop heat exchangers
In case of heating: C ¼ 0.02834, m ¼ 0.8018, n ¼ 0.4454. separately instead of extrapolating the data obtained on either
The proposed correlation exhibits a R2 value of 99.91% for constant heat flux or isothermal heat source and heat sinks.
cooling and 99.28% for heating reflecting a good fit with the It may be noted that for a very small temperature difference of
simulated data. The proposed correlation has been validated for the 8 K between CHX and HHX, Re is of the order of 105 which also
subcritical and supercritical CO2, air, and water for the following contributes to enhanced heat transfer.
ranges of Re and Pr: Based on the generated data, a new correlation has been
developed for Re in terms of modified Gr for turbulent flow. The
7 000  Re  360 000; 0:78  Pr  6:1
proposed correlation is functionally similar to the Vijayan correla-
tion reported earlier [20] with differing coefficients and indices.
Correlation has been developed employing a software Engineering
3.7. Validation of obtained results with published experimental Equations Solver (EES) V8.325. The proposed correlation exhibits
data a R2 value of 94.2%.
Vijayan correlation [21]:
In order to validate the results obtained from the CFD analysis, 1
they are compared with the experimental data reported earlier by Re ¼ 1:96ðGrm d=Lt Þ2:75 (7)
Vijayan [20], which was obtained on a water based NCL with an
electrical heater (with constant wall heat flux) in place of a hot heat Proposed correlation :
exchanger. Hence, comparison is made in terms of the non- 1
Re ¼ 1:77ðGrm d=Lt Þ2:76 ð7 000  Re  360 000Þ (8)
2700
Above correlation would be useful for designing the natural
90 bar
Fluent (Cooling) circulation loop with end heat exchanger of any dimensions within
Heat transfer coefficient (W/m K)

2400 Fluent (Heating) the above mentioned range of Reynolds number. Correlation has
2

Dittus-Boelter (Cooling)
been developed for a steady, single-phase natural circulation loop.
Dittus-Boelter (Heating)
2100 Gnielinsky
All Properties used in the correlation are taken at average
temperature of the loop.
1800
4. Conclusions

1500
In the present study, results obtained from CFD analysis of
a carbon dioxide based, single-phase (subcritical and supercritical)
1200 NCL with end heat exchangers have been presented. Steady flow
simulation of the 3-dimensional model of a rectangular NCL has
900 been implemented. The following conclusions are made after the
discussion of the different cases:
600
310 320 330 340 350 (i) Results show that due to the presence of bends and strong
Temperature (K) buoyancy effects near pseudo-critical zone, fluid parameters
such as local velocity and temperature vary in all three
Fig. 11. Comparison of heat transfer coefficient with different correlations at 90 bar. dimensions.
 A.K. Yadav et al. / Applied Thermal Engineering 36 (2012) 288e295 295

(ii) Results show that as operating temperature of loop fluid [6] K. Ochsner, Carbon dioxide heat pipe in conjunction with a ground source
heat pump (GSHP), Appl. Therm. Eng. 28 (2008) 2077e2082.
increases mass flow rate as well as heat transfer rate increase
[7] D.E. Kim, M.H. Kim, J.E. Cha, S.O. Kim, Numerical investigation on thermal-
up to the pseudo-critical point and starts decreasing after this hydraulic performance of new printed circuit heat exchanger model, Nucl.
point. Eng. Des. 238 (2008) 3269e3276.
(iii) Heat transfer rate in the supercritical case is higher than that [8] D.B. Kreitlow, G.M. Reistad, Themosyphon models for downhole heat
exchanger application in shallow geothermal systems, J. Heat Transfer 100
for the subcritical scenario. To achieve the maximum heat (1978) 713e719.
transfer rate, operating condition should be chosen near the [9] K.E. Torrance, Openeloop themosyphons with geological application, J. Heat
pseudo-critical region in supercritical zone. Transfer 100 (1979) 677e683.
[10] K. Kiran Kumar, M. Ram Gopal, Steady-state analysis of CO2 based natural
(iv) As operating pressure increases peak of the maximum heat circulation loops with end heat exchangers, Appl. Therm. Eng. 29 (2009)
transfer coefficient shifts towards the higher operating 1893e1903.
temperature. [11] X. Zhang, L. Chen, H. Yamaguchi, Natural convective flow and heat transfer of
supercritical CO2 in a rectangular circulation loop, Int. J. Heat Mass Transf. 53
(v) Reynolds number for CO2 reaches a high value of the order of (2010) 4112e4122.
105 even for a small temperature difference of 8 K and oper- [12] L. Chen, X. Zhang, H. Yamaguchi, Z. (Simon) Liu, Effect of heat transfer on the
ating pressure of 60 bar which indicates efficient heat transfer. instabilities and transitions of supercritical CO2 flow in a natural circulation
loop, Int. J. Heat Mass Transf. 53 (2010) 4101e4111.
(vi) New correlations are proposed for Re in terms of modified Gr,
[13] S.S. Pitla, E.A. Groll, S. Ramadhyani, New correlation to predict the heat
and Nu in terms of Re and Pr. It is expected that these corre- transfer coefficient during in-tube cooling of turbulent supercritical CO2, Int. J.
lations will be useful in the design and analysis of NCLs with Refrig. 25 (2002) 887e895.
[14] Z. Du, W. Lin, A. Gu, Numerical investigation of cooling heat transfer to
end heat exchangers.
supercritical CO2 in a horizontal circular tube, J. Supercrit. Fluids 55 (2010)
116e121.
Acknowledgements [15] J. Yang, Y. Oka, Y. Ishiwatari, J. Liu, J. Yoo, Numerical investigation of heat
transfer in upward flow of supercritical water in circular tubes and tight fuel
rod bundles, Nucl. Eng. Des. 237 (2007) 420e430.
The present work was carried out under a project sponsored by [16] P.F. Lisboa, J. Fernandes, P.C. Simoes, J.P.B. Mota, E. Saatdjian, Computational-
Extramural Research Division, Council of Scientific and Industrial fluid-dynamics study of a Kenics static mixer as a heat exchanger for super-
critical carbon dioxide, J. Supercrit. Fluids 55 (2010) 107e115.
Research (CSIR), Government of India. The financial support [17] B.E. Launder, D.B. Spalding, The numerical computation of turbulent flows,
provided by CSIR is gratefully acknowledged. Comput. Methods Appl. Mech. Eng. 3 (1974) 269e289.
[18] P.K. Vijayan, H. Austregesilo, Scaling laws for single-phase natural circulation
loops, Nucl. Eng. Des. 152 (1994) 331e347.
References [19] H.H. Bau, K.E. Torrance, Transient and steady behavior of an open, symmet-
rically heated, free convection loop, Int. J. Heat Mass Transf. 24 (1981)
[1] K. Kiran Kumar, M. Ram Gopal, Carbon dioxide as secondary fluid in natural 597e609.
circulation loops, Proc. IMechE, Part E: J. Process Mech. Engg. 223 (2009) [20] P.K. Vijayan, Experimental observations on the general trends of the steady
189e194. state and stability behaviour of single-phase natural circulation loops, Nucl.
[2] V. Dostal, P. Hejzlar, M.J. Driscoll, The supercritical carbon dioxide power Eng. Des. 215 (2002) 139e152.
cycle: comparison to other advanced power cycles, Nucl. Technol. 154 (2006) [21] S.K. Mousavian, M. Misale, F. D’Auria, M.A. Salehi, Transient and stability
283e301. analysis in single-phase natural circulation, Ann. Nucl. Energy 31 (2004)
[3] P. Bondioli, C. Mariani, E. Mossa, A. Fedelli, A. Muller, Lampante olive oil 1177e1198.
refining with supercritical carbon dioxide, J. Am. Oil Chem. Soc. 69 (1992) [22] NIST Standard Reference Database-REFPROP, Version 8.0 (2006).
477e480. [23] F.W. Dittus, L.M.K. Boelter, Heat Transfer in Automobile Radiators of Tubular
[4] F.C.V.N. Fourie, C.E. Schwarz, J.H. Knoetze, Phase equilibria of alcohols in Type, 2, University of California Publications in Engineering, 1930, 443e461.
supercritical fluids Part I. The effect of the position of the hydroxyl group for [24] V. Gnielinsky, New equation for heat and mass transfer in turbulent pipe and
linear C8 alcohols in supercritical carbon dioxide, J. Supercrit. Fluids 47 (2008) channel flow, Int. Chem. Eng. 16 (1976) 359e368.
161e167. [25] G.K. Filonenko, Hydraulic resistance in pipes. in: E. Schlunder, et al. (Eds.),
[5] H. Yamaguchi, X.R. Zhang, K. Fujima, Basic study on new cryogenic refriger- Heat Exchanger Design Handbook, vol. 1. Teploenergetika, Hemisphere Pub.
ation using CO2 solidegas two phase flow, Int. J. Refrig. 31 (2008) 404e410. Corp., 1954, pp. 40e44.