# Numerical Simulation and Stability Study of Natural Convection in an Inclined Rectangular Cavity

Authors: Hua-Shu Dou; Gang Jiang; Chengwang Lei
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This paper examines the process of instability of natural convection in an inclined cavity based on numerical simulations. The energy gradient method is employed to analyze the physics of the flow instability in natural convection. It is found that the maximum value of the energy gradient function in the flow field correlates well with the location where flow instability occurs. Meanwhile, the effects of the flow time, the plate length, and the inclination angle on the instability have also been discussed. It is observed that the locations of instabilities migrate right as the flow time increased. With the increase of plate length, the onset time of the instability on the top wall of the cavity decreases gradually and the locations of instabilities move to the right side. Furthermore, the locations of instability move left with the increase of the inclination angle in a certain range. However, these positions move right as the accumulation of the heat flux is restrained in the lower left corner of the cavity once the inclination angle exceeds a certain range.


---

## Body

## 1.  Introduction

Transient natural convection flows in a cavity are common in industrial applications such as in heat exchangers, solar collectors, and nuclear reactors and in our daily life such as in light emitting diode (LED) street lights, computers, and mobile phones. Actually, some engineering problems are related to the cases with an inclined cavity. Yet, natural convection adjacent to an inclined plate has received less attention than the classic cases of vertical and horizontal plates. Natural convection heat transfer is regarded as one of the three basic forms of heat transfer, and there exists some difference of heat transfer efficiency between different manners of natural convection. The achieved results suggested that turbulent flow has a more powerful ability to transfer heat in natural convection compared to laminar flow.Patterson and Imberger [1] carried out extensive investigations on the transient behavior of natural convection of a two-dimensional rectangular cavity in which the two opposing vertical sidewalls are simultaneously heated and cooled by an equal amount. The authors present several flow regimes of the flow development of the boundary layer which were named as “conduction regime,” “stable convection regime,” and “unstable convection regime,” respectively. These studies are based on the relative values of the Rayleigh number Ra, the Prandtl number Pr, and the aspect ratio of the cavity. It was found that the flow of natural convection in an inclined cavity loses its stability by forming longitudinal vortices. Sparrow and Husar [2] made experiments on inclined plates to reveal the presence of cellular secondary flows superposed upon the natural convection main flow and believed that these longitudinal vortices were the first stage of the laminar-turbulent transition process. Haaland and Sparrow [3] tried to use linear stability theory to predict a critical angle (or range of angles) at which the disturbances change their character from travelling waves to longitudinal vortices. They finally found that it is difficult to predict the critical angle theoretically by treating each type of disturbance separately.Lloyd and Sparrow [4] made some investigations with the aim of establishing the relationship between the inclination angle and the nature of the instability. They tried to determine accurate quantitative information on the angular dependence of the Rayleigh number for instability. They found that waves were the mode of instability for inclination angles of less than 14 deg (relative to the vertical). Further, when the inclination angle is beyond 17 deg, the instability was characterized by longitudinal vortices. The range between 14 deg and 17 deg was a zone of continuous transition, with the two modes of instability coexisting.Ganesan and Palani [5] proposed to study the natural convection effects on impulsively started inclined plate with heat and mass transfer by an implicit finite difference scheme of the Crank-Nicolson type. In order to access the accuracy of the numerical results, they compared their study with available exact solution of Moutsoglou and Chen [6] and achieved a good agreement. In addition, they observed that local wall shear stress decreases as the angle of inclination decreases.Xu et al. [7] researched transient natural convection flows around a thin fin on the sidewall of differentially heated cavity, which illustrated that the fluid boundary layer adjacent to a vertical thermal wall included three sublayers, and these sublayers were determined by different dynamic and energy balances. When the initial time t1 is sufficiently small, the balance was determined by viscous term and buoyancy term which yields the inner viscous layer. Subsequently, as the flow time increases to t2, the balance is still determined by viscous term and buoyancy term, which yields a viscous layer within the thermal boundary layer but outside the inner viscous layer. At the same time, with the formation of inner viscous layer and viscous layer, there is a balance between the conduction term and the unsteady term which yields the thermal boundary layer. Saha et al. [8] studied natural convection of an inclined flat plate under a sudden cooling condition, and they found that the cold boundary layer adjacent to the plate is potentially unstable, if the Rayleigh number Ra exceeded a critical value. In other words, the boundary layer would be always stable if Ra was below a threshold.Said et al. [9] made some numerical investigations on turbulent natural convection in a parallel-walled channel which is inclined with respect to gravity. They found that the channel overall average Nusselt number was reduced as the inclination angle was increased. The rate of reduction in the overall Nusselt number decreases as the Rayleigh number increases. Additionally, they observed that the local Nusselt number was much higher along the lower wall of the horizontal channel where cold air entered in comparison to the upper wall where hot air existed at both channel openings.Lin [10] presented a numerical experiment for the onset and its linear development of longitudinal vortices in natural convection over inclined plates. It was observed that the critical Grashof number increases with the increase of the inclination angle while the effect of inclination angle on the Nusselt number is less pronounced when the value of the inclined angle increases.Iyer and Kelly [11] did not agree with the conclusions made by Haaland and Sparrow [3]. They found that these experiments in [3] were not sensitive enough to detect the first instabilities predicted by theoretical analysis. Iyer and Kelly used a spatial linear stability analysis with the parallel flow assumption to examine the formation and growth of both wave instabilities and longitudinal vortices and attempted to find a correlation between experimental and theoretical results by finding the total amplification between the earliest disturbances and the observed disturbances.Thus, there are many factors to affect the natural convection heat transfer efficiency, in particular, the difference of laminar flow and turbulence. In order to achieve a flow field where the flow is in turbulent state, it is helpful to ascertain the locations where instabilities could occur and where instability would occur firstly.After almost 20 years of work, Dou et al. [12–18] suggested a new approach to analyze flow instability and turbulence transition based on “energy gradient method.” This approach is different from the linear stability theory, the weak nonlinear stability theory, the secondary instability theory, and the energy method. This approach explains the mechanism of flow instability from physics and derives the criteria of turbulence transition. The theoretical results are in agreement with the experimental data of the pipe Poiseuille flow, plane Poiseuille flow, plane Couette flow, Taylor-Couette flow, boundary layer flow, and so on.In this study, numerical simulation is used to obtain the flow field at various geometrical and flow parameters. Then, the energy gradient method is used to investigate the physical mechanism of flow instability in natural convection. The paper is divided into two sections. In the first section, the same computational geometries and numerical scheme of Saha et al. [8] are used to validate the numerical method used in this study. In the second section, energy gradient method is briefly introduced to calculate the value of K in the whole flow field and make some further investigations on the effects of the flow time, plate length, and the inclination angle on flow instability of natural convection. This is the first time for the energy gradient method to be used in natural convection.


## 2. Computational Geometry and Numerical Procedures

### 2.1. Computational Geometry

The computational geometry is shown in Figure1. We consider the top wall as a cooled inclined flat plate where the temperature is fixed at Tc. Initially, the fluid temperature in the domain is T0 which is higher than Tc. The four side walls of the domain are rigid and nonslip. Except for the top plate, all the three other walls of the rectangle cavity are assumed to be insulated. The length and the width of the cavity are defined as L and W, respectively. The inclination angle of the cavity is expressed with θ.Figure 1
Schematic of the inclined cavity.


### 2.2. Numerical Procedures

#### 2.2.1. Governing Equations

The development of natural convection adjacent to an inclined cavity is governed by the following two-dimensional Navier-Stokes and energy equations, and these equations are based on the Boussinesq approximation:(1)∂u∂x+∂v∂y=0,∂u∂t+u∂u∂x+v∂u∂y=-1ρ∂p∂x+ν(∂2u∂x2+∂2u∂y2)+gβsinθ(T-T0),∂v∂t+u∂v∂x+v∂v∂y=-1ρ∂p∂y+ν(∂2v∂x2+∂2v∂y2)+gβcosθ(T-T0),∂T∂t+u∂T∂x+v∂T∂y=k(∂2T∂x2+∂2T∂y2),
where x and y are the horizontal and vertical coordinates with origin at the lower left corner of the cavity, t is the time, T is the temperature, p is the pressure, u and v are the velocity components in the x and y directions, g is the acceleration due to gravity, β is the coefficient of thermal expansion,ρ is the fluid density, k is the thermal diffusivity, and υ is the kinematic viscosity.


#### 2.2.2. Numerical Scheme

The governing equations (1) are implicitly solved using a finite-volume SIMPLE scheme, with the QUICK scheme approximating the advection term. The diffusion terms are discretized using central differencing with second-order accurate. A second-order implicit time- marching scheme will be used for the unsteady term. The discretized equations are iterated with specified underrelaxation factors. The boundary condition of the left wall and the right wall is∂T/∂x=0, u=v=0; the boundary condition of the bottom is:∂T/∂y=0, u=v=0. In addition, we should note that the flow is unsteady, and all the plots depend on the time.


## 2.1. Computational Geometry

The computational geometry is shown in Figure1. We consider the top wall as a cooled inclined flat plate where the temperature is fixed at Tc. Initially, the fluid temperature in the domain is T0 which is higher than Tc. The four side walls of the domain are rigid and nonslip. Except for the top plate, all the three other walls of the rectangle cavity are assumed to be insulated. The length and the width of the cavity are defined as L and W, respectively. The inclination angle of the cavity is expressed with θ.Figure 1
Schematic of the inclined cavity.


## 2.2. Numerical Procedures

### 2.2.1. Governing Equations

The development of natural convection adjacent to an inclined cavity is governed by the following two-dimensional Navier-Stokes and energy equations, and these equations are based on the Boussinesq approximation:(1)∂u∂x+∂v∂y=0,∂u∂t+u∂u∂x+v∂u∂y=-1ρ∂p∂x+ν(∂2u∂x2+∂2u∂y2)+gβsinθ(T-T0),∂v∂t+u∂v∂x+v∂v∂y=-1ρ∂p∂y+ν(∂2v∂x2+∂2v∂y2)+gβcosθ(T-T0),∂T∂t+u∂T∂x+v∂T∂y=k(∂2T∂x2+∂2T∂y2),
where x and y are the horizontal and vertical coordinates with origin at the lower left corner of the cavity, t is the time, T is the temperature, p is the pressure, u and v are the velocity components in the x and y directions, g is the acceleration due to gravity, β is the coefficient of thermal expansion,ρ is the fluid density, k is the thermal diffusivity, and υ is the kinematic viscosity.


### 2.2.2. Numerical Scheme

The governing equations (1) are implicitly solved using a finite-volume SIMPLE scheme, with the QUICK scheme approximating the advection term. The diffusion terms are discretized using central differencing with second-order accurate. A second-order implicit time- marching scheme will be used for the unsteady term. The discretized equations are iterated with specified underrelaxation factors. The boundary condition of the left wall and the right wall is∂T/∂x=0, u=v=0; the boundary condition of the bottom is:∂T/∂y=0, u=v=0. In addition, we should note that the flow is unsteady, and all the plots depend on the time.


## 2.2.1. Governing Equations

The development of natural convection adjacent to an inclined cavity is governed by the following two-dimensional Navier-Stokes and energy equations, and these equations are based on the Boussinesq approximation:(1)∂u∂x+∂v∂y=0,∂u∂t+u∂u∂x+v∂u∂y=-1ρ∂p∂x+ν(∂2u∂x2+∂2u∂y2)+gβsinθ(T-T0),∂v∂t+u∂v∂x+v∂v∂y=-1ρ∂p∂y+ν(∂2v∂x2+∂2v∂y2)+gβcosθ(T-T0),∂T∂t+u∂T∂x+v∂T∂y=k(∂2T∂x2+∂2T∂y2),
where x and y are the horizontal and vertical coordinates with origin at the lower left corner of the cavity, t is the time, T is the temperature, p is the pressure, u and v are the velocity components in the x and y directions, g is the acceleration due to gravity, β is the coefficient of thermal expansion,ρ is the fluid density, k is the thermal diffusivity, and υ is the kinematic viscosity.


## 2.2.2. Numerical Scheme

The governing equations (1) are implicitly solved using a finite-volume SIMPLE scheme, with the QUICK scheme approximating the advection term. The diffusion terms are discretized using central differencing with second-order accurate. A second-order implicit time- marching scheme will be used for the unsteady term. The discretized equations are iterated with specified underrelaxation factors. The boundary condition of the left wall and the right wall is∂T/∂x=0, u=v=0; the boundary condition of the bottom is:∂T/∂y=0, u=v=0. In addition, we should note that the flow is unsteady, and all the plots depend on the time.


## 3. Criteria of Instability Based on the Rayleigh Number

The natural convection boundary layer adjacent to an inclined cavity is subjected to sudden cooling boundary condition which yields a thermal boundary layer along the wall. At the same time with the formation of the thermal boundary layer, the viscous boundary layer is formed which is determined by the balance between viscous term and inertial term. Meanwhile, the velocity inside the boundary layer develops, governed by the balance of viscous and inertial terms with the buoyancy term. It is known from [7, 19] that the thickness of thermal boundary layer, the velocity inside the thermal boundary layer, and the thickness of the viscous boundary layer are related to the balance between conduction and advection terms.We can find in [8] that, when the balance between conduction and advection terms plays a leading role in exchange of momentum and energy, the steady-state scales of the boundary layer (ts) can be achieved:
(2)ts~(1+Pr)1/2(1+tan2θ)1/2tanθ·Ra1/2Pr1/2(L2sin2θk),
where Ra is the Rayleigh number, Ra=gβΔT(Lsinθ)3/υk and Pr is the Prandtl number, Pr=υ/k. In the present case, the thermal boundary layer is bounded by a rigid surface of the plate and a cold air layer, which is equivalent to the free-rigid boundary configuration [20–24], in which there is a critical Rayleigh numberRac=1106.5. Also, there exists a critical time scale tB for the onset of thermal layer instability at a given Ra. If t>tB, the instability will set in before the growth of the thermal boundary layer completes. On the other hand, if t<tB, the instability will never occur no matter how much time it iterates. The critical time scale tB is described as follows [8]:
(3)tB=(RacRa)2/3L2sin2θk.
The ratio between the growth time of thermal boundary layer ts and the critical time scale tB described as follows [8],
(4)tstB=[(1+Pr)3(1+tan2θ)3RaPr3tan6θ•Rac4]1/6.
Absolutely, if Ra>Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3, that is, ts>tB, the instability will set in before the growth of the thermal boundary layer completes, or else, if Ra<Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3, that is, ts<tB, the instability will never occur.


## 4. Validation of Numerical Methods

### 4.1. Numerical Scheme Test

In order to verify the accuracy of the numerical scheme used in this study, the computational geometry of Saha et al. [8] is considered. In the following numerical simulations, the condition Pr<1 should be satisfied. If Pr>1, Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3 is very large. Thus, it is difficult to observe instability if the value of Ra is low.Figure2 shows the simulated temperature contours with three different cavity dimensions at 1000 s, which is taken from Saha et al. [8]. The lengths are 5.4 m, 10.8 m, and 16.2 m, respectively, and the width is 0.6 m.Temperature contours of three different geometries in [8].
(a)
Length = 5.4 m
(b)
Length = 10.8 m
(c)
Length = 16.2 mFigure3 presents the current results of numerical simulations at 1000 s for the cavity sizes in Figure 2. Comparing Figure 3 with Figure 2, it can be found that the present simulations are in agreement with those in [8], which demonstrates that the numerical scheme is reliable and accurate.Temperature contours of three different geometries (present calculation).
(a)
L = 5.4 m
(b)
L = 10.8 m
(c)
L = 16.2 m


### 4.2. Grid Independence Test

Grid independence will be examined with three different mesh sizes for the same domain. The three mesh sizes are as follows (Figure4): (a) 100×300 with 30000 cells, 60400 faces, and 30401 nodes; (b) 150×450 with 67500 cells, 135600 faces, and 60801 nodes; (c) 200×600 with 120000 cells, 240800 faces, and 120801 nodes.Comparison of the calculated results with three different mesh sizes: (a)Δx=0.054, (b) Δx=0.036, and (c) Δx=0.027.
(a)
100 * 300t = 300 s
(b)
150 * 450t = 300 s
(c)
200 * 600t = 300 sThe simulation results of temperature contours at 300 s with these three meshes are shown in Figure4. The temperature versus time with three different mesh sizes is shown in Figure 5, which is recorded at the same monitor point (2, 1). It can be seen from Figures 4 and 5 that mesh convergence has been achieved with these meshes.Figure 5
Temperature profile calculated with three different grid sizes.


### 4.3. Test of Simulation Results

Four different configurations are listed in Table1, and substituting these corresponding data into (2) and (3), we will get ts and tB, respectively. All the numerical models in Figure 5 satisfy the following conditions: Pr=0.72, θ=5.71°, and k=2.04×10-5. Figures 6(a2), 6(b2), 6(c2), and 6(d2) are partial enlarged drawings of the middle portion of the plate in Figures 6(a1), 6(b1), 6(c1), and 6(d1), respectively.Table 1
Values of configuration dimensions, Ra,ts, and tB shown in Figure 6.
Runs
Configuration dimension
Ra
t
s
t
B
(a1)
W
(
0.5
m
)
×
L
(
5
m
)
5.33
×
1
0
4
815.6 s
915.5 s
(b1)
W
(
0.5
m
)
×
L
(
5
m
)
5.33
×
1
0
5
257.9 s
197.2 s
(c1)
W
(
0.5
m
)
×
L
(
2
m
)
6.83
×
1
0
3
364.7 s
576.8 s
(d1)
W
(
0.5
m
)
×
L
(
15
m
)
2.88
×
1
0
6
998.9 s
576.8 sFigure 6
Four groups of geometries with different Ra: (a1)Ra=5.335e4, (b1) Ra=5.335e5, (c1) Ra=6.83e3, and (d1) Ra=2.88e6.By comparing these four groups of pictures in Figure6, some conclusions may be written as follows.(1)
It is observed that there is an “end effect” at the top right corner of the enclosure in each case, and this phenomenon is affected by the configuration and inclination angle of the cavity.(2)
The iterative time reachest=ts both in Figures 6(a1) and 6(c1), and it is very clear that the flow is still stable which can be seen from the partial enlarger of the middle portion of the plate.(3)
The oscillation of the curves from the isotherms in Figures6(b1) and 6(d1) indicates the onset of instability. It is clearly seen that the thermal boundary layer travels in waveform rather than in a smooth manner.(4)
It can be found thatts<tB is satisfied for the simulated results shown in Figures 6(a1) and 6(c1), and it is observed that the thermal boundary layer is stable at these cases. While, in Figures 6(b1) and 6(d1), ts>tB, the stability criterion is violated. Thus, the flow instability occurs at a proper time. These results are in good agreement with the predictions of criteria of instability based on Ra.Nevertheless, we observe some unique phenomena from these previous instability cases. Firstly, instabilities do not occur in the whole flow field. Secondly, there exists a distinctive time difference of instabilities at different locations. In other words, it is expected to employ a theory to predict the locations where instabilities could occur and positions where instability would take place firstly. In the following study, we will briefly introduce the energy gradient method to analyze the above phenomena observed.


## 4.1. Numerical Scheme Test

In order to verify the accuracy of the numerical scheme used in this study, the computational geometry of Saha et al. [8] is considered. In the following numerical simulations, the condition Pr<1 should be satisfied. If Pr>1, Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3 is very large. Thus, it is difficult to observe instability if the value of Ra is low.Figure2 shows the simulated temperature contours with three different cavity dimensions at 1000 s, which is taken from Saha et al. [8]. The lengths are 5.4 m, 10.8 m, and 16.2 m, respectively, and the width is 0.6 m.Temperature contours of three different geometries in [8].
(a)
Length = 5.4 m
(b)
Length = 10.8 m
(c)
Length = 16.2 mFigure3 presents the current results of numerical simulations at 1000 s for the cavity sizes in Figure 2. Comparing Figure 3 with Figure 2, it can be found that the present simulations are in agreement with those in [8], which demonstrates that the numerical scheme is reliable and accurate.Temperature contours of three different geometries (present calculation).
(a)
L = 5.4 m
(b)
L = 10.8 m
(c)
L = 16.2 m


## 4.2. Grid Independence Test

Grid independence will be examined with three different mesh sizes for the same domain. The three mesh sizes are as follows (Figure4): (a) 100×300 with 30000 cells, 60400 faces, and 30401 nodes; (b) 150×450 with 67500 cells, 135600 faces, and 60801 nodes; (c) 200×600 with 120000 cells, 240800 faces, and 120801 nodes.Comparison of the calculated results with three different mesh sizes: (a)Δx=0.054, (b) Δx=0.036, and (c) Δx=0.027.
(a)
100 * 300t = 300 s
(b)
150 * 450t = 300 s
(c)
200 * 600t = 300 sThe simulation results of temperature contours at 300 s with these three meshes are shown in Figure4. The temperature versus time with three different mesh sizes is shown in Figure 5, which is recorded at the same monitor point (2, 1). It can be seen from Figures 4 and 5 that mesh convergence has been achieved with these meshes.Figure 5
Temperature profile calculated with three different grid sizes.


## 4.3. Test of Simulation Results

Four different configurations are listed in Table1, and substituting these corresponding data into (2) and (3), we will get ts and tB, respectively. All the numerical models in Figure 5 satisfy the following conditions: Pr=0.72, θ=5.71°, and k=2.04×10-5. Figures 6(a2), 6(b2), 6(c2), and 6(d2) are partial enlarged drawings of the middle portion of the plate in Figures 6(a1), 6(b1), 6(c1), and 6(d1), respectively.Table 1
Values of configuration dimensions, Ra,ts, and tB shown in Figure 6.
Runs
Configuration dimension
Ra
t
s
t
B
(a1)
W
(
0.5
m
)
×
L
(
5
m
)
5.33
×
1
0
4
815.6 s
915.5 s
(b1)
W
(
0.5
m
)
×
L
(
5
m
)
5.33
×
1
0
5
257.9 s
197.2 s
(c1)
W
(
0.5
m
)
×
L
(
2
m
)
6.83
×
1
0
3
364.7 s
576.8 s
(d1)
W
(
0.5
m
)
×
L
(
15
m
)
2.88
×
1
0
6
998.9 s
576.8 sFigure 6
Four groups of geometries with different Ra: (a1)Ra=5.335e4, (b1) Ra=5.335e5, (c1) Ra=6.83e3, and (d1) Ra=2.88e6.By comparing these four groups of pictures in Figure6, some conclusions may be written as follows.(1)
It is observed that there is an “end effect” at the top right corner of the enclosure in each case, and this phenomenon is affected by the configuration and inclination angle of the cavity.(2)
The iterative time reachest=ts both in Figures 6(a1) and 6(c1), and it is very clear that the flow is still stable which can be seen from the partial enlarger of the middle portion of the plate.(3)
The oscillation of the curves from the isotherms in Figures6(b1) and 6(d1) indicates the onset of instability. It is clearly seen that the thermal boundary layer travels in waveform rather than in a smooth manner.(4)
It can be found thatts<tB is satisfied for the simulated results shown in Figures 6(a1) and 6(c1), and it is observed that the thermal boundary layer is stable at these cases. While, in Figures 6(b1) and 6(d1), ts>tB, the stability criterion is violated. Thus, the flow instability occurs at a proper time. These results are in good agreement with the predictions of criteria of instability based on Ra.Nevertheless, we observe some unique phenomena from these previous instability cases. Firstly, instabilities do not occur in the whole flow field. Secondly, there exists a distinctive time difference of instabilities at different locations. In other words, it is expected to employ a theory to predict the locations where instabilities could occur and positions where instability would take place firstly. In the following study, we will briefly introduce the energy gradient method to analyze the above phenomena observed.


## 5. Application of Energy Gradient Method

### 5.1. Energy Gradient Method

From the classical theory of the Brownian motion, the fluid particles exchange energy and momentum all the time via collisions. The fluid particle will collide with other particles in transverse directions as it flows along its streamline, and this particle would obtain energy expressed asΔE after many cycles; at the same time, the particle would drop energy due to viscosity along the streamline; with the same periods, the energy loss expressed as ΔH would be considerable. Consequently, there exists a critical value of the ratio of ΔE and ΔH, above which the particle would leave its equilibrium by moving to a new streamline with higher energy or lower energy and below which the particle would not leave its streamline for its oscillation would be balanced by the viscosity along the streamline. Making reference to [13–18], we can express the criteria of instability as follows:
(5)F=ΔEΔH=((∂E/∂n)(2A-/π))((∂H/∂s)(π/ωd)u)=2π2KA-ωdu=2π2Kvm′u<Const,
where
(6)K=∂E/∂n∂H/∂s.
Here, F is a function of coordinates which expresses the ratio of the energy gained in a half period by the particle and the energy loss due to viscosity in the half period. K is a dimensionless field variable (function) and expresses the ratio of transversal energy gradient and the rate of the energy loss along the streamline. E=p+1/2ρV2 is the kinetic energy per unit volumetric fluid, s is along the streamwise direction, and n is along the transverse direction. H is the loss of the total mechanical energy per unit volumetric fluid along the streamline for finite length, which can be calculated from the Navier-Stokes equations. Further, ρ is the fluid density, u is the streamwise velocity of main flow, A- is the amplitude of the disturbance distance, ωd is the frequency of the disturbance, and vm′=A-ωd is the amplitude of the disturbance of velocity.Equation (5) represents the criteria of instability; that is, if F exceeds its threshold, the fluid particle would lose its stability or the instability would never set in. Equation (6) represents characteristic of instability; that is, it shows the most dangerous positions in a flow field, and it indicates that instability would occur firstly at the position with the maximum value of K once instability sets in.


### 5.2. Criterion of Instability Based on Energy Gradient in Natural Convection

Dou and Phan-Thien [25] proposed an energy gradient theory which describes the rules of fluid material stability from the viewpoint of energy field and can be considered as a supplement to the Newtonian mechanics. They claimed that the instability of natural convection could not be resolved by Newton’s three laws, for the reason that a material system moving in some cases is not simply due to the role of forces. This method does not attribute the Rayleigh-Benard problem to forces, but to energy gradient. It postulates that when the fluid is placed on a horizontal plate and it is heated from below, the fluid density in the bottom becomes low which leads to energy gradient ∂E/∂y>0 along y-coordinate. Only when ∂E/∂y is larger than a critical value, will the flow become unstable, and then fluid cells of vorticities will be formed. This conclusion is in accordance with the former criteria of instability.When the fluid is placed on an inclined plate or a box, the criterion of natural convection can be written as (Figure7)
(7)K=(∂E∂x)2+(∂E∂y)2.
In present study, neglecting the influence of the gravity, we get E~p0.Figure 7
Calculation of the value ofK.Besides the ability to predict whether instability could occur in natural convection [25], energy gradient method has another two functions in natural convection when instabilities had occurred. It can predict the locations where instabilities could occur and the position where instability would take place firstly via the value of K. In the proposed method, the flow is expected to be more unstable in the area with high value of K than that in the area with low value of K. This is the distinctive difference between energy gradient theory and the criteria of instability based on the Rayleigh number and is also the focus of investigation in this paper.We should clarify that the value ofK derived from energy gradient method can only be applied in situations where instabilities could occur. Thus, all the subsequent numerical simulations should satisfy the condition that Ra is larger than Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3, in order to study the cases of instability occurrence.


## 5.1. Energy Gradient Method

From the classical theory of the Brownian motion, the fluid particles exchange energy and momentum all the time via collisions. The fluid particle will collide with other particles in transverse directions as it flows along its streamline, and this particle would obtain energy expressed asΔE after many cycles; at the same time, the particle would drop energy due to viscosity along the streamline; with the same periods, the energy loss expressed as ΔH would be considerable. Consequently, there exists a critical value of the ratio of ΔE and ΔH, above which the particle would leave its equilibrium by moving to a new streamline with higher energy or lower energy and below which the particle would not leave its streamline for its oscillation would be balanced by the viscosity along the streamline. Making reference to [13–18], we can express the criteria of instability as follows:
(5)F=ΔEΔH=((∂E/∂n)(2A-/π))((∂H/∂s)(π/ωd)u)=2π2KA-ωdu=2π2Kvm′u<Const,
where
(6)K=∂E/∂n∂H/∂s.
Here, F is a function of coordinates which expresses the ratio of the energy gained in a half period by the particle and the energy loss due to viscosity in the half period. K is a dimensionless field variable (function) and expresses the ratio of transversal energy gradient and the rate of the energy loss along the streamline. E=p+1/2ρV2 is the kinetic energy per unit volumetric fluid, s is along the streamwise direction, and n is along the transverse direction. H is the loss of the total mechanical energy per unit volumetric fluid along the streamline for finite length, which can be calculated from the Navier-Stokes equations. Further, ρ is the fluid density, u is the streamwise velocity of main flow, A- is the amplitude of the disturbance distance, ωd is the frequency of the disturbance, and vm′=A-ωd is the amplitude of the disturbance of velocity.Equation (5) represents the criteria of instability; that is, if F exceeds its threshold, the fluid particle would lose its stability or the instability would never set in. Equation (6) represents characteristic of instability; that is, it shows the most dangerous positions in a flow field, and it indicates that instability would occur firstly at the position with the maximum value of K once instability sets in.


## 5.2. Criterion of Instability Based on Energy Gradient in Natural Convection

Dou and Phan-Thien [25] proposed an energy gradient theory which describes the rules of fluid material stability from the viewpoint of energy field and can be considered as a supplement to the Newtonian mechanics. They claimed that the instability of natural convection could not be resolved by Newton’s three laws, for the reason that a material system moving in some cases is not simply due to the role of forces. This method does not attribute the Rayleigh-Benard problem to forces, but to energy gradient. It postulates that when the fluid is placed on a horizontal plate and it is heated from below, the fluid density in the bottom becomes low which leads to energy gradient ∂E/∂y>0 along y-coordinate. Only when ∂E/∂y is larger than a critical value, will the flow become unstable, and then fluid cells of vorticities will be formed. This conclusion is in accordance with the former criteria of instability.When the fluid is placed on an inclined plate or a box, the criterion of natural convection can be written as (Figure7)
(7)K=(∂E∂x)2+(∂E∂y)2.
In present study, neglecting the influence of the gravity, we get E~p0.Figure 7
Calculation of the value ofK.Besides the ability to predict whether instability could occur in natural convection [25], energy gradient method has another two functions in natural convection when instabilities had occurred. It can predict the locations where instabilities could occur and the position where instability would take place firstly via the value of K. In the proposed method, the flow is expected to be more unstable in the area with high value of K than that in the area with low value of K. This is the distinctive difference between energy gradient theory and the criteria of instability based on the Rayleigh number and is also the focus of investigation in this paper.We should clarify that the value ofK derived from energy gradient method can only be applied in situations where instabilities could occur. Thus, all the subsequent numerical simulations should satisfy the condition that Ra is larger than Pr3tan6θ•Rac4/(1+Pr)3(1+tan2θ)3, in order to study the cases of instability occurrence.


## 6. Results and Discussions

### 6.1. Results with New Geometries

We will use the same numerical scheme as that described previously to simulate the natural convection in a cavity. The numerical results and analysis are discussed as follows. In order to avoid the influence of “end effect” on numerical results, we choose another computational domain as shown in Figure8. The aspect ratio of the domain is asL(16.2m)×W(2m). The calculated results at Ra=2.1×108 are shown in Figure 8. According to the criterion based on the Rayleigh number, the flow is unstable. Figure 8(a) shows the total pressure contours. It is found from this picture that there are two areas of unstable region in the domain. This phenomenon is explained as follows. As the boundary conditions are defined, the temperature on the top wall is lower than the temperature of the fluid, which leads the fluid with higher temperature in the cavity to move upwards and gather together. Consequently, pressure difference develops gradually due to the movement of the fluid in the cavity, and the irregular phenomena of the total pressure contours indicate that the flow of the fluid is unstable.Calculated results of various parameters: (a) total pressure contours, (b)x-velocity contours, (c) y-velocity contours, (d) temperature contours, and (e) K contours.
(a)
Total pressure
(b)
x-velocity
(c)
y-velocity
(d)
Temperature
(e)
Value ofKFigure8(b) shows the velocity contours along x-coordinate, the straight line m-m splits the cavity into a left region and a right region, and the minus symbol of velocity means particles move along negative direction of x-coordinate. It is noticeable that the velocity is distinctive in the left area, while the velocity in the right area tends to be 0. However, we attribute this phenomenon to inclination of the cavity. For a horizontal plate, the velocity profile is symmetrical in the whole domain. Furthermore, we observe that the negative speed appears mainly close to the top wall. When the heated fluid expands and moves upwards to the top wall, it can be seen from Figure 8(a) that the pressure close to the top wall is much higher than any other area, and this higher pressure would drive the fluid near the top wall to travel to the left side; therefore, the negative velocity is formed. At the same time, we notice that there exists a continuous area above the top wall where the value of x-velocity (u) is large. In [26], we obtain that the thermal boundary layer, containing lots of heat flux, travels along the bottom wall in the process of natural convection; hence, the flow, we observed, is actually the heat flux; furthermore, the movement direction of the heat flux accords well with the value of x-velocity (u).Figure8(c) shows the velocity contours along y-coordinate; the straight line m-m splits the cavity into a left region and a right region. Similar to Figure 8(b), speed is fierce in the left district and is faint in the right district. Here, we get some new discoveries. Firstly, in comparison to Figure 8(b), we can find that the district with high value of y-velocity (v) is in good accordance with the district with high value of x-velocity (u); however, the areas with high value of v are obviously much larger than the areas with the high value of u; this phenomenon confirms that lots of fluid particles move upwards in this cavity, while, at the same time with the fluid moving upwards, the fluid also travels along the horizontal direction. Secondly, we find that the minus velocity along the negative vertical direction still exists in this cavity. For the temperature on the top wall is lower than that in the cavity, which leads to a large density of fluid near the top wall, then the fluid with large density moves downwards due to the buoyancy. Meanwhile, it is known from Figure 8(a) that the pressure near the top wall is much larger than that of the bottom, and the pressure difference would drive the fluid to move downwards. Consequently, the negative speed is produced under the action both of the buoyancy and the pressure difference.Figure8(d) shows the temperature contours in this cavity; the straight line m-m splits the cavity into a left region and a right region. Firstly, it is surveyed that the locations of instabilities concentrate mainly in the left region close to the top wall, while instabilities seldom occur in the right region of this cavity, and this is due to the inclination of the cavity resulting in lots of heat flux concentrated on the top wall in the left region. In addition, we find a unique location of instability above the bottom. It is known from [26] that the movement of thermal boundary layer leads to the accumulation of heat flux above the bottom; with reference to Figure 8(a), we know that there exists an obvious pressure difference at the location where heat flux is accumulated which leads to the fact that the fluid loses its stability and moves upwards. At last, we observe that the locations of instabilities accord well with the locations with high value of velocity. This illustrates that the instabilities result in the dramatic change of velocity.Figure8(e) shows the contours of the value of K; the straight line m-m splits the cavity into a left region and a right region. It is easy to get some similar and new discoveries by comparing Figure 8(d) to Figure 8(e). First, the area with high value of K is mainly concentrated in the left side of this cavity, and this is similar to the previous discoveries. Second, the area with high value of K is in excellent agreement with the locations of instabilities. This result can verify the accurate prediction of energy gradient method. Third, the manifestation of instability is the formation of vorticity, and positions with the formation of vorticities accord well with the locations with high value of K; however, there exists a region with low value of K inside the vorticity, as the blue area shows in the contours of the value of K in Figure 8(e). Fourth, there exists an area with high value of K on the left wall, yet no vorticity is formed. We can observe from the partial enlargement of isotherms in Figure 8(d) that the thermal boundary layer has lost its stability for it travels in a waveform. The reason why it does not form a vorticity is that the thermal boundary travels along the wall, and the heat flux is accumulated above the bottom which leads to the onset of instability. Except for the above discoveries, there exit two problems to be solved: the first one is that we cannot determine the critical value of K which can predict whether instability could occur in natural convection. The other is that we do not comprehend so far why the area with high value of K and the area with the low value of K coexist at the locations of instabilities.


### 6.2. Effect of Flow Time and Discussions

In order to investigate the effect of the flow time, we choose the same geometry and numerical scheme to achieve some numerical simulations. Figures9 and 10 show the contours of temperature and the value of K at 270 s and 300 s, respectively. The basic parameters are showed in the corresponding figures. Comparing Figures 9 and 10 with Figures 8(d) and 8(e), we can observe the following characteristics. First, the regions with high value of K coincide with the locations of instabilities. Second, with the flow time accumulated, the positions and areas of instabilities are increasing. Third, the area with high value of K and the area with the low value of K coexist at the locations of instabilities, and the reason is still unknown till now. Fourth, the locations of instabilities spread right along the top wall with the increase of the flow time. This is because the accumulation of heat flux above the bottom wall is restricted; in other words, the movement of the thermal boundary layer on the top wall is restrained. Thus, lots of heat flux is accumulated on the top wall near the right side which leads to the occurrence of instabilities. Fifth, the fluid in the cavity tends to flow in a turbulent manner. It is suggested in energy gradient method that the amplification of local turbulence phenomenon would lead to that the whole domain loses its stability. In the current case, the natural convection cannot be restricted so that the fluid tends to be turbulence, and this phenomenon is in good agreement with energy gradient method.Calculated results at 270 s: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKCalculated results at 300 s: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the flow time affects numerical results to a considerable degree. Two features are distinctive particularly: the locations of instabilities migrate right on the top wall and the fluid in the whole domain tends to flow in a turbulent manner as the flow time increases.


### 6.3. Effect of Plate Length and Discussions

Figure11 shows the contours of temperature and the value of K, respectively, with a plate length of 24.3 m. The basic parameters are showed in the corresponding figures. Comparing Figures 8(d) and 8(e) with Figures 11(a) and 11(b), we can also make some conclusions. First, the regions with high value of K coincide with the locations of instabilities. Second, the area with high value of K and the area with the low value of K coexist at the locations of instabilities. The above two viewpoints are similar to the former conclusions. Third, the onset time of instability is inversely proportional to the length of the plate; that is, the onset time of instability decreases with the increase of L. Corcione [27] stated that the heat transfer rate from any heated or cooled boundary surface of the enclosure increases as the Rayleigh number increases. This phenomenon is in good accordance with the results in [27]. In all the related cases, Ra is actually proportional to the length of the cavity L. Consequently, when we substitute the larger parameter of L into (2), it is easy to get a shorter time of ts which means that the flow is easier to lose its stability. Fourth, the locations of instabilities migrate right as the plate length increases. Because the migration length of the thermal boundary on the top wall is prolonged, the heat flux transfer would be limited as long as the instabilities set in, which in turn leads to the accumulation of heat flux on the top wall near the right side, and the heat flux loses its stability.Calculated results with length of 24.3 m: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the length of the cavity has an obvious effect on numerical simulation. Two features are distinctive particularly: the onset time of instabilities will get shorter and the locations of instabilities on the top wall will migrate right with the increase of the plate length.


### 6.4. Effect of Inclination Angle and Discussions

Figures12 and 13 show the contours of temperature and the value of K with inclination angle θ=10° and θ=15°, respectively. The basic parameters are showed in the corresponding figures. Kurian et al. [28] studied laminar natural convection inside inclined cylinders of unity aspect ratio and demonstrated that there exits a threshold of inclination angle. As the inclination angle increased from 0 deg to its threshold, convection increased to a maximum. However, when the inclination angle was greater than its critical value, the convection effects and the dimensionless axial temperature gradient decreased relatively small with the increase of inclination angle. Comparing these figures with Figures 8(d) and 8(e), we can obtain some similar and different observations. First, the regions with high value of K still coincide with the locations of instabilities. Second, the area with high value of K and the area with the low value of K coexist at the locations of instabilities. The above two observations are similar to the former conclusions, while the following three conclusions are different. Third, the intensity of flow instability in the left region is much stronger than that in the right region as the inclination angle increases. Upton and Watt [29] made experimental study in an inclined rectangular enclosure, and the results showed that the angle of inclination has a significant effect on the flow and heat transfer in natural convection in an enclosure. Buoyancy in the intrusion layer was found to be the main factor determining the character of these flows. When the angle of the inclined plate increases, more components of buoyancy will be imposed on the top wall which leads to that more heat flux will be accumulated in the left region. Hence, it is easy to observe that the intensity of flow instability is much stronger in the left region. This result is in good accordance with the conclusion stated by Upton and Watt [29]. Fourth, the locations of instabilities would spread left along the top wall with the increase of inclination angle in a certain range. Due to the same reason as stated previously, more heat flux will be accumulated in the left region; thus, it is easy to observe that more instabilities would occur in the left region as the inclination angle increases in a certain range which in turn hinders the movement of the heat flux. Fifth, the locations of instabilities would migrate right once the inclination angle exceeds a certain range. Although the increased component of buoyancy along the top wall can accelerate the movement of heat flux to the left side, the width of the cavity and the higher pressure difference will restrict the accumulation of heat transfer severely in the lower left corner. As a result, more heat flux will be accumulated on the top wall near the right side, and it is easy to observe lots of instabilities on the right side of the top wall. At the same time, we are surprised to find that the fourth and fifth viewpoints accord well with the results of Kurian et al. [28] to some extent.Calculated results with 10 degrees of angle: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKCalculated results with 15 degrees of angle: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the inclination angle of the cavity has a distinctive affection on numerical simulations. Three features are distinctive particularly. First, the intensity of flow instability in the left region is much stronger than that in the right region with the increase of inclination angle. In addition, the locations of instabilities will migrate left with the increase of inclination angle in a certain range. At last, the locations of instabilities will migrate right once the inclination angle exceeds a certain range.


## 6.1. Results with New Geometries

We will use the same numerical scheme as that described previously to simulate the natural convection in a cavity. The numerical results and analysis are discussed as follows. In order to avoid the influence of “end effect” on numerical results, we choose another computational domain as shown in Figure8. The aspect ratio of the domain is asL(16.2m)×W(2m). The calculated results at Ra=2.1×108 are shown in Figure 8. According to the criterion based on the Rayleigh number, the flow is unstable. Figure 8(a) shows the total pressure contours. It is found from this picture that there are two areas of unstable region in the domain. This phenomenon is explained as follows. As the boundary conditions are defined, the temperature on the top wall is lower than the temperature of the fluid, which leads the fluid with higher temperature in the cavity to move upwards and gather together. Consequently, pressure difference develops gradually due to the movement of the fluid in the cavity, and the irregular phenomena of the total pressure contours indicate that the flow of the fluid is unstable.Calculated results of various parameters: (a) total pressure contours, (b)x-velocity contours, (c) y-velocity contours, (d) temperature contours, and (e) K contours.
(a)
Total pressure
(b)
x-velocity
(c)
y-velocity
(d)
Temperature
(e)
Value ofKFigure8(b) shows the velocity contours along x-coordinate, the straight line m-m splits the cavity into a left region and a right region, and the minus symbol of velocity means particles move along negative direction of x-coordinate. It is noticeable that the velocity is distinctive in the left area, while the velocity in the right area tends to be 0. However, we attribute this phenomenon to inclination of the cavity. For a horizontal plate, the velocity profile is symmetrical in the whole domain. Furthermore, we observe that the negative speed appears mainly close to the top wall. When the heated fluid expands and moves upwards to the top wall, it can be seen from Figure 8(a) that the pressure close to the top wall is much higher than any other area, and this higher pressure would drive the fluid near the top wall to travel to the left side; therefore, the negative velocity is formed. At the same time, we notice that there exists a continuous area above the top wall where the value of x-velocity (u) is large. In [26], we obtain that the thermal boundary layer, containing lots of heat flux, travels along the bottom wall in the process of natural convection; hence, the flow, we observed, is actually the heat flux; furthermore, the movement direction of the heat flux accords well with the value of x-velocity (u).Figure8(c) shows the velocity contours along y-coordinate; the straight line m-m splits the cavity into a left region and a right region. Similar to Figure 8(b), speed is fierce in the left district and is faint in the right district. Here, we get some new discoveries. Firstly, in comparison to Figure 8(b), we can find that the district with high value of y-velocity (v) is in good accordance with the district with high value of x-velocity (u); however, the areas with high value of v are obviously much larger than the areas with the high value of u; this phenomenon confirms that lots of fluid particles move upwards in this cavity, while, at the same time with the fluid moving upwards, the fluid also travels along the horizontal direction. Secondly, we find that the minus velocity along the negative vertical direction still exists in this cavity. For the temperature on the top wall is lower than that in the cavity, which leads to a large density of fluid near the top wall, then the fluid with large density moves downwards due to the buoyancy. Meanwhile, it is known from Figure 8(a) that the pressure near the top wall is much larger than that of the bottom, and the pressure difference would drive the fluid to move downwards. Consequently, the negative speed is produced under the action both of the buoyancy and the pressure difference.Figure8(d) shows the temperature contours in this cavity; the straight line m-m splits the cavity into a left region and a right region. Firstly, it is surveyed that the locations of instabilities concentrate mainly in the left region close to the top wall, while instabilities seldom occur in the right region of this cavity, and this is due to the inclination of the cavity resulting in lots of heat flux concentrated on the top wall in the left region. In addition, we find a unique location of instability above the bottom. It is known from [26] that the movement of thermal boundary layer leads to the accumulation of heat flux above the bottom; with reference to Figure 8(a), we know that there exists an obvious pressure difference at the location where heat flux is accumulated which leads to the fact that the fluid loses its stability and moves upwards. At last, we observe that the locations of instabilities accord well with the locations with high value of velocity. This illustrates that the instabilities result in the dramatic change of velocity.Figure8(e) shows the contours of the value of K; the straight line m-m splits the cavity into a left region and a right region. It is easy to get some similar and new discoveries by comparing Figure 8(d) to Figure 8(e). First, the area with high value of K is mainly concentrated in the left side of this cavity, and this is similar to the previous discoveries. Second, the area with high value of K is in excellent agreement with the locations of instabilities. This result can verify the accurate prediction of energy gradient method. Third, the manifestation of instability is the formation of vorticity, and positions with the formation of vorticities accord well with the locations with high value of K; however, there exists a region with low value of K inside the vorticity, as the blue area shows in the contours of the value of K in Figure 8(e). Fourth, there exists an area with high value of K on the left wall, yet no vorticity is formed. We can observe from the partial enlargement of isotherms in Figure 8(d) that the thermal boundary layer has lost its stability for it travels in a waveform. The reason why it does not form a vorticity is that the thermal boundary travels along the wall, and the heat flux is accumulated above the bottom which leads to the onset of instability. Except for the above discoveries, there exit two problems to be solved: the first one is that we cannot determine the critical value of K which can predict whether instability could occur in natural convection. The other is that we do not comprehend so far why the area with high value of K and the area with the low value of K coexist at the locations of instabilities.


## 6.2. Effect of Flow Time and Discussions

In order to investigate the effect of the flow time, we choose the same geometry and numerical scheme to achieve some numerical simulations. Figures9 and 10 show the contours of temperature and the value of K at 270 s and 300 s, respectively. The basic parameters are showed in the corresponding figures. Comparing Figures 9 and 10 with Figures 8(d) and 8(e), we can observe the following characteristics. First, the regions with high value of K coincide with the locations of instabilities. Second, with the flow time accumulated, the positions and areas of instabilities are increasing. Third, the area with high value of K and the area with the low value of K coexist at the locations of instabilities, and the reason is still unknown till now. Fourth, the locations of instabilities spread right along the top wall with the increase of the flow time. This is because the accumulation of heat flux above the bottom wall is restricted; in other words, the movement of the thermal boundary layer on the top wall is restrained. Thus, lots of heat flux is accumulated on the top wall near the right side which leads to the occurrence of instabilities. Fifth, the fluid in the cavity tends to flow in a turbulent manner. It is suggested in energy gradient method that the amplification of local turbulence phenomenon would lead to that the whole domain loses its stability. In the current case, the natural convection cannot be restricted so that the fluid tends to be turbulence, and this phenomenon is in good agreement with energy gradient method.Calculated results at 270 s: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKCalculated results at 300 s: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the flow time affects numerical results to a considerable degree. Two features are distinctive particularly: the locations of instabilities migrate right on the top wall and the fluid in the whole domain tends to flow in a turbulent manner as the flow time increases.


## 6.3. Effect of Plate Length and Discussions

Figure11 shows the contours of temperature and the value of K, respectively, with a plate length of 24.3 m. The basic parameters are showed in the corresponding figures. Comparing Figures 8(d) and 8(e) with Figures 11(a) and 11(b), we can also make some conclusions. First, the regions with high value of K coincide with the locations of instabilities. Second, the area with high value of K and the area with the low value of K coexist at the locations of instabilities. The above two viewpoints are similar to the former conclusions. Third, the onset time of instability is inversely proportional to the length of the plate; that is, the onset time of instability decreases with the increase of L. Corcione [27] stated that the heat transfer rate from any heated or cooled boundary surface of the enclosure increases as the Rayleigh number increases. This phenomenon is in good accordance with the results in [27]. In all the related cases, Ra is actually proportional to the length of the cavity L. Consequently, when we substitute the larger parameter of L into (2), it is easy to get a shorter time of ts which means that the flow is easier to lose its stability. Fourth, the locations of instabilities migrate right as the plate length increases. Because the migration length of the thermal boundary on the top wall is prolonged, the heat flux transfer would be limited as long as the instabilities set in, which in turn leads to the accumulation of heat flux on the top wall near the right side, and the heat flux loses its stability.Calculated results with length of 24.3 m: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the length of the cavity has an obvious effect on numerical simulation. Two features are distinctive particularly: the onset time of instabilities will get shorter and the locations of instabilities on the top wall will migrate right with the increase of the plate length.


## 6.4. Effect of Inclination Angle and Discussions

Figures12 and 13 show the contours of temperature and the value of K with inclination angle θ=10° and θ=15°, respectively. The basic parameters are showed in the corresponding figures. Kurian et al. [28] studied laminar natural convection inside inclined cylinders of unity aspect ratio and demonstrated that there exits a threshold of inclination angle. As the inclination angle increased from 0 deg to its threshold, convection increased to a maximum. However, when the inclination angle was greater than its critical value, the convection effects and the dimensionless axial temperature gradient decreased relatively small with the increase of inclination angle. Comparing these figures with Figures 8(d) and 8(e), we can obtain some similar and different observations. First, the regions with high value of K still coincide with the locations of instabilities. Second, the area with high value of K and the area with the low value of K coexist at the locations of instabilities. The above two observations are similar to the former conclusions, while the following three conclusions are different. Third, the intensity of flow instability in the left region is much stronger than that in the right region as the inclination angle increases. Upton and Watt [29] made experimental study in an inclined rectangular enclosure, and the results showed that the angle of inclination has a significant effect on the flow and heat transfer in natural convection in an enclosure. Buoyancy in the intrusion layer was found to be the main factor determining the character of these flows. When the angle of the inclined plate increases, more components of buoyancy will be imposed on the top wall which leads to that more heat flux will be accumulated in the left region. Hence, it is easy to observe that the intensity of flow instability is much stronger in the left region. This result is in good accordance with the conclusion stated by Upton and Watt [29]. Fourth, the locations of instabilities would spread left along the top wall with the increase of inclination angle in a certain range. Due to the same reason as stated previously, more heat flux will be accumulated in the left region; thus, it is easy to observe that more instabilities would occur in the left region as the inclination angle increases in a certain range which in turn hinders the movement of the heat flux. Fifth, the locations of instabilities would migrate right once the inclination angle exceeds a certain range. Although the increased component of buoyancy along the top wall can accelerate the movement of heat flux to the left side, the width of the cavity and the higher pressure difference will restrict the accumulation of heat transfer severely in the lower left corner. As a result, more heat flux will be accumulated on the top wall near the right side, and it is easy to observe lots of instabilities on the right side of the top wall. At the same time, we are surprised to find that the fourth and fifth viewpoints accord well with the results of Kurian et al. [28] to some extent.Calculated results with 10 degrees of angle: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKCalculated results with 15 degrees of angle: (a) temperature contours and (b)K contours.
(a)
Temperature
(b)
Value ofKIn summary, the inclination angle of the cavity has a distinctive affection on numerical simulations. Three features are distinctive particularly. First, the intensity of flow instability in the left region is much stronger than that in the right region with the increase of inclination angle. In addition, the locations of instabilities will migrate left with the increase of inclination angle in a certain range. At last, the locations of instabilities will migrate right once the inclination angle exceeds a certain range.


## 7. Conclusions

Numerical simulations on natural convection in an inclined cavity have been carried out using the unsteady Navier-Stokes equations for various parameters such as the domain length and the inclined angle. The main conclusions are as follows.(1)
The energy gradient method is successfully employed to study the instability of the thermal boundary layer. It is found that instability occurs firstly at the position where the energy gradient function gets its maximumKmax. The regions with high value of K coincide with the locations of instabilities. These observations accord well with the energy gradient theory.(2)
The fluid within the whole domain tends to flow in a turbulent state with the increase of flow time if Ra is larger than its critical value.(3)
With the increase of plate length, the onset time of instability decreases and the locations of instabilities migrate right along the top wall.(4)
The intensity of instabilities in the left region of the studied inclined cavity is much stronger than that in the right region. When the inclined angle increases in a certain range, the locations of instabilities along the top wall migrate left.(5)
Once the inclined angle exceeds a certain range, the positions of instabilities would move right since the accumulation of the heat flux in the lower left corner of the cavity is restrained.



---

*Source: 198695-2013-03-07.xml*