Onset of Natural Convection and Transition Laminar-Oscillatory Convection Flow in Rayleigh-Bénard Configuration

Onset of Natural Convection and Transition Laminar-Oscillatory Convection Flow in Rayleigh-Bénard Configuration

Said BouabdallahBadia Ghernaout Mohamed Teggar Ahmed Benchatti F-Zohra Benarab 

Laboratory of mechanics, University Amar Telidji, Laghouat 03000, Algeria

Corresponding Author Email: 
fibonsaid@gmail.com
Page: 
151-157
|
DOI: 
https://doi.org/10.18280/ijht.340122
Received: 
| |
Accepted: 
| | Citation

OPEN ACCESS

Abstract: 

A numerical study of Rayleigh-Bénard convection in a rectangular cavity has been presented. The onset of natural convection and the transition from laminar to oscillatory convection were considered in this study. The finite volume method was used to solve numerically the governing equations of the phenomenon and the pressure-velocity coupling was matched by SIMPLER Algorithm of Patankar. This study was carried out for Rayleigh number varied from 103 to 106 in order to control the first value of critical Rayleigh number, Rac-1 corresponding the onset of convection, and that for different aspect ratios of the cavity. The threshold transition regime of laminar-oscillatory convection is depicted by the second critical Rayleigh number Rac-2. Also, the different modes of bifurcation of convection were determined and discussed.

Keywords: 

Rayleigh-Bénard convection, Natural convection, Oscillatory flow, Bifurcation.

1. Introduction

Due to its practical importance in many general science and engineering applications, Rayleigh–Bénard convection has been the subject of several theoretical, experimental, and numerical studies. Since, Rayleigh–Bénard convection presents the evolution from the stationary state to the fully developed turbulent regime with many different flow patterns and sequences of bifurcations; it is widely investigated as the problems of different transition mechanisms in hydrodynamics [1-2]. Most of the published works covering natural convection in enclosures that exist today can be classified into two categories: differentially heated enclosures [3-6] and enclosures heated from below and cooled from above (Rayleigh Bénard problems) [7]. Benchmark solutions related to differentially heated enclosures (first group) can be found in several numerical investigations [8]. However, numerical benchmark solutions related to the simplest case of 2D Rayleigh-Bénard convection are less encountered in the literature.

Some recent development in turbulent Rayleigh-Bénard convection was carried out by Lappa [9]. The case of two dimensional square enclosures was considered, the author began by the bifurcation and the symmetry breaking system. Four cases of possible symmetries were found (one cell, two horizontal cells, two vertical cells, and four cells). At the beginning of bifurcation, the one cell case was got. For the time dependence, other forms can appear following the oscillatory mode. Venturi et al. [10] studied the stochastic bifurcation and stability of the natural convection of Rayleigh-Bénard by different stochastic modelling approaches. They focused on fluid in the supercritical state and depicted the critical Rayleigh number at 2585 in a square cavity. They focused also on the sensibility of the initial conditions and bifurcation in steady state.

To controlling the amplitude of bifurcated Mas et al. [11] studied the bifurcation and stability of the Boussinesq equations and the onset of the Rayleigh-Bénard convection. A nonlinear theory for this problem was established in this paper using a new concept of bifurcation called attractor bifurcation. Kao et al. [12] investigated the 2D Rayleigh-Bénard convection using a simple LB model with the Boussinesq approximation, from the threshold of the primary instability with a theoretical value of critical Rayleigh number Rac = 1707:76 to the regime near the flow bifurcation to the secondary instability. Their results show that the periodic unsteady flows take place at certain Prandtl numbers with an appropriate Rayleigh number. Furthermore, the Nusselt number is found to be relatively insensitive to the Prandtl number in the ranges 0.71 < Pr < 70 and Ra < 105. Finally, the relationship between the Nusselt number and the Rayleigh number is also investigated.  Angeli et al. [13] considered a system of an air-filled square-sectioned 2D enclosure containing a horizontal heated cylinder. The results were shown in respect to the variation of Rayleigh number and the aspect ratio. The first bifurcation of the low-Ra fixed-point solution is tracked for each A-value, and chaotic flow features are detailed for the case A = 2.5. The analysis of global heat transfer data showed that Nu and Ra relationship is sensibly influenced by the transition mode associated with each value. A correlating equation for the average Nusselt number on the cylinder, derived for the subcritical case, was found to be valid up to slightly supercritical Ra value. Raji et al. [14] presented a numerical results of the natural convection in a square cavity filled with air, the temperature of the lower horizontal surface is kept constant(hot), while that of the upper surface is maintained at a cold temperature, the remaining upright walls are considered adiabatic. The Rayleigh number was considered in the range 103 < Ra < 7×106. Three different solutions were obtained (single-cell flow, two-cell vertical flow and horizontal flow bi-cellular).

Recently, Ghernaout et al. [15] presented the double diffusive natural convection in binary mixture under the effect of external magnetic field for steady and oscillatory state. The obtained results show a strong dependence of the structure of thermal and solutal effect to the buoyancy ratio and the oscillatory double diffusive flow occurs form a periodic time evolution where the phenomenon change around in each period time. A critical thermal Rayleigh number RaTCr and corresponding dominated frequency for the onset of oscillatory double diffusive convection were determined according the buoyancy ratio for different values of Hartman number (0, 25, 50, and 100).

The aim of this study is to propose 2D numerical solutions related to the natural convection in enclosure. We have interested to the determination of the onset of natural convection and the transition from laminar to oscillatory convection in Rayleigh–Bénard configuration. Also, the different modes of bifurcation of convection were determined and discussed.

2. Mathematical Formulation

The cavity which is heated from below and cooled from above corresponds to the configuration of the Rayleigh-Bénard dealing with the stability and motion of a fluid confined between two horizontal plates that are maintained at uniform temperatures (Fig. 1). We consider that the flow is incompressible and satisfied the Boussinseq approximation. To give the conservation equation in dimensionless form, we have used the dimensionless variables respectively:$X=\frac{x}{H}, Y=\frac{y}{H}, U=\frac{u}{v / H}, V=\frac{v}{v / H}$$P=\frac{p}{\rho(\alpha / H)^{2}}, \quad \theta=\frac{T-T_{f}}{T_{C}-T_{f}}, \tau=\frac{t}{H^{2} / \alpha}$ for cartesian coordinates, velocity components, the pressure, temperature and the time respectively. The dimensionless equation are given as follows:

  • Mass conservation equation

$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$(1)

  • X-momentum equation

$\frac{\partial U}{\partial \tau}+U \frac{\partial U}{\partial X}+V \frac{\partial U}{\partial Y}=-\frac{\partial P}{\partial X}+\left(\frac{\partial^{2} U}{\partial X^{2}}+\frac{\partial^{2} U}{\partial Y^{2}}\right)$(2)

  • Y-momentum equation

$\frac{\partial V}{\partial \tau}+U \frac{\partial V}{\partial X}+V \frac{\partial V}{\partial Y}=-\frac{\partial P}{\partial Y}+\frac{R a}{\operatorname{Pr}} \theta+\left(\frac{\partial^{2} V}{\partial X^{2}}+\frac{\partial^{2} V}{\partial Y^{2}}\right)$(3)

  • Energy equation

$\frac{\partial \theta}{\partial \tau}+U \frac{\partial \theta}{\partial X}+V \frac{\partial \theta}{\partial Y}=\frac{1}{\operatorname{Pr}}\left(\frac{\partial^{2} \theta}{\partial X^{2}}+\frac{\partial^{2} \theta}{\partial Y^{2}}\right)$(4)

The flow parameters appear in Eqs(1-4) are: Rayleigh number and Prandtl number defined respectively as: $R a=\frac{g \beta H^{3} \Delta T}{v \alpha}$ and $\operatorname{Pr}=\frac{v}{\alpha}$.

Figure 1. Schematic of Rayleigh-Bénard probleme.

The initial and boundary conditions are given as follows:

 

At $\tau=0$:  and $\theta=1$(5a)

At $\tau>0$:

$Y=0 : 0 \leq X \leq 1, U=V=0, \quad \theta=1$(5b)

$Y=1 : 0 \leq X \leq 1, U=V=0, \theta=0$(5c)

$X=0,1 : 0 \leq Y \leq 1, U=V=0, \partial \theta / \partial X=0$(5d)

3. Numerical Solution

Equations (1-4) associated with initial and boundary conditions Eqs(5a-d) have been solved by a finite-volume method [16]. The Quick scheme is used for discretization of convective and diffusive terms. The couple’s velocity-pressure is solved by SIMPLER Algorithm. The obtained algebraic equations are solved by the line-by-line tri-diagonal matrix algorithm (TDMA). The convergence is declared when the maximum relative change between two consecutive iteration levels is less than 10-4.

The grid independency of the numerical solution was established on a careful analysis of three grids (60×60), (80×80) and (100×100). We have used a refined grid in the lower and upper wall, due to the existence of a strong temperature and velocity gradients near this walls, refining value is equal to 1.05. We assume for this purpose that the case of a steady natural convection flow with Ra = 104. The Table 1 shows the grids effect on the deferent flow parameter (Ѱ, U, V, and Nu).  In table 2, we have compared our results with those of Ben Cheikh [4]. A Good agreement has been obtained with two results. In addition, we found that our results are in accordance to those availibale in the literature (Venturi [10] and Gelfgat [17]). The latter used the same of grid size. Shishkina [18] determined the minimum number of nodes to have an accurate result. This grid is considered to show the best compromise between computational time and the results precision.

Table 1. Flow Parameter: Ѱmax, Numoy, $\left|V_{\max }\right|$, $\left|U_{\max }\right|$,for different grid size (Ra = 104)

Grid size

Ѱmax

Numoy

$\left|V_{\max }\right|$

$\left|U_{\max }\right|$

60×60

321.79

2.1

26.54

25.7

80×80

316.61

2.02

27.81

26.07

100×100

 316.43

1.98

27.29

26.18

Table 2. Comparison of our results, with those of Ben Cheikh et al. [4] for Ra =104

Our results

Results of [4]

Relative error %

$\left|V_{\max }\right|$

$\left|U_{\max }\right|$

$\left|V_{\max }\right|$

$\left|U_{\max }\right|$

$\left|V_{\max }\right|$

$\left|U_{\max }\right|$

26.54

25.7

26.36

25.22

0.68

1.9

4. Results

All results presented herein are given in the dimensionless form. Comparisons with other studies available in literature have been carried out. A comparison was made with Rayleigh-Bénard convection studied by Turan [19] and the Benchmark results [11] presented in Table 3. We can see a slight relative difference, less than 1.5%, between the two results.

Another comparison was made with the works of Val Davis [8], presented on Table 4. In this case the heat gradient is horizontal. The results compared show a good agreement of the results with a small difference for all flow parameters.

4.1 Onset of natural convection

The heat convection mode is taking place by energy accumulation then the flow motion. The beginning of the convection in the case of Rayleigh-Bénard convection, then we will proceed to the consideration of the phenomenon by studying the threshold of the convection onset in the case of Rayleigh-Bénard convection, and then calculate the critical value of Rayleigh number in order to see the effect of the aspect ratio on this critical value. We will compare it with the theoretical result; and then make a close approach to the infinite cavity. Then and there, we will make the extension of the cavity.

A linear stability analysis of the Boussinesq equations about the linearly conducting profile between two infinite horizontal no-slip plates shows that the critical Rayleigh number Rc-1 ≈ 1707.76. The value of Rc-1 is independent of the Prandtl number [11], this critical value is necessary to determinate the mode of heat transfer. It represents the threshold of the onset of convection and the first bifurcation flow structure. We examine the effect of the aspect ratio on the first Rayleigh critical number (Rac-1). Four aspect ratios of the Rayleigh-Bénard cavity were examined (A = 1, 2, 4 and 8). In Figure 2, we presented the average Nusselt number along the hot wall in function of Rayleigh number. The value of Rac-1 for various aspect ratios was given in Table 5. We note that although the aspect ratio has an influence on the value of Rac-1, therefore the fluid motion is depending to the geometry of the cavity. In the literature, we found that for the infinite cavity case Rac-1 = 1708 but for a square cavity (A = 1) the value of critical Rayleigh is much greater (Rac-1 = 2585.01). In Table 6, we compared our results of Rac-1 for a different aspect ratio with the results of Gelfgat [17]. A good agreement was obtained.

4.2 Onset of oscillatory natural convection

At the beginning of the convection the flow is laminar and does not depend on time. But, an oscillating state appears later and the flow becomes time dependent. We will confirm that the oscillations are physical in origin, and then make an evaluation for different probes into the cavity.

Figure 2. Average Nusselt Number in function of Rayleigh number for aspect ratios A = 1, 2, 4 and 8

According to linear theory, a significant change in convection will take place as soon as the critical Rayleigh number is exceeded and Rayleigh-Bénard problem will become nonlinear. In our study we found that the time dependent flow begins at around 16 Rac-1, it is the oscillatory mode. The critical value of this mode is Rac2= 4.58×104. In order to view the physical oscillatory flow, we have chosen many arbitrary probes in the cavity. We note that the steady state is obtained up to the value of Ra = 4×104 (Figure 3a). However, past this value (Rac2 = 4.58×104), the flow becomes unsteady and not periodic (Figure 3b).

Table 3. Comparison of our results with those of Turan [19] and Benchmark [11]

Ra

Our results

Results of Turan [19]

Benchmark Results [11]

Relative difference with [19] (%)

103

1

1 .0004

1

0.039

104

2.2

2.1581

2.154

1.49

105

3.9

3.9103

3.907

0.26

106

6.4

6.3092

6.36

1.439

 
Table 4. Comparison of our results with those of Val Davis [8]

Ra

$\Psi$

$\Psi$ [8]

Nu

Nu [8]

103

1.171

1.170

1.118

1.117

104

5.061

5.059

2.25

2.240

105

9.147

9.059

4.53

4.505

106

16.22

16.24

8.9

8.810

 
Table 5. Critical Raleigh number for different aspect ratios

Aspect ratio

1

2

4

8

Rac-1

2 585.01

2 013.5

1810

1752

To control the nature of the oscillations (physics and not numerical), a probe was fixed in our cavity with the same flow parameter and we reduced the time step to the half, the figure 4 shows that there is no influence of the time step on the oscillation amplitude. Consequently, these instabilities are physical and not numerical. Another test of oscillatory mode is the phase portrait as showed in figure 5. The phase portraits reflect the change of hydrodynamic and thermal parameters between them. So, for a periodic oscillatory regime at given point in the flow, these changes are closed circles that reflects the periodicity of the flow regime. Moreover, for no periodically oscillations these changes appear become the endpoints or other unordered structure.  This technique has been used successfully by Ghernaout et al. [15].

(a) Ra = 4× 104

(b) Rac2 = 4.58×104

Figure 3. Time-dependent of dimensionless vertical velocity at the probe P1 for two value of Rayleigh number

Oscillatory flow regime occurs for the critical value of Rayleigh number Ra = 4.58×104. We note that steady state is obtained up to the value of Rayleigh number Ra = 4×104. In order to illustrate the temporal evolution of the horizontal, vertical velocity component, and temperature, we chose the same probes (p1, p6, p7, p8) in the cavity to look into the field, and the results are shown in figures 6a-c.

It is clear that these profiles are oscillatory and periodic, so the flow regime is unstable. We observe that the amplitude of these oscillations changes from one point to another in the enclosure. The difference in the degree of oscillation depends on the location of probes in relation to adiabatic walls to the hot wall and the cold wall. For example, probe 1 (Table 6) is in the middle of the enclosure, the probes 6, 7 and 8. Probes 2, 3, 4 and 5 are used to explain the wave propagation flow.

Table 6. Physical location of measurement probes

 

P1

P2

P3

P4

X

0.5

0.75

0.75

0.25

Y

0.5

0.75

0.25

0.25

 

P5

P6

P7

P8

X

0.25

0.2

0.1

0.6

Y

0.75

0.1

0.4

0.9

Figure 4. Time dependent of vertical velocity component for two time setp in P1 (choosing arbitrary) and for Rac2 = 4.58 ×104

Figure 5. Phase portrait of horizontal velocity component in function of temperature at point P2 (choosing arbitrary), for Rac2 = 4.58 ×104

The amplitude of the temperature is greater on the hot wall, as the particles back into the wall; it means average amplitude of the adiabatic wall. But, there is small amplitude for the center point and the cold wall. It is clearly observed on the probes 2, 3, 4 and 5, the same amplitude and the same structure of the wave, except in phase shift that is due to the wave propagation.

(a)Dimensionless horizontal velocity component U

(b)Dimensionless vertical velocity component V

(c)Dimensionless temperature

Figure 6. Time-dependent of dimensionless horizontal (a), vertical velocity component (b) and dimensionless temperature (c), in probes P1, P2, P3, P4, P5, P6, P7 and P8 (choosing arbitrary), for Rac2 = 4,58× 104

To better understand the oscillatory flow, we have presented the time-dependent flow on a period of time for Rac2 = 4.58×104 (Figure 7) and illustrated the structure of the flow by the streamlines evolution over time (τa, τb, τc, τd, τe and τf). It was found that the flow remains unstable and unicellular streamlines which looks oscillatory (Figures 8a-g). We can see two cells counter-rotating, the left one is bigger than the other, later they will have the same shape in the symmetry case. And then the right one becomes bigger. For the isotherms, we can see that the head of the mushroom move slowly between the left side and the central axe of the cavity. In the next time period it will move from the central axe to the right side (Figure 9a-g).

Figure 7. Time-evolution of the vertical velocity component in τ =1, for Rac2 = 4.58 × 104

(g) At time tg.

Figure 8. Streamlines of oscillatory convection (Rac2 = 4.58× 104), at several time (τa, τb, τc, τd, τe and τf)

4.3 Flow bifurcation

The bifurcation flow in the cavity and the different forms that may have: bifurcated cells, unicellular, two vertical cells, two horizontal cells are presented in this section. This study was done for values of Rac-1 up to 106. We can observe different bifurcation flow structures shown in Table 7.

Figure 10 shows the logarithmic bifurcation diagram of Rayleigh-Bénard convection. The figure also summarizes the streamlines and isothermal line according the Rayleigh number. The first bifurcation is observed for Rac-1, which represents the transition from conduction to convection (onset of convection). This bifurcation is represented by the one principal cell in the middle of the cavity and with a clockwise flow direction, followed by another bifurcation that appears for the value of Rac-1= 1.9×104. In this case, the basic unit cell of our convection (UC) is divided in two vertical cells (2VC), one cell with a clockwise direction and the other has a counter-clockwise, to move on impulse isothermal lines. In this case the particles do not lie near the adiabatic wall but in the center of the hot wall and the cold particles come down almost close of the two adiabatic walls. It will create a mushroom shape. Rayleigh-Bénard convection is very sensitive area which has a lot of bifurcations, as a result the two cells bifurcation, at Ra = 5×104, another single bifurcation cell appears similar to the first. Finally, for the value Ra = 4.2×105, there are two cells. But in this case, they are horizontal (2HC), hot and cold particles; go up and down on one and the same adiabatic wall.

(g) At time tg.

Figure 9. Isothermal lines of oscillatory convection (Rac2 = 4.58× 104), at several time (τa, τb, τc, τd, τe and τf)

Table 7. Flow Structure for different Rayleigh numbers

Ra

Cell number and disposition

Flow regime

[2.5×103, 1.8×104]

Unicellular (UC)

Steady state

[1.9×104, 4.57×104]

Two vertical cells (2VC)

Steady state

[4.58×104, 4.99×104]

Two vertical cells (2VC)

periodic oscillatory

[5×104, 4.1×105]

Unicellular (UC)

Steady state

[4.2×105, 106]

Two horizontal cells (2HC)

periodic oscillatory

Figure 10. Logarithmic bifurcation diagram of Rayleigh-Bénard convection with streamlines and isothermal line

5. Conclusions

The study has been made ​​for the Rayleigh number varied from 2.5×103 to 106 in order to define the Rac-1 corresponding to the onset of natural convection for different aspect ratios of the cavity. The transition threshold regime laminar-oscillatory convection which is defined by Rac2 is determined. In the range studied there were five modes of bifircation of Rayleigh-Bénard convection: unicellular (UC), two vertical cells (2VC) in steady state, two vertical cells (2VC) in oscillatory periodic state, and two horizontal cells (2HC) in oscillatory periodic state.

Acknowledgments

This paper was supported by the Research Team: Heat Transfer and Renewable Energy (HTRE) attached to the Mechanical Laboratory (LME), University of Laghouat that approved in 2009 under the number: 222 in 07/13/2009.

Nomenclature

A

Aspect ratio

g

k

gravitational acceleration, m.s-2

thermal conductivity, W.m-1. K-1

H

 dimensional height of the enclosure  (m)

Nu

Nusselt number

P

 dimensionless pressure 

Pr

Prandtl number

Ra

Rayleigh number

Rac-1

first critical Rayleigh number (onset of convection)

Rac-2

second critical Rayleigh number (oscillatory convective flow)

t

dimensional time (s)

T

dimensional temperature (K)

u,v

velocity components (m.s-1)

U,V

dimensionless velocity components

x,y

dimensional coordinate space (m)

X,Y

dimensionless coordinate space

W

dimensional width of the enclosure (m)

Greek symbols

 

α

thermal diffusivity (m2.s-1)

β

coefficient of thermal expansion at constant pressure (K-1)

λ

thermal conductivity of the fluid (w.m-1K-)

υ

kinematics viscosity (m2.s-1)

μ

dynamic viscosity (kg.s-1.m-1)

ρ

density (kg.m3)

τ

dimensionless time

θ

dimensionless temperature

  References

[1] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Oxford: Clarendon Press, 1961. DOI: 10.1063/1.3058072.

[2] R. M. Clever and F. H. Busse, “Transition to time-dependent convection”, J. Fluid Mech., vol. 65, no. 4, pp. 625-645, 1974. DOI: 10.1017/S0022112074001571.

[3] Dr. Kays A. Al-Tae'y, “Natural convection in a laterally and volumetrically heated square enclosure with time-periodic boundary conditions”. Int. J. Heat and Technology, vol. 27, no. 2, pp.75-83, 2009.

[4] N. B. Cheikh, B. B. Beya and T. Lili, “Aspect ratio effect on natural convection flow in a cavity submitted to a periodical temperature boundary”, J. Heat Transfer. vol. 129, pp. 1060-1068, 2007. DOI: 10.1115/1.2728908.

[5] V.A.F. Costa, M.S.A. Oliveira and A.C.M. Sousa, “Control of laminar natural convection in differentially heated square enclosures using solid inserts at the corners”, Int. J. Heat Mass Transfer. vol. 46, pp. 3529-3537, 2003. DOI: 10.1016/S0017-9310 (03) 00141-8U.

[6] G. Accary and I. Raspo, “A 3D finite volume method for the prediction of a supercritical fluid buoyant flow in a differentially heated cavity”, Int. J. Heat Mass Transfer, vol. 35, pp. 1316-1331, 2006. DOI: 10.1016/j.compfluid.2005.05.004.

[7] S. M. Aminossadati, B. Ghasemi, “Computational modelling of heat transfer in a partially partitioned enclosure”, Int. J. Heat and Technology, vol. 26, no. 2, pp. 61-67, 2008.

[8] G. De Vahl Davis, “Natural convection of air in a square cavity: a benchmark numerical solution”, Int. J. Numer. Methods Fluids, vol. 3, pp. 249-264, 1983.  DOI: 10.1002/fld.1650030305.

[9] M. Lappa, “Some considerations about the symmetry and evolution of chaotic Rayleigh–Bénard convection: The flywheel mechanism and the “wind” of turbulence”, C. R. Mécanique, vol. 339, pp. 261-268, 2011. DOI: 10.1016/j.crme.2011.05.002.

[10] D. Venturi, “Stochastic bifurcation analysis of Rayleigh–Bénard convection”, J. Fluid Mech., vol. 650, pp. 391-413, 2010. DOI: 10.1017/S0022112009993685.

[11] T. Mas, “Dynamic bifurcation and stability in the Rayleigh-Benard convection”, Commun. Math. Sci. vol. 2 no. 2, pp. 159-183, 2004.

[12] P. H. Kao, “Simulating oscillatory flows in Rayleigh–Bénard convection using the Lattice Boltzmann Method”, Int. J. of Heat and Mass Trans., vol. 50, 2007. DOI: 10.1016/j.ijheatmasstransfer.2007.01.035.

[13] D. Angeli, A. Pagano, M. Corticelli, A. Fichera, Giovanni S. Barozzi, “Bifurcations of natural convection flows from an enclosed cylindrical heat source,” Frontiers in Heat and Mass Trans, 2, 2011. DOI: 10.5098/hmt.v2.2.3003.

[14] A. Raji, M. Hasnaoui, M. Firdouss and C. Ouardi, “Natural convection heat transfer enhancement in a square cavity periodically cooled form above”, Numerical Heat Transfer, Part A, vol. 63 no. 7, pp. 511-533, 2013. DOI: 10.1080/10407782.2013.733248.

[15] B. Ghernaout, S. Bouabdallah, M. Teggar and H. Benniche, “Double diffusive natural convection in binary mixture under the effect of external magnetic field: steady and oscillatory state”, Int. J. Heat and Technology, vol. 33, no. 4, pp. 11-18, 2015. DOI: 10.18280/ijht.330402.

[16] S. Patankar, Numerical Heat Transfer and Fluid Flow, New York: Mc-Graw-Hill, 1980.

[17] A. Y. Gelfgat, “Different modes of Rayleigh-Bénard instability in two- and three-dimensional rectangular enclosures”, J. Comp. Phys, vol. 156, pp. 300-324, 1999. DOI: 10.1006/jcph.1999.6363.

[18] O. Shishkina, “Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution,” New J. Phys., vol. 12, pp. 1-9, 2010. DOI: 10.1088/1367-2630/12/7/075022.

[19] O. Turan, “Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure,” Journal of Non-Newtonian Fluid Mech., vol. 171, pp.83-96, 2012. DOI: 10.1016/j.jnnfm.2012.01.006.