Analysis of Mixed Convection Flow Characteristics in Square Cavity with Uniformly Heated Bottom Wall by Finite Volume Method

Analysis of Mixed Convection Flow Characteristics in Square Cavity with Uniformly Heated Bottom Wall by Finite Volume Method

Robins Aikkara

Dept. of Mechanical Engineering LBS College of Engineering Kasaragod, India

Aboobacker Kadengal

Dept. of Mechanical Engineering LBS College of Engineering Kasaragod, India

Abstract Mixed convection from a uniform heat source on the bottom of a square cavity is studied numerically. Two- dimensional forms of non-dimensional Navier-Stokes equations are solved by using control volume based finite volume technique. Three typical values of the Reynolds numbers are chosen as Re = 1, 10, and 100 and steady, laminar results are obtained in the values of Richardson number as Ri = 0, 1 and 10 and the values of Prandtl numbers as Pr = 0.1, 0.71, 1 and 10. The parametric studies for a wide range of governing parameters show consistent performance of the present numerical approach to obtain as stream functions and temperature profiles. Heat transfer rates at the heated walls are presented based on the value of Re and Pr. The computational results indicate that the heat transfer is strongly affected by Reynolds number and Richardson number. In the present investigation, bottom wall is uniformly heated while the two vertical walls are maintained at constant cold temperature and the top wall is well insulated. A complete study on the effect of Ri shows that the strength of circulation increases with the increase in the value of Ri irrespective of Re and Pr. As the value of Ri increases, there occurs a transition from conduction to convection dominated flow at Ri =1. A detailed analysis of flow pattern shows that the natural or forced convection is based on both the parameters Ri and Pr.

Keywords Mixed convection; Square cavity; Lid driven cavity; Uniform heating; Reynolds number; Richardson number and Prandtl number.

1. INTRODUCTION

   In a mixed convection, both natural convection and forced convection participate in the heat transfer process. The bulk fluid flow direction can be any of the three possible directions in a horizontal channel, forward, backward or upward. The forced flow can be in the same direction as the flow created by natural convection, and this flow condition is called assisting mixed convection. Whereas, for the other case, forced flow direction is in an opposing direction to the flow that is created by buoyancy, and this flow condition is referred to as opposing mixed convection.

   Since the early work by Burggraf [1], the lid-driven cavity flow is considered as the classical test problem for the assessment of numerical methods and the validation of Navier-Stokes codes. Highly-accurate solutions for the lid- driven cavity flow are computed by a Chebyshev collocation method is done by O. Botella [2]. Accuracy of the solution is

   achieved by using a subtraction method of the leading terms of the asymptotic expansion of the solution of the Navier- Stokes equations in the vicinity of the corners, where the velocity is discontinuous. Critical comparison with former numerical experiments confirms the high-accuracy of the method, and extensive results for the flow at Reynolds number Re = 1000 are presented. The Charles-Henri Bruneau [3], numerically simulate of the 2D lid-driven cavity flow are performed for a wide range of Reynolds numbers. Accurate benchmark results are provided for steady solutions as well as for periodic solutions around the critical Reynolds number.

   Moallemi and Jang [4] analyzed the effects of Prandtl number (Pr) on laminar mixed convection heat transfer in a lid driven cavity. They performed the numerical simulations for two-dimensional laminar flow (100 Re 2200) and studied the effects of small to moderate Prandtl numbers (0.01 Pr 50) on the flow and heat transfer characteristics in a square cavity for various values of Richardson number (Ri). The temperature and flow fields in the cavity show the strong influence of Prandtl number, Pr. The local and average Nusselt numbers are also reported for various values of Re, Ri and Pr. Mixed convection heat transfer in a lid driven cavity was also investigated by Prasad and Koseff [5]. They performed a series of experiments in a cavity filled with water and measured the heat flux at different locations over the hot cavity floor for a range of Re and Gr. Their results indicate that the overall (i.e. area averaged) heat transfer rate is a very weak function of Gr for 2200 Re 12000.

   In the present investigation, bottom wall is uniformly heated while the two vertical walls are maintained at constant cold temperature and the top wall is well insulated. A complete study on the effect of Gr shows that the strength of circulation increases with the increase in the value of Ri irrespective of Re and Pr. As the value of Ri increases, there occurs a transition from conduction to convection dominated flow. A detailed analysis of flow pattern shows that the natural or forced convection is based on both the parameters Ri and Pr.

2. MATHEMATICAL FORMULATION

   A two-dimensional square cavity is considered for the present study with the physical dimension as shown in Fig.1. The bottom wall of the cavity is maintained at a uniform temperature and the upper wall is insulated. The two vertical walls are maintained at cold temperature. It may be noted that the bottom wall is maintained at a higher temperature to induce buoyancy effect. The top wall is assumed to slide from left to right with a constant speed U0. The flow is assumed to be laminar and the fluid properties are assumed to be constant except for the density variation which is modeled according to Boussinesq approximation while viscous dissipation effects are considered to be negligible. The viscous incompressible flow and the temperature distribution inside the cavity are governed by the Navier Stokes and the energy equations, respectively. The aim of the current work is to investigate the steady state solutions and hence, we have considered the time independent differential governing equations. Similar procedure was also followed in the recent work on mixed convection. A number of earlier works was based on steady state solutions which were obtained via steady mathematical model. The governing equations are non-dimensionalized to yield

   (1)

   (2)

   (3)

   (4)

   Fig.1. Schematic diagram of a physical system

   The transformed boundary conditions are: U (X, 1) = 1,

   U (X, 0) = U (0, Y) = U (1, Y) = 0

   V (X, 0) = V (X, 1) = V (0, Y) = V (1, Y) = 0 (X, 0) = 1

   (0, Y) = (1, Y) = 0

   The dimensionless variables and parameters are defined as follows:

   , ,

   , ,

   , ,

   Here x and y are the distances measured along the horizontal and vertical directions, respectively; u and v are the velocity components in x and y directions, respectively; T denotes the temperature; p is the pressure and is the density; Th and Tc are the temperature at the hot and cold walls, respectively; L is the length of the side of the square cavity; X and Y are dimensionless coordinates varying along horizontal and vertical directions, respectively; U0 is the velocity of the upper wall; U and V are dimensionless velocity components in the X and Y directions, respectively; is the dimensionless temperature; P is the dimensionless pressure; Gr, Re and Pr are Grashof, Reynolds and Prandtl number, respectively.

3. SOLUTION METHODOLOGY

   The momentum and energy balance equations [Eqs.

   (2) (4)] are the combinations of a system of equations which have been solved using the finite volume method. The continuiy equation [Eq. (1)] has been used as a constraint due to mass conservation and this constraint may be used to obtain the pressure distribution. In order to solve Eqs. (2) (3), we use the finite volume discretisation procedure, if the nonlinearity in the momentum equations appears to be a difficulty in iteration procedure. Starting with a guessed velocity field, we could iteratively solve the momentum equations to arrive at the converged solution for the velocity components. The real difficulty in the calculation of the velocity field lies in the unknown pressure field. The pressure gradient forms a part of the source term for a momentum equation. For a given pressure field, there is no particular difficulty in solving the momentum equation. The pressure field is indirectly specified via the continuity equation [Eq. (1)]. When the correct pressure field is substituted into the momentum equations, the resulting velocity field satisfies the continuity equation.

   The discretised form of X-momentum [Eq. (2)] is written as

   (5)

   where Ae = Y

   The discretised form of Y-momentum [Eq. (3)] is written as

   (6)

   where An = X.

   The momentum equations can be solved only when the pressure field is given. Unless the correct pressure field is employed the resulting velocity field will not satisfy the continuity equations. Such an imperfect velocity field based on a guessed pressure field P* will be denoted by U*, V*. This

   starred velocity field will result from the solution of the following discretisation equations:

   (7)

   (8)

   In these equations [Eq. (7), (8)], the velocity components and pressure have been given the superscript *.

   The guessed pressure P* and resulting starred velocity field will progressively get closer to satisfying the continuity equation. Let us propose that the correct pressure P is obtained from

   (9)

   where will be called the pressure correction. The corresponding velocity , can be introduced in a similar manner:

   , (10)

   The equation (12) will be used to correct the momentum equations:

   (11)

   (12)

   The computational procedure is similar to the one described by Baliga and Patankar (1983), and Gresho et al. (1984). The resulting system of the coupled equations (1-4) with the associated boundary conditions have been solved numerically using control volume based finite volume method. The computational domain consists of 50×50 main grid points which correspond to 50×40 U and V staggered grid points. The control volume based finite volume method provides the smooth solutions at the interior domain including the corner regions. To ensure the convergence of the numerical solution to the exact solution, the grid sizes have been optimized and the results presented here are independent of grid sizes.

   Grid refinement tests have been performed for the case Re = 100 and Gr = 1000 using eight uniform grids 15×15, 25×25, 40×40, 50×50, 60×60, 75×75, 90×90 and

   100×100. Results show that when we change the mesh size from a grid of 25×25to a grid of 50× 50, the maximum value of stream function contour (max) and the maximum temperature contour (max) undergoes an increase of only 0.5% and 0.25%, respectively; then, because of calculation cost, the 50×50 grid is retained.

   The computer code has been validated with the solutions are available in the literatures. There are some possibilities of validating the numerical code. One possibility is to compare the numerical results obtained by our code with benchmarks available in the literature according to different works. Another option is to simulate a similar problem investigated by other authors with well accepted available

   results. The code validated by using lid-driven cavity benchmark problem by O. Botella and R. Peyret (1998).

   The computational model is validated for mixed convection heat transfer by comparing the results of correlation on mixed convection in lid-driven cavity with uniform heat flux in bottom wall performed by Tanmay Basak and I. Pop (2009). The right, top and left walls are insulated.

4. RESULT AND DISCUSSION

   The computational domain consists of 50×50 main grid points which correspond to 50×40 U and V staggered grid points. Numerical solutions are obtained for various values of Ri = 0 – 10, Pr = 0.0110 and Re = 1103 with uniform heating of the bottom wall where the two vertical walls are cooled and the top wall is well insulated with a horizontal velocity, U=1. The jump discontinuity in Dirichlet type of wall boundary conditions at the corner point (see Fig.1) corresponds to computational singularity. To ensure the convergence of the numerical solution to the exact solution, the grid sizes have been optimized and the results presented here are independent of grid sizes

   1. Effect of Richardson number

      Stream function, Temperature,

      (a)

      (b)

      (c)

      Fig.2. Stream function and temperature contours with Pr =0.1, Re =1:

      1. Ri =0, (b) Ri =1, and (c) Ri =10

         Fig. 2- 5 illustrate the stream function and isothermal contours with the constant temperature heated bottom wall for Ri = 0 10, Pr = 0.1 10 and Re = 1. Fig.2

         1. shows that the effect of lid-driven flow predominates the forced convection for Gr = 0, Pr = 0.1 and Re = 1. The value of Ri = 0, the effect of buoyancy gradually weaker as compared to the lid driven force. All amount of fluid is pulled up towards the left corner due to drag force created by the motion of the upper lid. The clockwise circulation is found to be more as compared to anticlockwise circulation. The secondary circulation occupies a very small portion of the bottom corners and primary circulation occupies the major portion of the cavity. As a result, the stream function contours near to the upper lid are not exactly oval shape is observed.

            The value of Reynolds number is very low and heat transfer is conduction dominated within the cavity. Due to dominant conduction mode, all the isotherms are smooth symmetric curve that span the entire cavity. Fig.2 (b) shows the contour plots for Ri = 1. That means strength of buoyancy same as compared to the lid driven force (Re = 1 and Gr = 1). The stream lines at the centre of the cavity are mostly elliptic in shape and far away from the centre are changes it shapes.

            Stream function, Temperature,

            (a)

            (b)

            (c)

            Fig.3. Stream function and temperature contours with Pr =0.7, Re =1:

            1. Ri =0, (b) Ri =1, and (c) Ri =10

         Fig.3 (a)-(c) show the stream function and temperature contour for various Ri with Pr = 0.7 and Re =1.

         Fig.3 (a) displays the effect of lid driven flow inside the cavity for Ri = 0, similar to Fig.2 (a), the strength of clockwise circulation stronger inside the cavity. The value of Ri = 0, strength of buoyancy is negligible in this diagram. It is interesting to observe that the temperature distribution are coupled with stream function for Pr = 0.7 and hence the isotherms gradually small tend to be asymmetric. Fig.3 (b) shows the contour plots for Ri = 1, similar to Fig.2 (b) the effect of lid driven flow and buoyancy effect is same. Fig.3

         1. illustrates contours for Ri = 10 the stream function circulations are observed to be symmetric.

            The effect of lid-driven flow totally vanishes for higher value of Ri as the buoyancy becomes the dominant force due to increase in the value of Gr. The temperature gradients near both the bottom and side wall tend to be significant leading to the development of thermal boundary layer.

            Fig.4 (a) – (c) display the effect of Ri varying within

            0 to 10 and Pr = 1. It is interesting to observe that the isotherms with 0.5 are coupled with stream functions and hence the isotherms gradually tend to be asymmetric. The stream function contours for Fig.4 (a) – (c) are similar to Fig.3 (a) – (c).

            Stream function, Temperature,

            (a)

            (b)

            (c)

            Fig.4. Stream function and temperature contours wit Pr =1, Re =1: (a) Ri =0, (b) Ri =1, and (c) Ri =10

            Figs. 5 (a) (c) illustrate the stream function and temperature contour for various Ri with Pr =10 and Re =1. Fig.5 (a) display the effect of Ri = 0, that means the effect of lid-driven flow dominates over the strength of buoyancy. Due to the low Reynolds number flow stream function circulation is symmetric about the vertical axis. It is observe that the temperature distribution are coupled with stream functions for Pr = 10 and hence the isotherms gradually tend to be asymmetric. Fig.5 (b) shows the contour plots for Ri = 1, similar to Fig.4 (c), the effect of lid driven flow gradually diminishes with the increase of Gr. Then the temperature contours start getting shifted towards left side walls. Fig.5 (c) illustrates contours for Ri = 10, here stream function circulation is observed to be symmetric. The effect of lid- driven flow totally vanishes for higher value of Ri as the buoyancy becomes the dominant force due to increase in the value of Gr.

            Stream function, Temperature,

            (a)

            (b)

            (c)

            Fig.5. Stream function and temperature contours with Pr =10, Re =1:

            1. Ri =0, (b) Ri =1, and (c) Ri =10

   2. Effect of Reynolds number

      Figs.6 and 7 display the stream function and isothermal contours for Re = 1, 10, 100 and 1000 with Ri = 1 corresponding to Pr = 0.7 (see Fig. 6) and 10 (see Fig. 7) with constant temperature heating of bottom wall. The increase in Re enhances effect of forced convection and suppresses the

      effect of natural convection. This can also be explained based on Richardson number Ri.

      Stream function, Temperature,

      (a)

      Fig.6. Stream function and temperature contours with Pr = 0.7, Ri =1: (a) Re =1, (b) Re =10, (c) Re =102, (d) Re =103

      Fig.6 (a)(d) distribution for Pr =0.7 and Ri =1 with various values of Re. Fig. 6(a) displays the distributions for Re =1 and Gr =1 which denotes a combined effect of buoyancy and lid-driven force inside the cavity. The clockwise circulation is found to be more as compared to anticlockwise circulation. The secondary circulation occupies a very small portion of the bottom corners and primary circulation occupies the major portion of the cavity. As a result, the stream function contours near to the upper lid are not exactly oval shape is observed. The value of Reynolds number is very low and heat transfer is conduction dominated

      within the cavity. Due to dominant conduction mode, all the isotherms are smooth symmetric curve that span the entire cavity. All the isotherms span the entire enclosure and they are symmetric with respect to vertical centre line. The conduction dominance may also be confirmed by the fact that the stream function and temperature contours are not coupled. Fig.6 (b) shows the effect for Re = 10 and Gr =100. As Re increases corresponding to Gr increases then the value of Ri as constant, the effect of buoyancy and lid driven forces are dominated. The secondary circulation occupies a very small portion of the bottom corners and the primary circulation occupies the major portion of the cavity. It is observe that the temperature distribution are coupled with stream function for Pr =0.7 and hence the isotherms gradually tend to be asymmetric. Fig.6 (c) isotherms gradually tend to be asymmetric.

      Stream function, Temperature,

      (a)

      (b)

      (c)

      (d)

      Fig.7. Stream function and temperature contours with Pr = 10, Ri =1: (a) Re =1, (b) Re =10, (c) Re =102, (d) Re =103

      Fig.6 (c) displays the distributions for Re = 100 and Gr = 104 which denotes a combined effect of buoyancy and lid-driven force inside the cavity. Due to the increase of Re, stream function circulation near the central region start getting shifted towards the right side wall. The circulation near the central region becomes stronger and consequently, the temperature contours with = 0.1-0.3 start getting shifted towards the side walls (see Fig.6 (c)). The presence of significant convection is also exhibited in other temperature contour lines which start getting deformed and lifted towards the top plate. Fig.6 (d) shows the effect for Re = 103 and Gr = 106, the effect of buoyancy as well as lid driven force also affect the circulation.

      The secondary circulation occupies a small portion of the bottom right corner and the primary circulation occupies the major portion of the cavity. It is interesting to observe that the strength of secondary circulation increases and strength of primary circulation decreases at the centre of the cavity. The temperature gradient near the bottom wall increases and the isotherms start shifting from the right cold wall. As the value of Re increases, it is observed that a large region near the right half becomes isothermally cooled and the effect of heating will be confined only near the bottom and left walls of the cavity.

      Fig.7 shows the distributions for higher Prandtl number fluid (Pr = 10) with Re = 1, 10, 102, 103 and Ri = 1. Fig.7 (a) shows contour plots for Ri = 1. That means strength of buoyancy same as compared to the lid driven force (Re = 1 and Gr = 1). The stream lines at the centre of the cavity are mostly elliptic in shape and far away from the centre are changes it shapes. . It is interesting to observe that the temperature distribution are coupled with stream function for Pr = 10 and hence the isotherms gradually tend to be asymmetric. The Fig.7 (b) displays the effect of lid driven flow inside the cavity for Re =10 and Gr =102, it stream function is similar to Fig.6 (b). The temperature contours with = 0.1-0.3 start getting shifted towards the side walls (see Fig. 7 (b)). The presence of significant convection is also exhibited in other temperature contour lines which start getting deformed and lifted towards the top plate. Fig.7 (c) shows the contour plots for Re = 102 and Gr =104. Due to the increase of Re, stream function circulation near the central region start getting shifted towards the right side wall. It is observed that a large region near the right half become isothermally cooled and the effect of heating will be confined only near the bottom and left walls of the cavity forming a strong thermal boundary layer attached to the bottom wall occupying only 30% of the cavity. Fig.7 (d) displays the effect of lid driven flow inside the cavity for Re =103 and Gr

      =106, it stream function is similar to Fig.6 (d). The anti- clockwise circulation occupies a small portion of the bottom right corner and the clockwise circulation occupies the major portion of the cavity. It is interesting to observe that the strength of anti-clockwise circulation increases and strength of clockwise circulation decreases at the centre of the cavity. A lager cold region near the right wall is observed as seen in Fig.7 (d) and the region of uniform heating of bottom wall.

   3. Effect of Prandtl number

   Stream function, Temperature,

   (a)

   (b)

   (c)

   (d)

   Fig.8. Stream function and temperature contours with

   Re =10, Ri =1 :(a) Pr = 0.1, (b) Pr = 0.7, (c) Pr = 1, and (d) Pr =10

   The effect of Pr for representative higher values of Re and Ri has been illustrated in Fig. 8 9. Although Ri is fixed in each case, the figures show that the flow inside the cavity is a function of Pr. It is interesting to observe that either conduction, or natural convection or forced convection is dominant for various values of Pr with Ri = 1 in presence of uniform heating of bottom wall (see Fig.8 (a)(c)). The isotherms illustrate conduction dominance at Pr = 0.1 with the maximum value of stream function being 0.18. It is interesting to observe that the circulation cells are dominant which signify that the natural convection, however weak, is gradually developed and the conduction dominant heating effect is observed due to less intense circulations (Fig.8 (a)).

   The flow circulations are qualitatively similar for Pr = 0.7. The mximum value of stream function for the circulation cell is 0.2 and the size of the circulation cells is larger than that with Pr = 0.1. Due to enhanced circulation and thermal mixing at the right half, the isotherms are pushed towards the left-side wall. The dominant natural convection is attributed to the asymmetric isotherms as seen in Fig.8 (b). Fig.8 (c) displays the distributions for Pr = 1, it is similar to Fig.8 (b).

   Stream function, Temperature,

   (a)

   (b)

   (c)

   (d)

   Fig.9. Stream function and temperature contours with

   Re =100, Ri =1 :(a) Pr = 0.1, (b) Pr = 0.7, (c) Pr = 1, and (d) Pr =10

   It is observe that temperature distributions are coupled with stream function for Pr = 1 and hence the isotherms gradually tend to be asymmetric. A thin boundary layer is found to occur along the top portion of the right wall whereas the isotherms are pushed near the left wall. Thus

   strong thermal boundary layers are also formed along the left and bottom walls.

   Fig.9 (a)(c) illustrate the stream function and temperature contours for Re=102 and Gr = 104 with various Pr in presence of uniform heating of the bottom wall. It may be noted that, Ri =1 for all the cases and it is interesting to observe that primary circulation cells almost occupy the entire cavity. It is also observed that secondary circulation cells appear near the corner of the bottom wall. Fig.9 (a) shows that the temperature distributions are weakly coupled with stream function for Pr = 0.1 and hence the isotherms tend to be symmetric resulting in conduction dominated flow. As Pr increases to 0.7 and 1 the isotherms are clustered towards the bottom and left walls and convection plays a dominant role in heat transfer. For further increase of Pr to 10, the compression of isotherms is more prominent and a strong thermal boundary layer develops near the bottom and left wall (see Fig.9 (d)). Due to the compression of isotherms towards the left wall, the heat transfer rates through the bottom and left walls are more as compared to the right wall.

5. CONCLUSIONS

A numerical investigation on mixed convection in a square cavity with various boundary conditions was carried out using a finite volume method. The prime objective of the investigation is to study the effect of uniform heating of the bottom, on the flow and heat transfer characteristics due to mixed convection in cavity. It is evident from Figs. 2 – 5 for fixed Re and Pr, the strength of circulation increases with the increase in Ri. As Ri increases, the effect of buoyancy increases leading to an increase in the strength of circulation. Due to increase in circulation strength, the isotherms are stretched along the side walls and heat is transferred mostly by convection for higher value of Pr. The effect of Re has also been studied in the present investigation for fixed value of Pr and Ri. It is observed that the effect of natural convection decreases and forced convection increases with the increase of Re. It has also been observed that for higher value of Pr, the effect of heating is more pronounced near the bottom and left walls as the formation of thermal boundary layers is restricted near the bottom and left wall for uniform heating cases. The heat transfer rate is very high at the edges of the bottom wall and it decreases at the center of the cavity.

Laminar convection in a two-dimensional, horizontally driven rectangular enclosure with a prescribed constant temperature heat source mounted on the bottom wall is simulated numerically in this work. Mixed convection arises as the buoyancy-induced cold flow from the source interacts with an externally induced cold air flow. The effects of different ventilation orientations are also described to figure out the best cooling performance. The heat transfer results explain the importance of the non-dimensional parameters like Reynolds number and Richardson number in the natural and mixed convection regime. The effects of these parameters on the flow fields are also investigated. The governing parameter affecting heat transfer is the Richardson number. For Ri > 1, the heat transfer is dominated by natural convection. When Ri < 1, the flow and heat transfer are dominated by forced convection. The mixed regime is obtained when Ri = 1.

ACKNOWLEDGMENT

The authors wish to acknowledge Department of Mechanical Engineering, LBS College of Engineering, Kasaragod, India, for support and technical help throughout this work.

REFERENCES

1. Burggraf, O.R., Analytical and numerical studies of the structure of steady separated flows, Journal of Fluid Mechanics, 1966, 24, 113-151

2. O. Botella and R. Peyret., Benchmark spectral results on the lid-driven cavity flow, Computers & Fluids Vol. 27, No. 4, pp. 421-433, 1998

3. Charles-Henri Bruneau and Mazen Saad., The 2D lid-driven cavity problem revisited, Computers & Fluids 35 (2006) 326348

4. M.K. Moallemi, K.S. Jang, Prandtl number effects on laminar mixed convection heat transfer in a lid-driven cavity, Int. J. Heat Mass Transfer 35 (1992) 18811892.

5. A.K. Prasad, J.R. Koseff, Combined forced and Natural convection heat transfer in a deep lid-driven cavity flow, Int. J. Heat Fluid Flow 17 (1996) 460467.

6. H.F. Oztop, I. Dagtekin, Mixed convection in two sided lid-driven differentially heated square cavity, Int. J. Heat Mass Transfer 47 (2004) 17611769.

7. E. Papanicolaou, Y. Jaluria, Mixed convection from an isolated heat source in a rectangular enclosure, Numer. Heat Transfer 18 (1990) 427461.

8. T.H. Hsu, P.T. Hsu, S.P. How, Mixed convection in a partially divided rectangular enclosure, Numer. Heat Transfer A 31 (1997) 655683.

9. T.H. Hsu, S.G. Wang, Mixed convection in a rectangular enclosure with discrete heat sources, Numer. Heat Transfer A 38 (2000) 627652.