Science and Technology of Nuclear Installations

Volume 2008, Article ID 895695, 10 pages

http://dx.doi.org/10.1155/2008/895695

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

^{1}INSSET/LETEM, Université de Picardie Jules Verne, 48 rue Raspail, BP 422, 02109 Saint-Quentin, , France^{2}Dipartimento di Ingegneria Industriale e Meccanica, Università di Catania, Viale A. Doria 6, 95125 Catania, , Italy

Received 16 May 2007; Accepted 23 November 2007

Academic Editor: John Cleveland

Copyright © 2008 Gilles Desrayaud and Alberto Fichera. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 () and cooled over the remaining half at a constant temperature () as shown in Figure 1. The inner (respectively outer) cylinder radius is (resp., ), the channel gap being defined as . 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: 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 and .

The nondimensional boundary conditions are the following, the angular coordinate being measured from the downward vertical (see Figure 1): The motionless and isothermal solution used as the initial guess for computations is given by The nondimensional parameters that govern the flow are the Rayleigh number Ra built on the channel gap , the Prandtl number , and the radius ratio which are defined by The dimensionless mass flow rate circulating inside the loop (which is also the dimensionless cross-average velocity) is defined as 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 at and . 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 control volumes in the 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 ratios and The working fluid was water with Prandtl number, , in all cases. Three videos accompany this article (see videos in Supplementary Material available online at http://dx.doi.org/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 Figure 2 for radius ratios and . These regimes are pseudoconductive (ps-), 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.

##### 3.1. Steady Flow

###### 3.1.1. Pseudoconductive Regime

At very low values of the Rayleigh number
(, ; see Figure 3), the stratified
temperature which is almost symmetrically distributed in both halves of the
annulus (
and 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., )
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.

###### 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. Figure 4 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 . 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 and 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.

###### 3.1.3. Secondary Flow Regime

At and (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)).

The dimensionless temperatures at the inner and outer walls are shown in Figure 7 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 times greater than the surface of the inner wall, the total heat flux transferred to the convective fluid flow is also times greater and there results a higher level of temperature at the outer wall. At (22 500, resp.), the average temperature of the outer wall is (resp., 0.527) while it is only 0.390 (resp., 0.394) at the inner wall. Following the main circulation from to , Figure 7 () 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 , 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.

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 Figure 2). 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.

##### 3.2. Oscillatory Regime

Time history of the mass flow rate is presented in Figure 8(a) from to 23 500 for a radius ratio, . At and for this radius ratio, the circulation flow is unidirectional (see Figure 5(b)) while secondary cells are depicted at (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 , 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 . At this Rayleigh number, the limit cycle illustrated in Figure 8(b) indicates that the motion is (almost) periodic with one frequency, . The phase portrait of 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

Figure 9(a) presenting a very long time history of the mass flow rate from to 47 000 for 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 and , an enlargement of the time history of the axial component of the velocity at the middle of the cooler 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 (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.

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 for with only one time series (horizontal axis), that of the mass flow rate and using two time step lags of 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.

#### 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.

*Nomenclature*

: | Gap width , m |

: | Acceleration due to gravity, m/s^{2} |

: | Prandtl number |

: | Specific heat flux at the heater, W/m^{2} |

: | Dimensionless mass flow rate |

: | Dimensionless radial coordinate |

: | Inner cylinder radius, m |

: | Outer cylinder radius, m |

: | Radius ratio |

Ra: | Rayleigh number () |

: | Dimensional temperature, K |

: | Dimensional cold temperature, K |

: | Dimensionless radial and axial velocity |

*Greek*: | Thermal diffusivity of the fluid, m^{2}/s |

: | Coefficient of thermal expansion of the fluid, K^{-1} |

: | Thermal conductivity of the fluid, W/(mK) |

: | Kinematic viscosity of the fluid, m^{2}/s |

: | Dimensionless time |

: | Angular coordinate measured from downward vertical |

: | Dimensionless temperature |

: | Temperature difference , K |

*Superscript**: | Variable with dimension |

*Subscripts*: | Cold exchanger |

: | Hot exchanger |

: | Inner cylinder |

: | Outer cylinder |

: | At the walls |

#### Acknowledgments

This work was supported by Research Grant no. 031265 from the French National Institute for Advances in Scientific Computations (IDRIS-Computer Center). The authors acknowledge the help of R. Spampinato with some of the computational work.

#### References

- Y. Zvirin, “A review of natural circulation loops in pressurized water reactors and other systems,”
*Nuclear Engineering and Design*, vol. 67, no. 2, pp. 203–225, 1982. View at Publisher · View at Google Scholar - J. Keller, “Periodic oscillations in a model of thermal convection,”
*Journal of Fluid Mechanics*, vol. 26, part 1, pp. 599–606, 1966. View at Publisher · View at Google Scholar - P. Welander, “On the oscillatory instability of a differentially heated fluid loop,”
*Journal of Fluid Mechanics*, vol. 29, part 1, pp. 17–30, 1967. View at Publisher · View at Google Scholar - W. R. V. Malkus, “Non-periodic convection at high and low Prandtl number,”
*Mémoires Société Royale de Sciences de Liège*, vol. 4, pp. 125–128, 1972. View at Google Scholar - R. Greif, “Natural circulation loops,”
*Journal of Heat Transfer*, vol. 110, pp. 1243–1258, 1988. View at Google Scholar - H. F. Creveling, J. F. de Paz, J. Y. Baladi, and R. J. Schoenhals, “Stability characteristics of a single-phase free convection loop,”
*Journal of Fluid Mechanics*, vol. 67, part 1, pp. 65–84, 1975. View at Publisher · View at Google Scholar - G. Cammarata, G. Desrayaud, A. Fichera, and A. Pagano, “An ordinary differential model for rectangular natural circulation loops,” in
*Proceedings of the 12th International Heat Transfer Conference*, vol. 2, pp. 273–278, Grenoble, France, August 2002. - A. Fichera and A. Pagano, “Modelling and control of rectangular natural circulation loops,”
*International Journal of Heat and Mass Transfer*, vol. 46, no. 13, pp. 2425–2444, 2003. View at Publisher · View at Google Scholar - M. Misale, M. Frogheri, P. Ruffino, and F. D'Auria, “Steady state and stability behavior of a single-phase natural circulation loop,” in
*Proceedings of the 11th International Heat Transfer Conference*, vol. 3, pp. 385–390, Kyongju, Korea, August 1998. - A. S. Lavine, R. Greif, and J. A. Humphrey, “Three-dimensional analysis of natural convection in a toroidal loop—the effect of Grashof number,”
*International Journal of Heat and Mass Transfer*, vol. 30, no. 2, pp. 251–262, 1987. View at Publisher · View at Google Scholar - G. Desrayaud, A. Fichera, and M. Marcoux, “Numerical investigation of natural circulation in a 2D-annular closed-loop thermosyphon,”
*International Journal of Heat and Fluid Flow*, vol. 27, no. 1, pp. 154–166, 2006. View at Publisher · View at Google Scholar - P. Cadiou, G. Desrayaud, and G. Lauriat, “Natural convection in a narrow horizontal annulus: the effects of thermal and hydrodynamic instabilities,”
*Journal of Heat Transfer*, vol. 120, no. 4, pp. 1019–1026, 1998. View at Google Scholar - O. Sano, “Cellular structure in a natural convection loop and its chaotic behavior. I . Experiment,”
*Fluid Dynamics Research*, vol. 8, no. 5-6, pp. 189–204, 1991. View at Publisher · View at Google Scholar - C. H. Stern, R. Greif, and J. A. Humphrey, “An experimental study of natural convection in a toroidal loop,”
*Journal of Heat Transfer*, vol. 110, pp. 877–884, 1988. View at Google Scholar - B. J. Huang and R. Zelaya, “Heat transfer behaviour of a rectangular thermosyphon,”
*Journal of Heat Transfer*, vol. 110, pp. 487–493, 1988. View at Google Scholar