Abstract: In this study, we investigated numerically the transient natural convection in a square cavity with two horizontal adiabatic sides and vertical walls composed of two regions of same size maintained at different temperatures. The flow has been assumed to be laminar and bi-dimensional. The governing equations written in dimensionless form and expressed in terms of stream function and vorticity, have been solved using the Alternating Direction Implicit (ADI) method and the GAUSS elimination method. Calculations were performed for air (Pr = 0.71), with a Rayleigh number varying from 2.5x105 to 3.7x106. We analysed the effect of the Rayleigh number on the route to the chaos of the system. The first transition has been found from steady-state to oscillatory flow and the second is a subharmonic bifurcation as the Rayleigh number is increased further. For sufficiently small Rayleigh numbers, present results show that the flow is characterized by four cells with horizontal and vertical symmetric axes. The attractor bifurcates from a stable fixed point to a limit cycle for a Rayleigh number varying from 2.5x105 to 2.51x105. A limit cycle settles from Ra = 3x105 and persists until Ra = 5x105. At a Rayleigh number of 2.5x105 the temporal evolution of the Nusselt number Nu(t) was stationary. As the Rayleigh number increases, the flow becomes unstable and bifurcates to a time periodic solution at a critical Rayleigh number between 2.5x105 and 2.51x105. After the first HOPF bifurcation at Ra = 2.51x105, the oscillatory flow undergoes several bifurcations and ultimately evolves into a chaotic flow.
INTRODUCTION
Natural convection is a recurrent phenomenon in the world around us and most of these natural convection flows are turbulent, especially those encountered in engineering applications. Unsteady and turbulent natural convections have thus attracted increasing interest over the last decades for twomain reasons: first, the desire to improve our phenomenological understanding of turbulent natural convection and second, the pressing need for numerical models capable of predicting the corresponding flow structures and related heat transfer in industrial applications. Natural convection flows in enclosures are usually subdivided into two main classes, those heated from below and those heated from the side. The configuration of the differentially heated cavity models in many engineering applications such as cooling of electronic components, nuclear reactor insulation, ventilation of rooms etc., involves building system containing multilayered walls, double windows and other air gaps in unventilated spaces, energy systems like solar collectors, energy storage devices, furnaces, heat exchangers, materials processing such as solidification phenomena and growing crystals. A literature survey over the last decades shows that the natural convection transfer equations are solved using the stream function and vorticity formulation and the implicit numerical method proposed by Douamna et al. (1999). The simple generalised differentially quadrature method has been used by Shu and Week (2002). In the past years, most studies relevant to natural convection in cavities have focused on laminar steady and transient flow regimes (Devahldavis, 1983), transition to unsteadiness and route to chaos (LeQuere and Behnia, 1998; Ndame, 1992). The conclusions are as follows: a laminar steady natural convection was observed for small values of the Rayleigh number. When the Rayleigh number is increasing above a critical value, unsteady convection is observed due to the loss of stability of steady solutions. Victorri and Blondeaux (1991) numerically studied the transition process that leads to the oscillatory flow over a wavy wall from a periodic behaviour to chaos using a FEIGENBAUM scenario. They found that by increasing the Reynolds number, the flow experiences an infinite sequence of period doubling (Pitchfork bifurcation) which takes place at successive critical values. Yoo and Han (2000) have studied numerically the bifurcation sequences to the chaos in a natural convection in a horizontal annulus cavity. For a Prandtl number Pr equals to 0.1, the HOPF bifurcation appears for a critical Rayleigh number equals to 1700. Skouta et al. (2001) showed that in an inclined cavity filled with air, the critical Rayleigh number is equal to 1.11x105 for the first bifurcation. Gollub and Besson (1980) identified experimentally four different routes to turbulence in Rayleigh-Benard convection: Transition to chaos via quasi-periodic state with two or three incommensurable frequencies, period-doubling and intermittency. The bifurcation sequence to temporal chaos (Schuster, 1984) has been an interesting subject in many areas of nonlinear systems and has also been studied widely in fluid dynamics. Many routes to chaos have been found theoretically and experimentally for free and forced convective flows. Lamsaadi et al. (2006) have studied numerically and analytically the flow and heat transfer characteristics and multiplicity of steady states for natural convection in a horizontal rectangular cavity filled with Non Newtonian power-law fluids and heated from all sides. This problem has been implemented in some industrial thermal processes.
Mukutmoni and Yang (1993) found two transitions which were documented numerically. The first transition was from steady-state to oscillatory flow and the second was a subharmonic bifurcation as the Rayleigh number was increasing further. This study has been done in a rectangular enclosure with insulated side walls. The aspect ratios were 3.5 and 2.1 and the boussinesq fluid was water with a Prandtl number of 2.5. Bratsun et al. (2003) have studied both experimentally and numerically the convective flow in a tall vertical slot with differently heated walls. The flow was investigated for the fluid with Prandtl number Pr = 26. They noted that it started with a plane parallel flow as primary solution, thus became unstable to two counters propagating waves. As the Grashof number was increasing even further, a chaotically oscillating cellular pattern consisting of the pieces of broken waves arised. The formation of a structure in the form of vertical rolls chaotically modulated along axes concluded this complicated picture. Kuang et al. (2006) have developed an algorithm for suppressing the chaotic oscillations in non linear dynamical systems. This algorithm used the Lyapunov-Krasovskii (LK) method. The main attention of the study of Yu and Leung (2003) was focused on the computation of the SNF (the simplest normal form) of HOPF bifurcation with perturbations parameters. Ishida et al. (2006) have proposed an accurate and simple method to evaluate the Lyapunov spectrum. This method is suitable for any discretization method that finally expresses a governing equation system in the form of an ordinary differential equation system. The method was applied to evaluate up to the second largest Lyapunov exponents for natural convection in a rectangular cavity with heated and cooled side walls.
Prasad and Das (2007) have presented a study of mixed convection inside a rectangular cavity. The integral form of the governing equations were solved numerically using finite-volume method. SIMPLE algorithm with higher-order upwinding scheme was used. The Grashof number has been given the values of 0, 104 and 106, while the aspect ratio (height/width) was fixed to 0.5, 1 and 2, keeping the Prandtl number (Pr) = 1. A HOPF bifurcation has been observed at Gr = 105.
The problem of natural convection in a cavity with adiabatic horizontals walls and verticals walls composed of two regions of same size maintained at different temperatures and with upper half temperatures less than the lower half temperature, has been studied, at our knowledge, only by Jahnke et al. (1998). For symmetric boundary conditions, the structure of the flows is, in some ways, similar to the Rayleigh-Benard problem in Hele-Shaw cells. Our main attention for this study was to describe the routes to the chaos and to investigate the periodic, the quasi periodic and the chaotic regimes. In the present case, we considered that the cavity is filled with air which a Prandtl number equals to 0.71 and an aspect of ratio of 1. Transfer equations were solved by an implicit finite difference method. Forward differences were used for the time derivatives and central differences for space derivatives. The velocity components were calculated with a central finite difference approximation of the stream function. The effects of Rayleigh numbers on the natural convection and the roads of the chaos that the system borrowed have been studied. We plotted the temporal evolution of the hot global Nusselt number, we represented the amplitude spectrum with the Fast Fourier Transform (FFT) and the attractors in a space trajectory.
We calculated the fractal dimension versus the space trajectory dimension. Chaos was determined by calculating the largest Lyapunov exponent. This exponent was estimated using a one dimensional time series. In this case, the average temperature was used as the time series.
PROBLEM FORMULATION
We considered an enclosure with a ratio
were H is the height and L the width. The horizontals walls are adiabatic and the verticals walls are composed of
| Fig. 1: | Physic model |
two regions of same size maintained at different temperatures. The upper half temperatures are less than the lower half temperature. Figure 1 shows the schema of the physical model along with the system of coordinates. The bottom wall coincides with the origin of a Cartesian coordinate system in which the x-axis is the horizontal distance from the left vertical wall and the y-axis is the vertical distance measured from the bottom wall that is in the opposite direction to the gravitational force.
In this enclosure filled with air, transfers are governed by the natural convection equations. The following assumptions are made:
| • | The thermo physical properties of the fluid are constant except for the density variation that induces the buoyancy forces, which follow the Boussinesq approximation. |
| • | The viscous dissipation in the energy equation is neglected. The flow is laminar, bi-dimensional and incompressible. Radiation exchanges between walls are neglected. Air is supposed to be an ideal gas |
Under these assumptions, the dimensionless equations using the stream function-vorticity formulation are as follows:
| (1) |
| (2) |
| (3) |
with
Where, Ω is the dimensionless vorticity, Ψ the dimensionless stream function, u the dimensionless velocity component along the x-direction and v the dimensionless velocity component along the y-direction. The boundary and initial conditions consist of:
Initial conditions
(a) t<=t0 t0 Corresponds to the time when the vertical walls are subjected to the temperatures Th (hot) and Tc (cold).
0<=x<=1, 0<=y<=1: T = 0 Ω = 0 Ψ = 0 u =
0 v = 0 |
Boundary conditions
The average Nusselt numbers are defined as follows.
NUMERICAL PROCEDURE
The energy, the momentum and the elliptic Poisson equations for the stream function associated to the boundary conditions are solved using the Alternating Direction Implicit method [ADI] and the Gauss algorithm method. The velocity field is calculated by central differences scheme. For the iterative method, the convergence criterion for Ψ is defined as:
while the convergence criterion for T or Ω is defined as
Where, f = T, Ω, . FK and fk+1 are the values of T or Ω at iterations k and (k+1), respectively.
FRACTAL DIMENSION
The signal sampling allows the reconstitution of not only the trajectory but
also a set of N points on this trajectory. These N points reveal the fractal
structure of the attractor based on their positions symbolized by the coordinate
vector
Where, H(x) represents the Heaviside function, which is equal to 0 for x < 0 and 1 for x > 0. The symbol ||…|| represents the Euclidian norm vector, i and j represent the point indexes and r is the radius of the hyper-sphere containing the N points.
For the N selected points, the function becomes
We sum for the index j, i.e., we investigate the j points relatively to the
C(r) is proportional to the average number of points from the set points that are located inside the hyper sphere centered in some i points of the set:
For a small value of the radius, this number is proportional to the dimension N(r). For an attractor, N(r) ~rd, where d is the fractal dimension; it is proportional to
For a small value of r, d is obtained from the slope of the characteristics maps log (N (r)) = f (log (r)).
Lyapunov exponent: A chaotic solution X (t) is represented by a strange attractor. The sensitivity of the initials conditions implicates a positive Lyapunov Exponent (Ly). This result shows that the next curves are divergent. For a periodic and quasi periodic solution, the Lyapunov exponent is negative. All the chaotic phenomenon are very sensitive to initials conditions.
The discrepancy Δx between two values Φ1 and Φ2 which are initially very close evolved with the time follows an exponent law:
| Ly | = | Largest Lyapunov exponent |
| t | = | Non dimensional time |
Φ1 and Φ2 are the same parameter calculated from the different initial condition but very imminent.
RESULTS AND DISCUSSION
Ours calculations were performed for the mesh 121*121 with a step time of 10-5. We represented the stream function lines in Fig. 2a-c and the isothermal lines in Fig. 3a-c for a Rayleigh number of 103, 104 and 105. The flow is quadricellular for a Rayleigh number equals to 103, the cells are symmetric over the vertical axis. The flow is anything in the vertical and horizontal median axes. The cell in the left inferior zone and those localised in the right superior zone turn in the trigonometric sense, the cell turns in the clockwise sense. Figure 3a shows that the heat transfer is purely conductive.
In Fig. 2b and c, there is a big cell at the centre of the cavity which turns in the trigonometric sense and drives the heat correctly through the cold wall. In Fig. 3b, the isotherm lines represent a thermal stratification while Fig. 3c shows that the isotherm lines walk along the wall of the cavity. A boundary layer appears and the isotherms are then crowded.
By increasing the Rayleigh number, the attractor reaches a stable state. In Fig. 4, the phase trajectory is a limit point at Ra = 2.5x105. This limit point losses its stability and allows the periodic solution to appear at Ra = 3x105 according to the theory of Floquet (Berge and Pomeau, 1998).
The solution is a stable attractor which corresponds to a limit cycle in Fig. 5. With a Rayleigh number of 5.105 there is doubling period apparition which corresponds to a sub harmonic bifurcation. At a Rayleigh number equals to 7x105 the space trajectory is a tore (Fig. 6) and a new bifurcation appears; thus, the system undergoes some instabilities which lead to the chaos. This phenomenon has been showed at the Rayleigh number of 106 as seen in Fig. 7 and the whiff appears at the Rayleigh number of 1.5x106 and this state is intermittent Fig. 9. The dynamic system evolve and attractor become strange as in Fig. 8 for a Ra = 3x106 because the attractor dimension is fractal. The obtained bifurcations follow the Ruelle and Takens scenario (Berge and Pomeau, 1998). In Table 1, we computed the relative errors between the theoretical attractor dimension and the values obtained from our code. The validation is thus acceptable since the discrepancies between the values are from 0.83 to 2.32%.
| Fig. 2: | Stream function lines |
| Fig. 3: | Isotherms lines |
| Fig. 4: | Limit point, (a) temporal evolution and (b) phase space |
| Fig. 5: |
Limit cycle, (a) temporal evolution, (b) phase space and (c) amplitude spectrum |
In Table 2, we show that the fractal dimension is varied with the Rayleigh number and the space trajectories dimension. We plotted the characteristics curves for the theoretical attractors and for the attractors of our system in Fig. 10 and 11, respectively. These curves confirm our results.
| Fig. 6: | Tore, (a) temporal evolution, (b) phase space and (c) amplitude spectrum |
| Fig. 7: | Chaos beginning, (a) temporal evolution, (b) phase space and (c) amplitude spectrum |
| Fig. 8: | Chaos, (a) temporal evalution of the hot global nusselt No., (b) phase space and (c) amplitude spectrum |
| Fig. 9: | Laminar windows, (a) temporal evolution, (b) phase space and (c) amplitude spectrum |
| Fig. 10: | Characteristic lines in (N(r)) = f (in (r)) |
| Table 1: | Theoretical attractor dimension |
| Table 2: | Attractor dimension |
Initial conditions sensibility: Figure 12 shows, for a Rayleigh number equals to 106, the temporal evolution of the average temperature calculated from different initial values; these values are very close (error of 10-3). We notice that the curves are first almost overlapping, then they separate and evolve differently.
Divergence grade of evolution: In Fig. 13, we show the discrepancy between Tave 1 and Tave 2. The curve is plotted with a logarithmic scale. We note that the discrepancy increases with time according with the exponent law. An estimation of the largest Lyapunov exponent is done by the average slope of the curve. The obtained value, 78.94, confirms the sensitivity of the initials conditions.
| Fig. 11: | Characteristic curves for different Rayleigh number |
| Fig. 12: | Temporal evolution of the average temperature |
| Fig. 13: | Divergence grade of evolution |
CONCLUSION
In conclusion, the unsteady natural convection study has showed a transition to the determinist chaos. We examined the dynamic phenomenon in a filled air enclosure. For the determination of the system state, we used the Fast Fourier Transform, the largest Lyapunov exponent and the fractal dimension attractor. We plotted the space trajectory for different Rayleigh numbers and then dimensioned the attractors. We noticed that the temporal evolution of the total Nusselt number is asymptotic for a Rayleigh number equals to 2.5x105.
The range of the first bifurcation has been localised between a Rayleigh number of 2.5x105 and 2.51x105. The system transits from the limit point to the limit cycle according to a HOPF bifurcation. At a Rayleigh number of 3x105, the flow is periodic. When the Rayleigh number increases, the attractor undergoes some period doubling and for a Rayleigh number equals to 5x105, this phenomenon appears in the temporal evolution. At a Rayleigh number of 7x105, the attractor becomes a tore in the space trajectory. The tore undergoes an instability and becomes chaotic at a Rayleigh number of 106. We obtained the laminar window for a Rayleigh number equals to 1.5x106. We noticed that the chaotic whiff are longest and approach.
Then the chaos installs definitively at a Rayleigh number equals to 3x106. The system states have been quantified from the space trajectories. We have elaborated a code based on the Grassberger and Procaccia formulation for the attractor dimension. Our code has been validated with acceptable relative errors. The attractor dimension values are conform to the space trajectory. The largest Lyapunov exponent has been calculated and concluded that the attractor is Strange.
NOMENCLATURE
| a | = | Thermal diffusivity (m2 sec-1) |
| A | = | Aspect ratio (height/length) |
| g | = | Acceleration of the gravity (m2 sec-1) |
| Gr | = | Grashof number |
| H | = | Cavity height (m) |
| L | = | Cavity length (m) |
| Ly | = | Lyapunov exponent |
| Nu | = | Nusselt No. |
| Pr | = | Prandtl No. = v/a |
| Ra | = | Rayleigh No. = g βH3 Δ T/γ.a |
| Racr | = | Critical Rayleigh number |
| t | = | Dimensionless time |
| t` | = | Time (s) |
| T | = | Dimensionless temperature |
| Tc | = | Cold wall temperature (°C) |
| Th | = | Hot wall temperature (°C) |
| Tave 1 | = | Average temperature for the first initial condition |
| Tave 2 | = | Average temperature for the second initial condition |
| ΔT | = | Temperature difference = Th-Tc (°C) |
| u | = | Dimensionless velocity component along x-direction |
| v | = | Dimensionless velocity component along y-direction |
| u` | = | Velocity component in the x-direction (m sec-1) |
| v` | = | Velocity component in the x-direction (m sec-1) |
| x | = | Dimensionless horizontal coordinate axis |
| y | = | Dimensionless vertical coordinate axis |
| x` | = | Horizontal coordinate axis (m) |
| y` | = | Vertical coordinate axis (m) |
| f | = | Signal frequency (sec-1) |
| Ψ | = | Dimensionless stream function |
| Ψ` | = | Stream function (m2 sec-1) |
| Ω | = | Dimensionless vorticity |
| Ω` | = | Vorticity (sec-1) |
| ς | = | Density (kg m-3) |
| γ | = | Kinematic viscosity (m2 sec-1) |
| β | = | Coefficient of thermal expansion (1/K) |