Numerical Study of the Routes toward Chaos of Natural Convection within an Inclined Enclosure

Laboratoire d’Energies Thermiques Renouvelables, Dé partement de physique, UFR Sciences Exactes et Appliquées, Université de Ouagadougou, 03 BP 702 1 Ouagadougou 03, Burkina Faso Laboratoire de Mahtématiques Appliquées, Départemen t d physique, UFR Sciences Exactes et Appliquées, Université de Ouagadougou, 03 BP 7021 O uagadougou 03, Burkina Faso Laboratoire de Mahtématiques et Physiques des Systè me (MEPS), Université de Perpignan, 52, avenue Paul Alduy 66860 Perpignan Cedex, France


Introduction
Natural convection heat transfer under different physical configurations has been studied extensively numerically and experimentally due to its wide applications such as thermal storage plants, thermal buildings [1][2], industrial processes such cooling of electronic fitting [3], stochastic climate dynamics [4].The spatio-temporal behaviour of fluid motion in cavities with heated and cooling vertical walls and horizontals walls adiabatic is qualitatively different depending on aspect ratio values and the Rayleigh number [5][6].Generally, when Rayleigh number is small, the flow is laminar and becomes more complicated [7] as Rayleigh number increases.In Rayleigh-Bénard convection, the flow shows, when the aspect ratio is small, a variety of routes to chaos such as quasi-periodic and intermittencies, depending on the Prandtl number [8].Spatio-temporal intermittencies appear with large aspect ratio [9][10][11].Yang [12] showed that more new and refined experimental and numerical techniques are needed to get better and more accurate resolutions of the spatial and temporal flow details.Four routes to chaos are frequently observed in low-dimensional chaotic systems: the quasi-periodic route via a sequence of Hopf bifurcations , the period-doubling cascade route [13][14], the direct transition to chaos via intermittency [15], and the homoclinic catastrophe.A numerical study of the bifurcation structure and stability of multiple solutions for laminar mixed convection in a rotating curved duct of square cross-section showed that flow and temperature fields on various solution branches are symmetric/asymmetric multi-cells patterns [16].Furthermore, stable steady 2-D solution, stable multi-cell solution, periodic oscillation, chaotic oscillation and symmetric -breaking oscillation led by sub harmonic bifurcation are displayed and the sub harmonic bifurcation is identified to be another route to the chaos.
A numerical simulation [17] of the two-dimensional thermal convection in a rectangular cavity showed that the pitchfork bifurcation is structurally unstable to the inclination angle of the cavity if the bifurcation flow is anti-symmetric, and is structurally stable if it is symmetric in the horizontal direction.The study of convection heat transfer in a square enclosure, with 45° angle inclination [18][19], showed that the larger the Rayleigh number is, the more sensitive the attractor becomes to time steps and mesh grids.The system transits from a fixed point toward chaos via a limit cycle, a period doubling cascade, periodic windows, and tangential bifurcation.
In this research, natural convection heat transfer in a cavity inclined with respect to the horizontal plane is studied numerically.The cavity is filled with air.The effect of the Raleigh number and the inclination angle of the cavity on flow and heat transfer and the route toward chaos are investigated for with aspect ratio equal to1.

Formulation of the Problem
A schematic representation of the system under study is shown in Fig. 1.We consider the unsteady natural convection heat transfer in a cubic cavity filled with air inclined with respect to the horizontal plane (Fig. 1).Initially, the fluid (air) in the cavity is at rest at a temperature T c and at time t 0 , the vertical sidewalls (AD and CD) are heated to T c while the sidewalls BC and AD are maintained at temperature T h ; so that T h >T c .It is assumed that flow and heat transfer are two-dimensional.The flow is described by Navier-Stokes and energy equations with the Boussinesq assumption [20].Navier-Stokes equation, transformed using the stream function formulation and the energy equation are written in non dimensional form in a Cartesian coordinate system: -Initial conditions Natural convection in the cavity is produced by temperature difference between hot and cold vertical walls.Initially, the cavity walls are supposed to be of the same temperature.Therefore, we used standard conditions where t 0 is the starting (initial) time.
-Boundary conditions -For 0 < y < 1; x = 0 and x = 1 In practical applications, a quantity of physical interest is to be determined, which is heat transfer at the walls.This may be obtained in term of local Nusselt number, Nu AB (x), Nu DC (x) respectively for heated sides walls and Nu AD (y), Nu BC (y) respectively for cooled sides walls, from the following relations: Governing Eqs.(1-6) are solved using the Alternating Direction Implicit Method (ADI) in a uniform grid and the Gauss elimination method [18].The convective terms are discretised by a centred scheme and the wall vorticity boundary condition is determined by Woods correlation [19].Eq. ( 3) is solved using a successive over relaxation method [18,21].For each time steps, convergence was declared when the following criterion is met.

Validation
In order to validate our numerical method, we propose to compare ours results with the most closely related numerical solution [18][19][20][21][22][23][24][25][26].We notice that for a 45° inclined square enclosure, with respect to the horizontal plane and heated from two opposite sides, ours results about the Nusselt numbers and the maximum stream function value align with the results by Skouta et al. [18] (Table 1).To ensure that the results from the numerical study are independent from the computational grid, a grid sensitivity analysis was carried out for several Rayleigh numbers.Computations are carried out with several mesh sizes, several step time and Rayleigh number values [27].

Results and Discussion
Results are in various forms, including T 2 torus, strange attractor, Poincaré section and Fourier frequency spectrum for two inclination angles: 25° and 65° and Rayleigh number varying between 10 6 and 7.5x10 6 .

Bifurcation limits point/limit cycle
For Rayleigh numbers inferior or equal to 1.95x10 6 , dynamic variables evolve as time increases toward an asymptotic limit of a steady state.For example, Fig. 2 exhibits for Ra=1.95x10 6 the evolution during the time of the Tm (Fig. 2a) and the trajectory in the plane phase (Tm, Psimax) (Fig. 2b).We notice that the system evolves toward an asymptotic limit of a steady state for both inclination angle (25° and 65°) used in this study.Results on the evolution at fixed point of the Nusselt average number as the Rayleigh number increases can be presented for the inclination angle of 25° by:  As the Rayleigh number increases, the flow intensity and the computing time to obtain the attractor increase.We report on Fig. 6 the evolution during time of Tm (Fig. 3a) and the trajectory in plane phase (Tm, Psim) (Fig. 3b).We get the value of this number for which the attractor becomes periodic as the Rayleigh number increases.This critical Rayleigh number is between 1.95x10 6 and 1.96x10 6 (Fig. 3).The plotting of the amplitude spectrum of dynamic variables obtained by Fast Fourier Transform corroborates the existence of a limit cycle (Fig. 3c).We notice that with Ra =1.958x10 6 around the bifurcation point, the frequency is constant (Fig. 4b) whereas the amplitudes increase linearly with (Ra-1.958x10 6) 1/2 (Fig. 4b), meaning that the bifurcation from the point attractor to the limit cycle is an overcritical Hopf bifurcation [21]

T 2 torus
Fig. 8 shows two frequencies whose irrational ratio Ra ≈ 3.5x10 6 .This implies a second value of critical Rayleigh number around 3.5x10 6 , for which the periodic regime undergoes a transition by destruction of the linearly stability toward a periodic cycle and T 2 torus.We notice for the two inclination angles, the stability of the quasi-periodicity of the two frequencies for Rayleigh numbers, varying between 3.5x10 6 and 3.64x10 6 .For example, Fig. 8 gives Ra=3.613x10 6 , the temporal signal of Nu h (Fig. 5a), the spectrum amplitude of Nu h (Fig. 5b), the trajectory in the plane phase (Fig. 5c) and the Poincaré section (Psimin, Psimax; Tm = 0.54) (Fig. 5d) which includes two heaps of points.The power spectra of Nuh, Tm, Psim and Psimax shows the presence of two new rays and the harmonics of which are a linearly combination put as m 1 f1+m 2 f2 in which m 1 and m 2 are natural integers.

Cycle fitted on a T 2 torus
We report on Fig. 6 the evolution, for Ra= 3.7x10 6 , during Nu h time (Fig. 6a), its power spectrum (Fig. 6b), the trajectory in the plane phase (Nu h , Tm) (Fig. 6c) and the Poincaré section (Psimin, Psimax; Tm = 0.56) (Fig. 6d).We notice that ratio fe/fb is rational for Rayleigh numbers near 3.65x10 6 and Fourier Spectrum rays are harmonics of the lower frequency.So, both frequencies hook and the flow periodic are characterised by a T 2 torus.

Destabilization of the system by appearance of a chaotic window
Increasing Ra from 4.10x10 6 to 4.68x10 6 leads to the onset of peaks for which the power spectrum increases to become dense and in chaotic state (Fig. 7).The chaotic state seems to be confirmed by the Poincaré section and by the positive Lyapanov exponent (Table 6).This result aligns with those by Skouta [18].Lets us now examine the impact of the increase of Rayleigh number on the natural convection flow.We notice that with Ra=4.7x10 6 and 25° inclination angle, both frequencies are as follow: f1=47 and f2=53, for which the ratio is irrational (Fig. 8).It means the apparition of T 2 torus.The T 2 torus is observed until Ra =5.6x10 6 for 25° inclination angle.The increasing in the inclination angle leads to an increase in the range of values of the Rayleigh number for which T 2 torus exists.We notice that for an inclination angle of 65°, T 2 Torus appears until Ra = 5.00x10 6 .

Strange attractor
For the both inclination angles, the attractor is, for Ra =5.7x10 6 or Ra =5.5x10 6 , a cycle fitted on T 2 torus.We notice that for an inclination angle of 25° and Ra = 5.9x10 6 , Ra= 6x10 6 the Tm temporal signal, its power spectrum and its great sensitive dependence on initial conditions as the trajectory in the plane phase (Tm, Psim) are characteristic of a chaotic attractor as illustrated in Fig. 9.However, the Tm temporal signal shows irregular amplitudes characterising the chaos state (Fig. 9a).The Tm power spectrum presents a noise (Fig. 9b).So, it is very difficult to determine the frequencies.The Poincaré section is formed of scattered points (Fig. 9d).Results reported above show that the attractor is strange.As Rayleigh number continues to increase, the chaotic regime expands.
We calculated the largest Lyapunov exponent values versus Ra using the relation below: where T 01 = 0 and T 02 = 10 -6 are initial temperatures.
The great positive Lyapanov exponent increasing with Rayleigh numbers (Table 2) shows the intensification of the chaos and that the attractor is strange.Fractal dimensions of attractors have been calculated using Saraille code FD3 [28-29] which uses the iterations of physical variables for example Nuch, Psimax, Tm and their successive derivatives with regard to time.This sets sizes, defines a pseudo-space of the phases of dimension n + 1, and where n is the number of derivatives.When n increases, FD3 gives values which can increase until the limit corresponding to the dimension d of the attractor with nearly 10% incertitude margin.Calculations give d > 2 for a number of Rayleigh above 6.10 6 .The use of FD3 confirms the appearance of a strange attractor when d is above 2.

Conclusion
Natural convection in an enclosure inclined to the horizontal plane was numerically investigated in this study in order to determine the routes to chaos.This enclosure is heated from two opposite side and cooled on the other two sides.Transfers equations are solved using the vorticity-stream function formulation, the alternating direction implicit method (ADI) and the Gauss elimination method.The impact of Rayleigh number and the inclination angle are examined in detail.Major results are: -The routes to the chaos for complementary angles (25° and 65°) are identical.
-With 25° inclination angle, the attractor is a fixed point for Rayleigh number less than 1.95x10 6 ; it transits from fixed point toward limit cycle via an overcritical Hopf bifurcation for a value of Rayleigh number between ranging 1.95x10 6 and 1.96x10 6 .-For the Rayleigh number varying between 5.00x10 6 and 3.64x10 6 , the limit cycle undergoes a second Hopf bifurcation which turns it into T 2 torus.A cycle aligned on T 2 torus appears after a hooking of frequencies which continues until Ra = 4.00x 10 6 .-A chaotic sequence seems to appear for Ra varying between 4.10x0 6 to 4.68x10 6 and a laminar flow appears from Ra =4.7x10 6 .-For Ra higher than 5.6x10 6 , a hooking of frequencies brings back the T 2 torus to a limit cycle fitted on it.-When Rayleigh number increases up to 5.9x10 6 , the attractor is chaotic again.

Fig. 1 .
Fig. 1.Physical layout of the system.A with:[x = x * /H], [y = y * /H], ω * .H 2 /a, ψ = ψ * /a and t = t*a/H 2where T represents the dimensionless temperature, ω, the vorticity and ψ the stream function.Ra is Rayleigh number and Pr Prandtl number.The associated initial and boundary conditions are: on the cooled wall BC The average Nusselt numbers AB Nu , DC Nu on the hot walls AB, DC respectively and AD Nu , BC Nu on the cooled walls AD, BC respectively are obtained by integration of the local Nusselt number through the corresponding surface and are given by ∫

Fig. 2 .
Fig. 2. Representation of the limit point attractor before limit cycle for Ra = 1,95.10 6.a): Temporal signal of the temperature Tm, b): Trajectory in the phase plane (Tm, Psimax).
NuH global Nusselt number on the hot sides AB and CD [ ∫ NuC global Nusselt number on the cooled sides AD and BC [ ∫ Nx node number in the [Ax) direction Ny node number in the [Ay) direction Pr Prandtl number: along x direction : u * .a/Hv dimensionless velocity component along y direction :v * .a/Hx, y velocity coordinates: x * /H, y * /H Greek symbols α inclination angle of the side wall AB with the horizontal plane (rad).β thermal expansion coefficient (K -1 ) λ thermal conductivity (W.m -1

Table 2 .
The largest Lyapunov exponent values for different Rayleigh numbers when tilt angle is 65°.