# Numerical Analysis of General Trends in Single-Phase Natural Circulation in a 2D-Annular Loop

Authors: Gilles Desrayaud; Alberto Fichera
Journal: Science and Technology of Nuclear Installations
 (2008)

---

## Abstract

The aim of this paper is to address fluid flow behavior of natural circulation in a 2D-annular loop filled with water. A two-dimensional, numerical analysis of natural convection in a 2D-annular closed-loop thermosyphon has been performed for various radius ratios from 1.2 to 2.0, the loop being heated at a constant flux over the bottom half and cooled at a constant temperature over the top half. It has been numerically shown that natural convection in a 2D-annular closed-loop thermosyphon is capable of showing pseudoconductive regime at pitchfork bifurcation, stationary convective regimes without and with recirculating regions occurring at the entrance of the exchangers, oscillatory convection at Hopf bifurcation and Lorenz-like chaotic flow. The complexity of the dynamic properties experimentally encountered in toroidal or rectangular loops is thus also found here.


---

## Body

## 1. Introduction

Natural circulation is an important mechanism, and the
knowledge of its behavior is of interest for cooling purposes in industrial
processes, including solar water heaters, geothermal processes, gas turbine
blade cooling, and as part of the emergency core cooling system in nuclear reactors.
Heat removal by natural circulation provides a practical means to cool and
depressurize the reactor primary coolant system after a reactor accident. But
in advanced nuclear reactor designs, heat removal systems driven by natural
convection are also a potentially important design feature. Indeed, some advanced
nuclear plant designs rely on natural circulation to remove core power under
normal operation (startup, normal power operation, and shutdown), and some
designs rely on natural circulation to provide cooling of the containment.
Because of their practical importance, thermosyphons have been the subject of a
large number of theoretical and experimental studies. A review of the wide
applications of natural circulation loops in engineering systems has been given
by Zvirin [1]. These also attract attention because of the variety of fluid motions and the complexity of the dynamic properties encountered, in spite of
the simplicity of their geometry. The research pioneered by Keller [2], Welander [3], and Malkus [4] has been reviewed by Greif [5]. The presence of
a reverse flow region was first qualitatively reported by Creveling et al. [6] who also first observed the Lorenz-like chaotic flow in their experiments (see
also [7–9]).Previous studies of toroidal loops have utilized a one-dimensional
approach by averaging the governing equations over the pipe cross-section,
which required a priori specifications of the friction and the heat transfer
coefficients. To the best of our knowledge, in the literature only one
numerical experiment on steady 3D flow in a toroidal loop exists (Lavine et al.
[10]).The purpose
of the present study is to investigate by direct numerical integration of the
governing equations the steady and unsteady motions in a thermosyphon of simple
well-defined geometry, a 2D-annular loop.


## 2. Governing Equations and Resolution

Consider a 2D-annular loop of laminar
fluid heated over one-half of its area at a uniform heat flux (q″) and cooled over the remaining half at a constant temperature (Tc)
as shown in Figure 1. The inner
(respectively outer) cylinder radius is ri* (resp., ro*), the channel gap being defined as a=ro*−ri*. The
classical governing (Navier-Stokes plus energy) equations for incompressible
fluid with the Boussinesq approximation used in cylindrical coordinates are not
written here due to lack of space but can be found in the work of Desrayaud et
al. [11]. The way in which the equations are nondimensionalized is the following
for the dimension, time, velocity, and temperature:
(1)r=r*−ri*a,τ=αt*a2,u=au*α,Θ=T−TcΔT,ΔT=aq″κ,
The thermophysical characteristics of the fluid are
its thermal diffusivity α, its
kinematic viscosity ν, and
its thermal conductivity κ. The sign β is the coefficient of thermal expansion and
the nondimensional radial and axial velocities are u and v.Figure 1
Coordinates and 2D-annular geometry.The nondimensional boundary conditions are the following,
the angular
coordinateθ being measured from the downward vertical
(see Figure 1):
(2)Heatingsectionat−π2≤θ≤π2,r=0and1,u=v=0,∂Θ∂r=−1,Coolingsectionatπ2≤θ≤3π2,r=0and1,u=v=0,Θ=0,2π-periodicityintheθdirection,0<r<1,u,v,Θ∣θ=0=u,v,Θ∣θ=2π.
The motionless and isothermal solution used as the
initial guess for computations is given by
(3)forτ=0at0≤r≤1,0≤θ≤2π,u=v=0,Θ=0.
The nondimensional parameters that govern the flow are the Rayleigh
number Ra built on the channel gap a, the Prandtl number Pr, and the radius ratio R which are defined by
(4)Ra=gβa3ΔTαv,Pr=vα,R=ro*ri*.
The dimensionless mass flow rate circulating inside
the loop (which is also the dimensionless cross-average velocity) is defined as Q=∫01v(r,θ)dr and is angle independent owing to the laminar fluid incompressibility.The control volume procedure is utilized to discretize
on a staggered, uniform cylindrical grid the nonlinear system of governing
equations and boundary conditions with the second-order centered scheme for the
convective terms. The SIMPLER algorithm is employed for
the velocity-pressure coupling. The momentum and energy equations are cast in
transient form and the time-integration is performed using an alternating
direction implicit (ADI) scheme. Boundary conditions of periodic type were applied atq=0 and
2π. This
consists of using the solutions calculated at a previous radial sweeping as the
boundary conditions in the angular sweeping of the ADI procedure. A uniform grid with 30×320 control
volumes in the r and θ directions was
used to obtain all the results presented in this paper. The numerical code has been
extensively validated on several cases (see [11, 12]).


## 3. Results and Discussion

Over 100 runs
have been carried out for steady flows and 50 for unsteady flows for
radius ratiosR=1.2,1.4,1.6,1.8, and 2.0. The working
fluid was water with Prandtl number, Pr=5, in all cases. Three videos accompany this article (see videos in  Supplementary Material available online at doi:10.1155/2008/895695). These have been compressed in mp4
format and as a result, they are somewhat inferior in quality; however, in the
resulting stylized representation, the important features show up rather well
(although the colors do not correspond to those in the illustrations in the
text).A map of the various regimes exemplifying the “route
to chaos” in such configuration is presented in Figure2 for radius ratios R=2.0,1.8,1.6,1.4, and 1.2. These regimes are pseudoconductive (ps-c), convective (conv), secondary motion (second),
oscillatory (oscill), and chaotic (reverse) regimes. It must be noted that
this map only gives general trends of the flow. No attempt to find accurate
transition Rayleigh numbers has been made because unstationary motions are very
CPU-time consuming. For instance, it is a hard task to properly capture by
brute force integration the Rayleigh number at which the monoperiodic motion
occurs because the Ra-gap is very narrow. Moreover, this regime can degenerate
after long time integration into chaotically reversing motion. In the
following, the results used to illustrate the various regimes are for different
radius ratios but, as shown in Figure 2, all these regimes have been found whatever the radius ratio is.Figure 2
Map of the various regimes: pseudoconductive (ps-c), convective (conv), secondary motion (second), oscillatory (oscill), and chaotic (reverse) regimes.


### 3.1. Steady Flow

#### 3.1.1. Pseudoconductive Regime

At very low values of the Rayleigh number
(Ra=10, R=1.8; see Figure 3), the stratified
temperature which is almost symmetrically distributed in both halves of the
annulus ([0,−π]
and [0,π] in
Figure 3(a)) reveals the importance of axial conduction, with convection
playing only a minor role in the heat transfer. Hence, the heat is mainly transferred by
conduction from the lower hot part to the upper cold part, and thus this regime can
be appropriately qualified as pseudoconductive. A global counter-clockwise
motion all around the loop is almost nonexistent and two very large recirculations
of weak motion occur (see Figure 3(b)). These two large but weak vortices
characterize the pseudoconductive regime, both being located at the heater, one
along the inner wall while the other is along the outer wall. These two
vortices occupied each almost one quarter of the loop, a weak global motion
flowing in-between the vortices. A weak global motion always sets in around the
loop, even at very low Rayleigh numbers, due to the cylindrical geometry of
the loop, and a purely conductive solution was never found
(i.e., with no global fluid motion). This configuration is like a
Rayleigh-Bénard one (heating at the bottom and cooling at the top part) but
contrarily to this well-known case for which no motion is expected until a
critical Rayleigh number is reached, the heated fluid immediately proceeds upward
along the inner bottom heated cylinder because no stable equilibrium can exist
along a convex wall. Furthermore, due to the symmetry of
the heating and cooling sections to the vertical axis, the global motion can be
clockwise or counter-clockwise at random. Two convective solutions of very weak
motion branch off symmetrically from the state of rest (i.e., Ra=0)
undergoing a pitchfork bifurcation [11]. Surprisingly, Figure 2 shows that
decreasing the radius ratio decreases the Rayleigh number at which the pseudoconductive
regime occurs. When the radius ratio decreases, one way in which this can occur
being a decrease in the channel gap between the walls, the size of the cells
becomes too large to permit the fluid to move, the vena contracta being too small. The cells are then pushed away by
the flow and thus disappear.Temperature and stream function fields in a pseudoconductive regime;Ra=10, R=1.8, and Pr=5. (a) Isotherms, from 0.5 (blue) to 8.5 (red) in steps of 0.5. (b) Stream function field, from −0.005 (blue) to 0.1 (red) in steps of 0.005 (Vortex in the red zone: 0.102, 0.104).
(a)(b)


#### 3.1.2. Convective Regime

For moderate values of the Rayleigh number, a
quasi-one-dimensional flow exists along the loop and the flow is steady without
undergoing any oscillatory process. Figure4 shows the temperature and stream
function fields on the two branches of the pitchfork bifurcation for a Rayleigh
number equal to 1000 and a radius ratio R=2.0. One of these two
solutions has been obtained by increasing the Rayleigh number from 0 to 1000
(see Figure 4(a)) while the other has been obtained by increasing the Rayleigh
number step by step (see Figure 4(b)). With increase in Ra, the fluid flow
undergoes a short-term oscillation before becoming stable. The temperature and
stream function fields for Ra=22000 and R=1.60
are shown in Figure 5, with the fluid moving around the loop in a
counter-clockwise direction. The streamlines that were concentric for low
values of the Rayleigh number are now slightly deformed at the entrance region
of the exchangers, while
being more distorted at the heater (see Figure 5(b)). This distortion is due
to a zone of a very strong temperature increase clearly visible at the entrance
of the heater along the outer wall (see Figure 5(a)). These distorted
streamlines are less pronounced at the cooler entrance due to the imposed
temperature at the walls.Temperature and stream function fields in a convective regime;Ra=1000, R=2.0, and Pr=5.
(a)
Counter-clockwise motion; isotherms, from 0.1 (blue) to 1.1 (red) in steps of 0.1(b)
Clockwise motion; stream function field, from 1 to 11 in steps of 1Steady motion in a convective regime;Ra=22000, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.55 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)


#### 3.1.3. Secondary Flow Regime

AtRa=22500 and R=1.60 (see
Figure 6), it should be noted that, in addition to the circulatory main flow,
two cells in the first and in the third quadrants can be seen. These two recirculating
regions occur near the outer wall of the entrance of the exchangers (see Figure 6(b)),
with the bigger vortex always being at the heater while the one at the cooler
is very weak. As a consequence, the zone of the strong temperature increase
that extended largely downstream along the outer walls is now reduced (see Figure 6(a)).Steady motion in a secondary flow regime;Ra=22500, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.65 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)The dimensionless temperatures at the inner and outer walls are shown in
Figure7 only along the heater, the cooler being at an imposed temperature. The
temperature of the outer wall is always largely greater than the temperature of
the inner wall. This is due to the boundary condition, imposed heat flux at the
heater, combined with the curvature. Since the surface of the outer wall is R
times greater than the surface of the inner wall, the total heat flux
transferred to the convective fluid flow is also R times greater and there
results a higher level of temperature at the outer wall. At
Ra=22000 (22 500, resp.), the average temperature of the
outer wall is Θ¯h=0.496 (resp., 0.527) while it is only 0.390 (resp., 0.394) at the inner
wall. Following the main circulation from θ=−90° to θ=90°, Figure 7 (Ra=22000) clearly shows at the entrance of
the heater a sharp temperature increase at the outer wall followed by a slight
decrease of the temperature. This inverse temperature gradient explains why the
fluid is able to proceed upstream along the outer wall at the entrance of the
heat exchanger. At Ra=22500, a very large
recirculation occurs at the entrance of the heater (see Figure 6(b)) amplifying the strong inverse
temperature gradient (see Figure 7). Moreover, no increase of the average
temperature at the inner wall is noticed from 22 000 to 22 500 while
a large one occurs at the outer wall.Figure 7
Inner and outer walls temperature at the heater;R=1.6, Pr=5.The
secondary cellular structures appear at values of the Rayleigh number that are all high since the radius
ratio is small, and are representative of the upper limit of the steady motion
(see Figure2). Decreasing the radius ratio has a very strong stabilizing effect, with
the promotion of the secondary cells thus becoming very difficult for narrow
channel gaps because a consequence of the occurrence of the cells is always a
decrease of the mass flow rate, as can be seen in Figure 8(a). There is a
competition between the strengthening of cells due to the temperature gradient
and the main circulation of the flow which contains and limits the cells’
growth.Oscillatory regime.
(a)
Time history of the mass flow rate;Ra=23500, R=1.6, and Pr=5(b)
Limit cycle;Ra=24600, R=1.6, and Pr=5


### 3.2. Oscillatory Regime

Time history of the mass flow rate is presented in Figure8(a) from
Ra=22000 to 23 500 for a radius ratio,
R=1.6. At Ra=22000 and for this radius ratio, the
circulation flow is unidirectional (see Figure 5(b)) while secondary cells are
depicted at Ra=22500 (see Figure 6(b)). After initial damped
oscillations, the mass flow rate abruptly decreases owing to the occurrence of
secondary cells along the outer wall at the entrance of the exchangers. This is
followed by small sustainable oscillations of the secondary cells on the spot, with their strength (and size) varying slightly. The flow
structure behavior of the stream function field during the same time history as
in Figure 8(a) is well documented in Video 1 (see Supplementary Material), clearly showing the occurrence of the secondary
cells followed by small oscillatory motion of these cells. At
Ra=23500, the amplitude of the periodic motion is very small
and the oscillations are very slowly damped (not distinguishable in Figure 8(a)
or in Video 1 (see Supplementary Material)).
This is why the limit cycle to which the system is attracted is drawn for
Ra=24600. At this Rayleigh number, the limit
cycle illustrated in Figure 8(b) indicates that the motion is (almost)
periodic with one frequency, f=6.66±0.39. The phase
portrait of [v(0.5,0),Θ(0.5,0)] is for a sampling interval 0.0025 during a time
interval 1 representing 6 cycles with almost 310 iterations per cycle.


### 3.3. Chaotically Reversing Motion

Figure9(a) presenting a very long time history of the mass flow rate
from Ra=45000 to 47 000 for R=1.4 clearly
shows the aperiodic behavior of the flow. The nondimensional time interval is
equal to 35 which represents 70 000 iterations. It must be remarked that
the flow reversal process is only reached after a long time equal to 15 during
which the flow rate oscillates with increasing amplitude until it reverses
direction, whereupon oscillations initiate in a new flow direction. This
variable clearly indicates the changes of flow direction with its change of
sign. These reversals in flow directions are also accompanied by large
temperature variations which can be detrimental to the safety operation of heat
removal systems. Figure 9(b) presents, for Ra=15000 and
R=2.0, an enlargement of the time history of the axial component of
the velocity at the middle of the cooler (0.5,π) during a
short period of time which presents three oscillations followed by four reverse
flows. Figures 10 and 11 give eight snapshots of the temperature and of the
stream function fields equally time-spaced (0.048) in clockwise progression
around the figure at the center during one reverse flow process. The time locations of these snapshots are also
indicated as filled red squares in Figure 9(b). The flow direction is clearly
distinguishable on the temperature fields, with the isotherms following the
flow direction (see Figure 10). Relatively to the flow field, the recirculating
regions at the entrance of the exchangers shift their location at each flow
reversal occupying the first and third quadrants for counter-clockwise motion (see Figure 11 snapshots 8,1) and the second and fourth quadrants for clockwise motion
(see Figure 11 snapshots 4,5). In the center of Figures 10 and 11, we also show the average
temperature and stream function fields which are obtained by adding all the
instantaneous fields at each time step during the two quasiperiodic flow
evolutions showing reverse flows during the time [35.672,36.440] (see Figure 9(b)).
The average fields are not exactly symmetrical with respect to the vertical
centerplane because the flow evolution is not periodic. Videos 2 and 3 (see Supplementary Material), respectively, present the
behavior of the temperature and stream function fields during the same time
history as in Figures 9(b) and usefully complement Figures 10 and 11. The first
three oscillations of the temperature and stream function fields and the way in
which the four reverse flow directions occur are clearly seen.Chaotically reversing regime.
(a)
Time history;Ra=47000, R=1.4, and Pr=5(b)
Enlargement of the time history;Ra=15000, R=2.0, and Pr=5Equally spaced snapshots of the temperature;Ra=15000, R=2.0, and Pr=5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Equally spaced snapshots of the stream function;Ra=15000, R=2.0, and Pr=5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Several experiments on thermal convection in a closed
loop have been made. Although they were conducted in toroidal
[6, 13, 14] or rectangular [7–9, 15] loops and not
in an annular loop, some qualitative comparisons can nevertheless be made.
First, Stern et al. [14], using the same boundary conditions as in the present
numerical study but for a toroidal loop, experimentally demonstrated the flow
reversal in the heating section to be stronger than the reversal in the cooling
section. Moreover, all these experiments
[6–9, 14, 15] showed that after steady state,
oscillatory motions occurred followed by reverse flows as in the present
numerical experiment. Finally, some of them also experienced Lorenz-like
behavior [6–9] as in Figure 12 where the
recognizable shape of the Lorenz attractor can be seen. The attractor projected
in this plane resembles a butterfly. This phase
portrait has been built at Ra=48000 for R=1.4
with only one time series (horizontal axis), that of the mass flow rate and
using two time step lags of Δτ=0.04
and 0.06 for the other axis and during a total dimensionless time equal to 10.
Only one iterative point (in red) in ten is used to draw the Lorenz attractor.
In Figure 12, there are only two excursions of the flow around the positive
wing (the motion thus being in the counter-clockwise direction) while most of
them swirl over the negative wing (the motion being in the clockwise
direction). No better description of the phenomenon can be given than the
following: “As the Lorenz attractor
is plotted, a strand will be drawn from one point, and will start weaving the
outline of the right butterfly wing. Then it swirls over to the left wing and
draws its center. The attractor will continue weaving back and forth between
the two wings, its motion seemingly random, its very action mirroring the chaos
which drives the process.” (http://zeuscat.com/andrew/chaos/lorenz.html).
For those who wish to observe real-time Lorenz attractor movies, many examples
can easily be found on the Web.Figure 12
Lorenz-like attractor;Ra=48000, R=1.4, and Pr=5.


## 3.1. Steady Flow

### 3.1.1. Pseudoconductive Regime

At very low values of the Rayleigh number
(Ra=10, R=1.8; see Figure 3), the stratified
temperature which is almost symmetrically distributed in both halves of the
annulus ([0,−π]
and [0,π] in
Figure 3(a)) reveals the importance of axial conduction, with convection
playing only a minor role in the heat transfer. Hence, the heat is mainly transferred by
conduction from the lower hot part to the upper cold part, and thus this regime can
be appropriately qualified as pseudoconductive. A global counter-clockwise
motion all around the loop is almost nonexistent and two very large recirculations
of weak motion occur (see Figure 3(b)). These two large but weak vortices
characterize the pseudoconductive regime, both being located at the heater, one
along the inner wall while the other is along the outer wall. These two
vortices occupied each almost one quarter of the loop, a weak global motion
flowing in-between the vortices. A weak global motion always sets in around the
loop, even at very low Rayleigh numbers, due to the cylindrical geometry of
the loop, and a purely conductive solution was never found
(i.e., with no global fluid motion). This configuration is like a
Rayleigh-Bénard one (heating at the bottom and cooling at the top part) but
contrarily to this well-known case for which no motion is expected until a
critical Rayleigh number is reached, the heated fluid immediately proceeds upward
along the inner bottom heated cylinder because no stable equilibrium can exist
along a convex wall. Furthermore, due to the symmetry of
the heating and cooling sections to the vertical axis, the global motion can be
clockwise or counter-clockwise at random. Two convective solutions of very weak
motion branch off symmetrically from the state of rest (i.e., Ra=0)
undergoing a pitchfork bifurcation [11]. Surprisingly, Figure 2 shows that
decreasing the radius ratio decreases the Rayleigh number at which the pseudoconductive
regime occurs. When the radius ratio decreases, one way in which this can occur
being a decrease in the channel gap between the walls, the size of the cells
becomes too large to permit the fluid to move, the vena contracta being too small. The cells are then pushed away by
the flow and thus disappear.Temperature and stream function fields in a pseudoconductive regime;Ra=10, R=1.8, and Pr=5. (a) Isotherms, from 0.5 (blue) to 8.5 (red) in steps of 0.5. (b) Stream function field, from −0.005 (blue) to 0.1 (red) in steps of 0.005 (Vortex in the red zone: 0.102, 0.104).
(a)(b)


### 3.1.2. Convective Regime

For moderate values of the Rayleigh number, a
quasi-one-dimensional flow exists along the loop and the flow is steady without
undergoing any oscillatory process. Figure4 shows the temperature and stream
function fields on the two branches of the pitchfork bifurcation for a Rayleigh
number equal to 1000 and a radius ratio R=2.0. One of these two
solutions has been obtained by increasing the Rayleigh number from 0 to 1000
(see Figure 4(a)) while the other has been obtained by increasing the Rayleigh
number step by step (see Figure 4(b)). With increase in Ra, the fluid flow
undergoes a short-term oscillation before becoming stable. The temperature and
stream function fields for Ra=22000 and R=1.60
are shown in Figure 5, with the fluid moving around the loop in a
counter-clockwise direction. The streamlines that were concentric for low
values of the Rayleigh number are now slightly deformed at the entrance region
of the exchangers, while
being more distorted at the heater (see Figure 5(b)). This distortion is due
to a zone of a very strong temperature increase clearly visible at the entrance
of the heater along the outer wall (see Figure 5(a)). These distorted
streamlines are less pronounced at the cooler entrance due to the imposed
temperature at the walls.Temperature and stream function fields in a convective regime;Ra=1000, R=2.0, and Pr=5.
(a)
Counter-clockwise motion; isotherms, from 0.1 (blue) to 1.1 (red) in steps of 0.1(b)
Clockwise motion; stream function field, from 1 to 11 in steps of 1Steady motion in a convective regime;Ra=22000, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.55 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)


### 3.1.3. Secondary Flow Regime

AtRa=22500 and R=1.60 (see
Figure 6), it should be noted that, in addition to the circulatory main flow,
two cells in the first and in the third quadrants can be seen. These two recirculating
regions occur near the outer wall of the entrance of the exchangers (see Figure 6(b)),
with the bigger vortex always being at the heater while the one at the cooler
is very weak. As a consequence, the zone of the strong temperature increase
that extended largely downstream along the outer walls is now reduced (see Figure 6(a)).Steady motion in a secondary flow regime;Ra=22500, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.65 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)The dimensionless temperatures at the inner and outer walls are shown in
Figure7 only along the heater, the cooler being at an imposed temperature. The
temperature of the outer wall is always largely greater than the temperature of
the inner wall. This is due to the boundary condition, imposed heat flux at the
heater, combined with the curvature. Since the surface of the outer wall is R
times greater than the surface of the inner wall, the total heat flux
transferred to the convective fluid flow is also R times greater and there
results a higher level of temperature at the outer wall. At
Ra=22000 (22 500, resp.), the average temperature of the
outer wall is Θ¯h=0.496 (resp., 0.527) while it is only 0.390 (resp., 0.394) at the inner
wall. Following the main circulation from θ=−90° to θ=90°, Figure 7 (Ra=22000) clearly shows at the entrance of
the heater a sharp temperature increase at the outer wall followed by a slight
decrease of the temperature. This inverse temperature gradient explains why the
fluid is able to proceed upstream along the outer wall at the entrance of the
heat exchanger. At Ra=22500, a very large
recirculation occurs at the entrance of the heater (see Figure 6(b)) amplifying the strong inverse
temperature gradient (see Figure 7). Moreover, no increase of the average
temperature at the inner wall is noticed from 22 000 to 22 500 while
a large one occurs at the outer wall.Figure 7
Inner and outer walls temperature at the heater;R=1.6, Pr=5.The
secondary cellular structures appear at values of the Rayleigh number that are all high since the radius
ratio is small, and are representative of the upper limit of the steady motion
(see Figure2). Decreasing the radius ratio has a very strong stabilizing effect, with
the promotion of the secondary cells thus becoming very difficult for narrow
channel gaps because a consequence of the occurrence of the cells is always a
decrease of the mass flow rate, as can be seen in Figure 8(a). There is a
competition between the strengthening of cells due to the temperature gradient
and the main circulation of the flow which contains and limits the cells’
growth.Oscillatory regime.
(a)
Time history of the mass flow rate;Ra=23500, R=1.6, and Pr=5(b)
Limit cycle;Ra=24600, R=1.6, and Pr=5


## 3.1.1. Pseudoconductive Regime

At very low values of the Rayleigh number
(Ra=10, R=1.8; see Figure 3), the stratified
temperature which is almost symmetrically distributed in both halves of the
annulus ([0,−π]
and [0,π] in
Figure 3(a)) reveals the importance of axial conduction, with convection
playing only a minor role in the heat transfer. Hence, the heat is mainly transferred by
conduction from the lower hot part to the upper cold part, and thus this regime can
be appropriately qualified as pseudoconductive. A global counter-clockwise
motion all around the loop is almost nonexistent and two very large recirculations
of weak motion occur (see Figure 3(b)). These two large but weak vortices
characterize the pseudoconductive regime, both being located at the heater, one
along the inner wall while the other is along the outer wall. These two
vortices occupied each almost one quarter of the loop, a weak global motion
flowing in-between the vortices. A weak global motion always sets in around the
loop, even at very low Rayleigh numbers, due to the cylindrical geometry of
the loop, and a purely conductive solution was never found
(i.e., with no global fluid motion). This configuration is like a
Rayleigh-Bénard one (heating at the bottom and cooling at the top part) but
contrarily to this well-known case for which no motion is expected until a
critical Rayleigh number is reached, the heated fluid immediately proceeds upward
along the inner bottom heated cylinder because no stable equilibrium can exist
along a convex wall. Furthermore, due to the symmetry of
the heating and cooling sections to the vertical axis, the global motion can be
clockwise or counter-clockwise at random. Two convective solutions of very weak
motion branch off symmetrically from the state of rest (i.e., Ra=0)
undergoing a pitchfork bifurcation [11]. Surprisingly, Figure 2 shows that
decreasing the radius ratio decreases the Rayleigh number at which the pseudoconductive
regime occurs. When the radius ratio decreases, one way in which this can occur
being a decrease in the channel gap between the walls, the size of the cells
becomes too large to permit the fluid to move, the vena contracta being too small. The cells are then pushed away by
the flow and thus disappear.Temperature and stream function fields in a pseudoconductive regime;Ra=10, R=1.8, and Pr=5. (a) Isotherms, from 0.5 (blue) to 8.5 (red) in steps of 0.5. (b) Stream function field, from −0.005 (blue) to 0.1 (red) in steps of 0.005 (Vortex in the red zone: 0.102, 0.104).
(a)(b)


## 3.1.2. Convective Regime

For moderate values of the Rayleigh number, a
quasi-one-dimensional flow exists along the loop and the flow is steady without
undergoing any oscillatory process. Figure4 shows the temperature and stream
function fields on the two branches of the pitchfork bifurcation for a Rayleigh
number equal to 1000 and a radius ratio R=2.0. One of these two
solutions has been obtained by increasing the Rayleigh number from 0 to 1000
(see Figure 4(a)) while the other has been obtained by increasing the Rayleigh
number step by step (see Figure 4(b)). With increase in Ra, the fluid flow
undergoes a short-term oscillation before becoming stable. The temperature and
stream function fields for Ra=22000 and R=1.60
are shown in Figure 5, with the fluid moving around the loop in a
counter-clockwise direction. The streamlines that were concentric for low
values of the Rayleigh number are now slightly deformed at the entrance region
of the exchangers, while
being more distorted at the heater (see Figure 5(b)). This distortion is due
to a zone of a very strong temperature increase clearly visible at the entrance
of the heater along the outer wall (see Figure 5(a)). These distorted
streamlines are less pronounced at the cooler entrance due to the imposed
temperature at the walls.Temperature and stream function fields in a convective regime;Ra=1000, R=2.0, and Pr=5.
(a)
Counter-clockwise motion; isotherms, from 0.1 (blue) to 1.1 (red) in steps of 0.1(b)
Clockwise motion; stream function field, from 1 to 11 in steps of 1Steady motion in a convective regime;Ra=22000, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.55 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)


## 3.1.3. Secondary Flow Regime

AtRa=22500 and R=1.60 (see
Figure 6), it should be noted that, in addition to the circulatory main flow,
two cells in the first and in the third quadrants can be seen. These two recirculating
regions occur near the outer wall of the entrance of the exchangers (see Figure 6(b)),
with the bigger vortex always being at the heater while the one at the cooler
is very weak. As a consequence, the zone of the strong temperature increase
that extended largely downstream along the outer walls is now reduced (see Figure 6(a)).Steady motion in a secondary flow regime;Ra=22500, R=1.6, and Pr=5. (a) Isotherms, from 0.05 (blue) to 0.65 (red) in steps of 0.05. (b) Stream function field, from 5 (blue) to 60 (red) in steps of 5.
(a)(b)The dimensionless temperatures at the inner and outer walls are shown in
Figure7 only along the heater, the cooler being at an imposed temperature. The
temperature of the outer wall is always largely greater than the temperature of
the inner wall. This is due to the boundary condition, imposed heat flux at the
heater, combined with the curvature. Since the surface of the outer wall is R
times greater than the surface of the inner wall, the total heat flux
transferred to the convective fluid flow is also R times greater and there
results a higher level of temperature at the outer wall. At
Ra=22000 (22 500, resp.), the average temperature of the
outer wall is Θ¯h=0.496 (resp., 0.527) while it is only 0.390 (resp., 0.394) at the inner
wall. Following the main circulation from θ=−90° to θ=90°, Figure 7 (Ra=22000) clearly shows at the entrance of
the heater a sharp temperature increase at the outer wall followed by a slight
decrease of the temperature. This inverse temperature gradient explains why the
fluid is able to proceed upstream along the outer wall at the entrance of the
heat exchanger. At Ra=22500, a very large
recirculation occurs at the entrance of the heater (see Figure 6(b)) amplifying the strong inverse
temperature gradient (see Figure 7). Moreover, no increase of the average
temperature at the inner wall is noticed from 22 000 to 22 500 while
a large one occurs at the outer wall.Figure 7
Inner and outer walls temperature at the heater;R=1.6, Pr=5.The
secondary cellular structures appear at values of the Rayleigh number that are all high since the radius
ratio is small, and are representative of the upper limit of the steady motion
(see Figure2). Decreasing the radius ratio has a very strong stabilizing effect, with
the promotion of the secondary cells thus becoming very difficult for narrow
channel gaps because a consequence of the occurrence of the cells is always a
decrease of the mass flow rate, as can be seen in Figure 8(a). There is a
competition between the strengthening of cells due to the temperature gradient
and the main circulation of the flow which contains and limits the cells’
growth.Oscillatory regime.
(a)
Time history of the mass flow rate;Ra=23500, R=1.6, and Pr=5(b)
Limit cycle;Ra=24600, R=1.6, and Pr=5


## 3.2. Oscillatory Regime

Time history of the mass flow rate is presented in Figure8(a) from
Ra=22000 to 23 500 for a radius ratio,
R=1.6. At Ra=22000 and for this radius ratio, the
circulation flow is unidirectional (see Figure 5(b)) while secondary cells are
depicted at Ra=22500 (see Figure 6(b)). After initial damped
oscillations, the mass flow rate abruptly decreases owing to the occurrence of
secondary cells along the outer wall at the entrance of the exchangers. This is
followed by small sustainable oscillations of the secondary cells on the spot, with their strength (and size) varying slightly. The flow
structure behavior of the stream function field during the same time history as
in Figure 8(a) is well documented in Video 1 (see Supplementary Material), clearly showing the occurrence of the secondary
cells followed by small oscillatory motion of these cells. At
Ra=23500, the amplitude of the periodic motion is very small
and the oscillations are very slowly damped (not distinguishable in Figure 8(a)
or in Video 1 (see Supplementary Material)).
This is why the limit cycle to which the system is attracted is drawn for
Ra=24600. At this Rayleigh number, the limit
cycle illustrated in Figure 8(b) indicates that the motion is (almost)
periodic with one frequency, f=6.66±0.39. The phase
portrait of [v(0.5,0),Θ(0.5,0)] is for a sampling interval 0.0025 during a time
interval 1 representing 6 cycles with almost 310 iterations per cycle.


## 3.3. Chaotically Reversing Motion

Figure9(a) presenting a very long time history of the mass flow rate
from Ra=45000 to 47 000 for R=1.4 clearly
shows the aperiodic behavior of the flow. The nondimensional time interval is
equal to 35 which represents 70 000 iterations. It must be remarked that
the flow reversal process is only reached after a long time equal to 15 during
which the flow rate oscillates with increasing amplitude until it reverses
direction, whereupon oscillations initiate in a new flow direction. This
variable clearly indicates the changes of flow direction with its change of
sign. These reversals in flow directions are also accompanied by large
temperature variations which can be detrimental to the safety operation of heat
removal systems. Figure 9(b) presents, for Ra=15000 and
R=2.0, an enlargement of the time history of the axial component of
the velocity at the middle of the cooler (0.5,π) during a
short period of time which presents three oscillations followed by four reverse
flows. Figures 10 and 11 give eight snapshots of the temperature and of the
stream function fields equally time-spaced (0.048) in clockwise progression
around the figure at the center during one reverse flow process. The time locations of these snapshots are also
indicated as filled red squares in Figure 9(b). The flow direction is clearly
distinguishable on the temperature fields, with the isotherms following the
flow direction (see Figure 10). Relatively to the flow field, the recirculating
regions at the entrance of the exchangers shift their location at each flow
reversal occupying the first and third quadrants for counter-clockwise motion (see Figure 11 snapshots 8,1) and the second and fourth quadrants for clockwise motion
(see Figure 11 snapshots 4,5). In the center of Figures 10 and 11, we also show the average
temperature and stream function fields which are obtained by adding all the
instantaneous fields at each time step during the two quasiperiodic flow
evolutions showing reverse flows during the time [35.672,36.440] (see Figure 9(b)).
The average fields are not exactly symmetrical with respect to the vertical
centerplane because the flow evolution is not periodic. Videos 2 and 3 (see Supplementary Material), respectively, present the
behavior of the temperature and stream function fields during the same time
history as in Figures 9(b) and usefully complement Figures 10 and 11. The first
three oscillations of the temperature and stream function fields and the way in
which the four reverse flow directions occur are clearly seen.Chaotically reversing regime.
(a)
Time history;Ra=47000, R=1.4, and Pr=5(b)
Enlargement of the time history;Ra=15000, R=2.0, and Pr=5Equally spaced snapshots of the temperature;Ra=15000, R=2.0, and Pr=5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Equally spaced snapshots of the stream function;Ra=15000, R=2.0, and Pr=5.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Several experiments on thermal convection in a closed
loop have been made. Although they were conducted in toroidal
[6, 13, 14] or rectangular [7–9, 15] loops and not
in an annular loop, some qualitative comparisons can nevertheless be made.
First, Stern et al. [14], using the same boundary conditions as in the present
numerical study but for a toroidal loop, experimentally demonstrated the flow
reversal in the heating section to be stronger than the reversal in the cooling
section. Moreover, all these experiments
[6–9, 14, 15] showed that after steady state,
oscillatory motions occurred followed by reverse flows as in the present
numerical experiment. Finally, some of them also experienced Lorenz-like
behavior [6–9] as in Figure 12 where the
recognizable shape of the Lorenz attractor can be seen. The attractor projected
in this plane resembles a butterfly. This phase
portrait has been built at Ra=48000 for R=1.4
with only one time series (horizontal axis), that of the mass flow rate and
using two time step lags of Δτ=0.04
and 0.06 for the other axis and during a total dimensionless time equal to 10.
Only one iterative point (in red) in ten is used to draw the Lorenz attractor.
In Figure 12, there are only two excursions of the flow around the positive
wing (the motion thus being in the counter-clockwise direction) while most of
them swirl over the negative wing (the motion being in the clockwise
direction). No better description of the phenomenon can be given than the
following: “As the Lorenz attractor
is plotted, a strand will be drawn from one point, and will start weaving the
outline of the right butterfly wing. Then it swirls over to the left wing and
draws its center. The attractor will continue weaving back and forth between
the two wings, its motion seemingly random, its very action mirroring the chaos
which drives the process.” (http://zeuscat.com/andrew/chaos/lorenz.html).
For those who wish to observe real-time Lorenz attractor movies, many examples
can easily be found on the Web.Figure 12
Lorenz-like attractor;Ra=48000, R=1.4, and Pr=5.


## 4. Conclusion

The behavior
of a natural circulation 2D-annular loop that is heated uniformly over the
lower half and cooled by maintaining a constant wall temperature over the upper
half has been numerically investigated. Detailed numerical results for the flow
in the thermosyphon are presented for various radius ratios and all these
findings have been demonstrated by brute force integration of the complete
Navier-Stokes equations.It
has been numerically demonstrated that the complexity of the dynamic properties
experimentally encountered in toroidal and rectangular loops is also found here
as follows:(i)
pseudoconductive
regime for which heat is mainly transferred by conduction from the heater to
the cooler;(ii)
steady flow
showing a global motion all around the loop with heat being convected by the
fluid from the heater to the cooler. The motion can be clockwise or
anticlockwise according to chance due to a Pitchfork bifurcation;(iii)
steady flow
with recirculating regions showing a slowing down of the mass flow rate due to
the occurrence of a vortex at the entrance of the heater reducing the “vena contracta;”(iv)
periodic
motion during which the cells oscillate on
the spot with a constant amplitude of very small magnitude;(v)
Lorenz-like chaotic
flow. The flow rate oscillates with
increasing amplitude until it eventually reverses direction, whereupon
oscillations initiate in a new flow direction.
This chaotic regime appears to be reached through a Ruelle-Takens
scenario involving a sequence of Hopf bifurcations: from a fixed point to a
periodic orbit, hence it is likely that the system possesses a strange
attractor.



---

*Source: 895695-2008-04-29.xml*