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.
(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. 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.
(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 1
(a)
(b)
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)).
(a)
(b)
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.
| (a) Time history of the mass flow rate; , , and |
| (b) Limit cycle; , , and |
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.
| (a) Time history; , , and |
| (b) Enlargement of the time history; , , and |
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(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 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/s2 |
| : | Prandtl number |
| : | Specific heat flux at the heater, W/m2 |
| : | 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 |
| : | Thermal diffusivity of the fluid, m2/s |
| : | Coefficient of thermal expansion of the fluid, K-1 |
| : | Thermal conductivity of the fluid, W/(mK) |
| : | Kinematic viscosity of the fluid, m2/s |
| : | Dimensionless time |
| : | Angular coordinate measured from downward vertical |
| : | Dimensionless temperature |
| : | Temperature difference , K |
| *: | Variable with dimension |
| : | 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.
Supplementary Materials
Video 1 shows the flow structure evolution of the stream function field from Ra = 22000 to Ra = 23500 at R = 1.6 during the same time history as in Figure 7a. Video 2 shows the flow structure evolution of the temperature field during the same time history as in Figure 9. Video 3 shows the flow structure evolution of the stream function field during the same time history as in Figure 10.