Numerical Investigation of Free Convection Heat Transfer from Two-Dimensional Rectangular Enclosure with Discrete Isothermal Heating from Bottom Side

Numerical Investigation of Free Convection Heat Transfer from Two-Dimensional Rectangular Enclosure with Discrete Isothermal Heating from Bottom Side

Mohammed Ali Mahmood Hussein Ameer Resan Kalash Ibtisam A. Al-BeldaweeLaith Jaafer Habeeb 

Refrigeration and Air Conditioning Engineering Department, Al-Rafidain University Collage, Baghdad 10, Iraq

Training and Workshop Center, University of Technology, Baghdad 10, Iraq

Electromechanical Engineering Department, University of Technology, Baghdad 10, Iraq

Corresponding Author Email:
11 March 2019
22 November 2019
5 December 2019
Available online: 
26 December 2019
| Citation

© 2019 IIETA. This article is published by IIETA and is licensed under the CC BY 4.0 license (



The thermal energy exchange by means of free convection between two sources of heat located in an enclosure with a rectangular shape on the bottom side with sawtooth periodic (repetition) upper wall has been studied numerically. The length of the two isothermal sources was about 20% of the total bottom surface length with a constant power. The impacts of the sawtooth waves and the separation of both heat sources in the heat exchange were analyzed during the numerical simulation. The working fluid was air with Prandtl number (0.72) and Rayleigh numbers that varied from 103 to 105. The finite volume method with SIMPLE algorithm was used to solve the governing partial differential equations (continuity, Navier Stokes and energy equations) and to analyze the flow and temperature fields inside the enclosure. The simulated data revealed that an increase in the distance between two heat sources results in a raise of the mean Nusselt Number for all the ranges of Rayleigh Numbers studied, while an increase in the sawtooth periodic results in a diminishing of the mean Nusselt Number for all the Rayleigh Numbers analyzed.


natural convection, rectangular enclosure, discrete heat sources, finite volume method

1. Introduction

The enclosures Natural convection with differentially heated from the sides and containing solid obstacles or porous media has involved an excessive deal of attention because of its fundamental nature and the wide range of applications in, for example, air conditioning indoor, electronics cooling, food processing, and steel and iron making [1]. Through The issue of heat, the transfer in natural convection is very important because of its many applications in many engineering fields. However, even though the fluid is not subjected to a forced motion when its kind of convection occurs, streams are produced because of the buoyancy effects that are stimulated in the free convection due to the density and temperature gradients. Since there are no external devices to force the motion of a fluid, free convection is a cheap and effective way of exchanging thermal energy, and that is why this study was focused on analyzing the heat transferred in enclosures. The surface of these enclosures was corrugated to improve the heat transfer in many engineering applications, such as the cooling of electronic devices and the design of heat exchangers, as well as in wavy solar collectors. Hasan et al. [2] developed a numerical simulation that allowed determining the free convective heat exchange in a 2D square enclosure with a wavy top surface. The study was carried out for several corrugation amplitudes and aspect ratios to evaluate the hydraulic and thermal behavior. The upper wall was heated under a constant heat flux and the side walls were maintained at the environmental temperature, while the bottom wall was subjected to a thermally insulated condition. The working fluid is air with Prandtl number (0.71). The governing equations were solved using the finite element technique and the stream line function for several magnitudes of Rayleigh numbers (103~ 106). It was found that the phenomenon of heat transfer was affected by the surface corrugation amplitude and the aspect ratio for a high value of Rayleigh number. Saha et al. [3] performed a numerical simulation for studying an inclined enclosure with corrugated side walls subjected to free convection. The bottom surface is subjected to a uniform heat flux, its length was around (20% – 80%) of the enclosure total length, while the upper surface was thermally insulated. The flowing fluid was assumed to be laminar, steady and incompressible. The finite element technique was performed to obtain an approximated solution for the Energy and the Navier Stokes equations. The effect of the inclination angle and the size of the discrete heat source on the heat transfer for various values of Grashof number (103~ 106) and Prandtl number (0.71) was studied. The authors presented their results in form of stream lines and isothermal profiles, and they concluded that an increase in the inclination angle resulted in an increase of the Nusselt Number for all the sizes of the heat sources studied. Das and Mahmud, [4] studied numerically the hydrodynamic and thermal behavior of a fluid in a square enclosure of sinusoidal corrugated upper and lower surfaces. The top surface was subjected to a constant cold temperature, whilst the bottom surface was heated under a constant temperature, and the remaining two side walls were thermally insulated. The approximated solution of the governing equations was reached by applying the finite volume method on a collected grid for the following range of Grashof number: 103~ 107; and for a Prandtl number of 1. The following amplitude/wavelength ratios were studied: 0.0, 0.05, 0.1 and 0.15, whilst the aspect ratio was maintained at 0.40. The authors concluded that the amplitude/wavelength ratio highly affects the flow behavior as well as the heat transfer rate. It was also inferred that the higher Grashof number amplifies the impact of the amplitude/wavelength ratio in the flow behavior and the heat transfer rate. Oztop et al. [5] evaluated the heat transfer by means of free convection in a wavy-walled enclosure with a discrete heat source, as well as the flow behavior inside of it by applying the finite volume method to solve the governing equations of interest. One of the two vertical walls of the enclosure was subjected to a constant cold temperature and the other was maintained at a hot and steady temperature, whereas the wavy top and bottom walls were subjected to a thermally insulated condition. The authors concluded that the effective parameters (the internal/external Rayleigh number ratios, as well as the wavy walls amplitude) has a significant impact on the heat transferred and on the fluid field. The maximum heat transfer rate occurs when the amplitude of wavy walls increases depending on the internal and external Rayleigh number. The Prandtl number becomes effective when it is bigger than one (Pr >1). Varol & Oztop [6] presented a numerical study of natural convection heat transfer and fluid flow in a shallow wavy enclosure. The bottom surface is corrugated as a sinusoidal function and heated at a constant temperature, while the top surface remains at a constant cold temperature, and the two vertical side walls were thermally insulated. The governing equations were performed and solved using CFD Research Corporation. Many parameters were studied during the simulation, such as aspect ratios, Rayleigh number and non-dimensional wave length. The simulated data were displayed as streamlines, isotherms and Nusselt number. The results revealed that the Nusselt number values increase with the aspect ratio, and the Rayleigh number increases with the decreasing of the non-dimensional wave length. Ahmed and Yovanovich [7] investigated numerically the laminar free convection heat transfer rate and the fluid flow in a two-dimensional model of a square enclosure with discrete heat sources. The source of thermal energy is placed at the center of the side walls, whilst the top and bottom surfaces are subjected to a thermally insulated condition. The governing equations were solved using the finite difference technique based on Marker and Cell (MAC) model. The working fluid was air with a Prandtl number of 0.72 and various ranges of Rayleigh numbers (0 ≤ Ra ≤ 106). The relationships between the Nusselt and Rayleigh numbers were based on the different scale lengths that were analyzed. Finally, correlations were developed for a range of Rayleigh numbers from 0 to 109. Saha et al. [8] carried out a numerical simulation in order to analyze the improvement of the heat transfer rate by means of free convection in a 2D rectangular enclosure. The bottom surface is subjected to a constant high temperature and the upper wall is considered adiabatic, while the two side walls are kept at a constant low temperature. The Finite Element Technique was applied in order to obtain an approximate solution of the Navier Stokes equations as well as of the energy equation using a triangular elements discretization scheme and a nonlinear coupled solution algorithm for all the variables involved in the flow domain. Wide ranges of Grashof number (103~ 106) were used during the simulation, and Prandtl number was taken as 0.71. The influences of aspect ratios (0.5 ~ 1) and the enclosure tilt angles (0º ~ 30º) on the fluid flow and the thermal characteristics were studied. The aim of this study is to show the ability of the numerical solution to handle the internal heat transfer problems. The results of the investigation are represented using isotherm plots and stream-lines. Also, this study presents the changes of the local Nusselt number as a function of different conditions of the heat sources. Ahmed et al. [9] analyzed, by means of numerical simulations, the impact of corrugated walls on the fluid flow behavior and in the thermal energy exchange by means of free convection in a 2D rectangular channel. Three heat sources were placed at the bottom wall of the channel, whereas the wavy shape was on the upper surface. The effect of changing the amplitude of the wave on the improvement of the heat transfer was studied. The finite volume method was used to solve the steady 2D incompressible continuity, Navier- Stokes and energy equations. Body fitted coordinates are employed to represent the complex wavy wall accurately and to adjust the velocity field. SIMPLE algorithm was used to satisfy the conservation of mass. Different ranges of wave amplitude (-0.5 ≤ A ≤ 0.5) and Reynolds number (50 ≤ Re ≤ 1000) were used for air with a Prandtl number of 0.7. The data demonstrated that the maximum improvement in the thermal energy exchange occurs when the wavy shape located on the top surface and a half wave over each source.  Al-Balushi et al. [10] numerically studied a square enclosure using the model of nonhomogeneous dynamic; they investigated the heat transfer of unsteady free convective inside. Nano fluids with magnetic nanoparticles was used as the working fluid. The enclosure cold wall was the top horizontal wall, the cavity between the two vertical walls was considered as thermally insulated, and a uniform temperature was maintained for the bottom wall of the cavity. The results showed that the volume fraction of nanoparticles and the thermal Rayleigh number intensify lead to an increase in the average Nusselt number. It was also indicated that for the nanoparticles blade shape, higher average Nusselt numbers were found. Harish [11] numerically investigated a partial square enclosure through studding the behavior of the thermal plume driven by the buoyancy. Isothermal protruding heat source block was maintained inside the enclosure. Unsteady Favre-averaged Navier-Stokes equation was used to model the turbulent flow with buoyancy modified Lam Bremhorst with low Reynolds number k−ϵ model. It was found that the location of the heat source has a strong effect on the rate of the vent bidirectional exchange, heat transfer characteristics and thermal mixing inside the enclosure. Guestal et al. [12] numerically studied the natural convection heat transfer in a cylindrical enclosure that was horizontally heated partially using two nanofluids. Using a constant temperature, the heat was applied to the lower part of the enclosure. Correlation equations were obtained for the average Nusselt number as a function of the volume fraction, the Rayleigh number, and the heated length for both types of nanofluid.

In this investigation, the focus was on the study of thermal energy exchanged released from two separate sources of heat, as well as the influence of the separation between the sources and the wave amplitude on the flow and heat transfer inside the enclosure of the corrugated ceiling.

2. Mathematical Modeling

2.1 The present study assumptions

(1) Incompressible and Newtonian fluid.

(2) 2D Laminar flow.

(3) The thermo-physical properties of the fluids are constant, except in the buoyancy term.

(4) The pressure changes of the fluid contained in the enclosure do not reach high values.

(5) The heat dissipation by viscous effects is neglected.

2.2 Governing equations

For the case of a 2D, laminar heat transfer by means of free convection inside an enclosure, as shown in Figure 1, the governing equations that describe the conservation of mass, momentum, and energy can be written as follows, respectively:

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$                                    (1)

$\rho \left( u\frac{\partial u}{\partial x}+v\frac{\partial v}{\partial y} \right)=-\frac{\partial P}{\partial x}+\mu \left( \frac{{{\partial }^{2}}u}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}u}{\partial {{y}^{2}}} \right)$            (2)

$\rho \left( u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} \right)=-\frac{\partial P}{\partial y}+\mu \left( \frac{{{\partial }^{2}}v}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}v}{\partial {{y}^{2}}} \right)+g\beta (T-{{T}_{c}})$      (3)

$\rho \left( u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y} \right)=\alpha \left( \frac{{{\partial }^{2}}T}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}T}{\partial {{y}^{2}}} \right)$                    (4)

Figure 1. Schematic diagram of the physical system

2.3 The governing equations non-dimensional form

The following are the normalized form of the governing equations that describe the conservation of mass, momentum and energy, respectively:

$\frac{\partial {{u}^{*}}}{\partial {{x}^{*}}}+\frac{\partial {{v}^{*}}}{\partial {{y}^{*}}}=0$                                   (5)

$\left( {{u}^{*}}\frac{\partial {{u}^{*}}}{\partial {{x}^{*}}}+{{v}^{*}}\frac{\partial {{v}^{*}}}{\partial {{y}^{*}}} \right)=-\frac{\partial {{P}^{*}}}{\partial {{x}^{*}}}+\Pr \left( \frac{{{\partial }^{2}}{{u}^{*}}}{\partial {{x}^{*}}^{^{2}}}+\frac{{{\partial }^{2}}{{u}^{*}}}{\partial {{y}^{*}}^{^{2}}} \right)$            (6)

$\left( {{u}^{*}}\frac{\partial {{v}^{*}}}{\partial {{x}^{*}}}+{{v}^{*}}\frac{\partial {{v}^{*}}}{\partial {{y}^{*}}} \right)=-\frac{\partial {{P}^{*}}}{\partial {{y}^{*}}}+\Pr \left( \frac{{{\partial }^{2}}{{v}^{*}}}{\partial {{x}^{*}}^{^{2}}}+\frac{{{\partial }^{2}}{{v}^{*}}}{\partial {{y}^{*}}^{^{2}}} \right)+Ra\Pr {{T}^{*}}$  (7)

   $\left( {{u}^{*}}\frac{\partial {{T}^{*}}}{\partial {{x}^{*}}}+{{v}^{*}}\frac{\partial {{T}^{*}}}{\partial {{y}^{*}}} \right)=\left( \frac{{{\partial }^{2}}{{T}^{*}}}{\partial {{x}^{*}}^{^{2}}}+\frac{{{\partial }^{2}}{{T}^{*}}}{\partial {{y}^{*}}^{^{2}}} \right)$            (8)


$\begin{align}  & {{T}^{*}}=\frac{T-{{T}_{c}}}{{{T}_{h}}-{{T}_{c}}},\quad {{P}^{*}}=\frac{P{{L}^{2}}}{\rho {{\alpha }^{2}}},\quad {{v}^{*}}=\frac{vL}{\alpha },\quad  \\ & {{u}^{*}}=\frac{uL}{\alpha },\quad {{y}^{*}}=\frac{y}{L},\quad {{x}^{*}}=\frac{x}{L} \\ \end{align}$

2.4 Boundary conditions

To solve the governing equations (mass, momentum and energy), the boundary conditions of these equations must be established, as shown in Figure 2. 

$\left. \begin{array}{*{35}{l}}   {{u}^{*}}={{v}^{*}}=0,{{T}^{*}}=0 & \text{ for }(A,B,C)  \\   {{u}^{*}}={{v}^{*}}=0,{{T}^{*}}=1 & \text{ for }(D,E)  \\   {{u}^{*}}={{v}^{*}}=0,\left( \frac{\partial {{T}^{*}}}{\partial {{y}^{*}}} \right)=0\quad \text{ } & \text{for }(F,G,H)  \\ \end{array} \right\}$      (9)

Figure 2. Boundary conditions

2.5 Average Nusselt number

The Nusselt number can be calculated by dividing the amount of heat transferred by convection over the amount of heat transferred by conduction, and these amounts are calculated near the heated walls of the enclosure. After compensation and simplification, one gets a local Nusselt number, as in the following equation: 

$Nu=\frac{hL}{k}=\frac{-k\left( \frac{\partial {{T}^{*}}}{\partial {{y}^{*}}} \right)L}{({{T}_{w}}-{{T}_{o}})k}=-\frac{\partial {{T}^{*}}}{\partial {{y}^{*}}}$                  (10)

The average Nusselt number is calculated as follows:

$N{{u}_{av}}=\frac{1}{L_{s}^{*}}\int\limits_{a}^{a+L_{s}^{*}}{Nud{{x}^{*}}}$                          (11)

3. Numerical Solution

The following is the Cartesian Coordinates form of the conservation equations [13]:

$ \frac{\partial ({{u}^{*}}{{\phi }^{*}})}{\partial {{x}^{*}}}+\frac{\partial ({{v}^{*}}{{\phi }^{*}})}{\partial {{y}^{*}}}=\frac{\partial }{\partial {{x}^{*}}}({{\Gamma }^{{{\phi }^{*}}}}\frac{\partial {{\phi }^{*}}}{\partial {{x}^{*}}})+\frac{\partial }{\partial {{y}^{*}}}({{\Gamma }^{{{\phi }^{*}}}}\frac{\partial {{\phi }^{*}}}{\partial {{y}^{*}}})+S_{\phi }^{*}$   (12)


Γ: Effective diffusion coefficient.

f: Geherao dependent variable.

Sf: Source term.

The discretization scheme used was based on the elliptic partial differential equations in order to produce curvilinear coordinates [14], i.e. a method of approximating the differential equations by a system of algebraic equations for the variables at some set of discrete locations in space and time. There are many approaches, but the most important of which are: Finite Element Method (FEM), Finite Volume Method (FVM), and Finite Difference Method (FDM). When the grid is very fine, the same solution was found from each type of these methods. However, for some classes of problems, some methods are more suitable than others. The use of the finite element method was through to be much better since it offers greater flexibility, especially for difficult geometry. Finite difference methods are widely used and are the more common, finite difference may be derived either by using the Taylor expansion polynomial approximation or by the use of the finite volume scheme. Eq. (12) is converted into a computational domain by using the transformation: ${{\eta }^{*}}={{\eta }^{*}}\left( {{x}^{*}},{{y}^{*}} \right),{{\xi }^{*}}={{\xi }^{*}}\left( {{x}^{*}},{{y}^{*}} \right)$ and the chain rule to convert the derivatives $\left( \frac{\partial }{\partial {{x}^{*}}},\frac{\partial }{\partial {{y}^{*}}} \right)$ [15], also the final form of the transformed equation can be written as:

$\frac{\partial }{\partial {{\xi }^{*}}}\left( G_{1}^{*}{{\phi }^{*}} \right)+\frac{\partial }{\partial {{\eta }^{*}}}\left( G_{2}^{*}{{\phi }^{*}} \right)=\frac{\partial }{\partial {{\xi }^{*}}}\left( {{\Gamma }^{{{\phi }^{*}}}}{{J}^{*}}\alpha \frac{\partial {{\phi }^{*}}}{\partial {{\xi }^{*}}} \right)+\frac{\partial }{\partial {{\eta }^{*}}}\left( {{\Gamma }^{{{\phi }^{*}}}}{{J}^{*}}{{\gamma }^{*}}\frac{\partial {{\phi }^{*}}}{\partial {{\eta }^{*}}} \right)+S_{\xi ,\eta }^{*}+\frac{S_{\phi }^{*}}{{{J}^{*}}}$ (13)


G1and G2: Contra-variant velocity components.

J : Transformation Jacobian.

α, β, γ: Transformation coefficients.

The transformation coefficients can be written as follows:

$\left. \begin{align}  & {{\alpha }^{*}}={{\left( \frac{\partial {{x}^{*}}}{\partial {{\eta }^{*}}} \right)}^{2}}+{{\left( \frac{\partial {{y}^{*}}}{\partial {{\eta }^{*}}} \right)}^{2}} \\ & \beta {}^{*}=\left( \frac{\partial {{x}^{*}}}{\partial {{\eta }^{*}}}\frac{\partial {{x}^{*}}}{\partial {{\xi }^{*}}} \right)+\left( \frac{\partial {{y}^{*}}}{\partial {{\eta }^{*}}}\frac{\partial {{y}^{*}}}{\partial {{\xi }^{*}}} \right) \\ & {{\gamma }^{*}}={{\left( \frac{\partial {{x}^{*}}}{\partial {{\xi }^{*}}} \right)}^{2}}+{{\left( \frac{\partial {{y}^{*}}}{\partial {{\xi }^{*}}} \right)}^{2}} \\ \end{align} \right\}$                   (14)

When Eq. (9) is integrated over a control volume of the computational domain, its convective term is discretized when a hybrid scheme is applied. On the other hand, the diffusion term is discretized with the usage of a central scheme. Also, the SIMPLE algorithm is applied to find an adjustment of the velocity field, with the principal aim of satisfying the law of mass conservation. Due to the fact that all the variables involved in the study are stored in the center of each cell of the discretized flow domain, the method of interpolation is applied in order to correct the pressure and avoid a decoupling between the pressure and the fluid velocity [16]. The solution of the discretized equations that result from this process is obtained by following an iterative procedure known as the Tri-Diagonal Matrix Algorithm. The TDMA (Tri-Diagonal Matrix Algorithm) enables the solution of a set of simultaneous, linear, algebraic equations. The tri-diagonal matrix algorithm (TDMA) is a simplified form of Gaussian elimination that can be used to solve the tri-diagonal system of equations. It consists of two parts, a forward elimination phase and a backward substitution phase. Examples of such matrices commonly arise from the natural cubic spline interpolation discretization and one-dimensional Poisson equation. In general, Thomas' algorithm is not stable, but in several special cases, it is so, such as when the matrix is symmetric positive definite or diagonally dominant (either by columns or by rows). In the general case, if stability is required, Gaussian elimination with partial pivoting (GEPP) is recommended instead. The convergence is reached when the residuals of each one of the governing equations reach a value inferior than 0.001.

The results of the current study were obtained by converting the previously mentioned equations into a computer program in the language (FORTRAN 90). FORTRAN came to take over this programming area in computationally intensive areas, such as computational chemistry, crystallography, computational physics, computational fluid dynamics, analysis of finite element, and prediction of the numerical weather. For computing of high-performance, it is a common language and it is used for programs that rank and benchmark the world's fastest supercomputers.

3.1 Finite-difference approximation error

Because the numerical solutions consist of an approximate solution of unsolvable partial differential equations, errors are introduced in the approximation of the derivative terms. First and second derivative terms are approximated through the use of a Taylor series expansion. Therefor,

$\frac{\partial \varphi }{\partial x}=\frac{{{\varphi }_{i+1}}-{{\varphi }_{i}}}{\Delta x}+E\left( \Delta x \right)$                          (15)

$\frac{{{\partial }^{2}}\varphi }{\partial {{x}^{2}}}=\frac{{{\varphi }_{i+1}}-2{{\varphi }_{i}}+{{\varphi }_{i-1}}}{{{\left( \Delta x \right)}^{2}}}+E\left( \Delta {{x}^{2}} \right)$       (16)

While Eq. (15) displays the first derivative via the forward difference method, a similar expression may be used for the backward difference the terms $E\left( \Delta x \right)$ and $E\left( \Delta {{x}^{2}} \right)$ represented the errors associated with each approximation. Thus, first and second-order accuracy is obtained for the first and second derivatives, respectively.

3.2 Grid independent solution test

Different tests were carried out for testing the independency solution of mesh. For this purpose, three grid systems of (80 × 80) nodes, (100 × 100) nodes and (120 × 120) nodes were adopted for specific cases. The variation trend in figure (3) implies that the discretized system with (100 × 100) nodes can represent a good quality grid for this problem, since there is no more noticeable improvement in the results with respect to grid number and calculation time enhancement, where the variation in local Nusselt number values is less than (1%) for the grid (120 × 120) nodes. 

Figure 3. The grid independent solution test for source (1) D*= 0.4, Ra=105

4. Results and Discussions

4.1 Effect of distance between discrete heat sources on (Nuav)

Figures (4) to (7) elucidate the streamline function by the left side and isotherm plots to the right side for working fluid at different distances between the heat sources and several wave amplitudes with wave number (N=3). Figures (8) to (9) eince that an increase in the distance between two heat sources results in a raise of the average Nusselt number to double for all studied cases and heat sources because of the increasing in temperature gradient at the inner and outer edges of the heat sources, as it seems from the isotherm plots in Figures (4) to (7), in which it can be appreciated that the thermal energy exchange rate raises its value. This leads to the conclusion that the distance between the heat sources is an important factor for controlling the field of flow within the enclosure and the process of heat transfer from thermal sources.

The results also showed that the average Nusselt number increases with increasing the Rayleigh number for all cases and for different distances as shown in Figures (8) to (9), because of the thermal energy exchange process between the sources of heat and the layers of fluid (air) touching them are by conduction at low Rayleigh number (103), as shown in Figure (4). Where, the isotherm lines are semi-circular and parallel, but it's when the Rayleigh number is increasing to (104), the heat transfer process is carried out by convection and the heat transfer rate becomes greater at Rayleigh number (105), as shown in the isotherm plots of Figure (5).

Figures (4) to (7) illustrate the thermal behavior for fluid (air) inside the enclosure. It can be noticed from the streamline plots shown to the right side formation of swirls over the discrete heat sources due to the fluid contact with the thermal sources, the fluid is heated and has less density, and because of the influence of buoyancy and gravitational force, there are two large vortices over each thermal source. These two vortices touch each other at the midpoint of the enclosure, the fluid flows above the first thermal source to the left of the space counterclockwise, and runs over the second thermal source to the right of the space clockwise. This behavior is the same for all cases and for all Rayleigh numbers (103-105). Note that the increase in the flow stream function values with the Rayleigh number increase is due to increase of the buoyancy force, as shown in the figures from (4) to (7). All swirls formed symmetrically around the center axis of the enclosure as well as the isotherm lines, which are distributed in parallel to each other and have a semi-circular shape around each thermal source as shown in the shapes to the left, and all the values of distance between the thermal sources and the cases shown in Figures (4) to (7) are due to symmetry space.

The results demonstrated that the flow force of the fluid increases with the decrease in the value of the distance between the two thermal sources, since the value of the flow stream function is greater for a distance of (D* = 0.4). At a distance of (D* = 0.6) and for all cases, the values of the average Nusselt number are large compared to the average Nusselt number at the distance (D* = 0.4) for all Rayleigh numbers. This increase is resulting from the approach of the two cold vertical walls leading to shortening the way of heat transfer, where the boundary layer thickness becomes thinner when approaching the heat sources from the vertical walls, so the average Nusselt number increases, as shown in Figures (8) and (9). It is also noted that it increases with the increase of the Rayleigh number due to the large amount of temperature gradient at the outer edge near the cold vertical wall where the differences of temperature between the vertical wall and the thermal source increases, leading to an increase in the rate of heat transfer for both two heat sources as well as a raise in the temperature gradient at the inner edge of both two heat sources. And, this is because of the distance between the sources, each thermal source acts as a source independent of the other, and the temperature gradient increases between the inner edge and the fluid, so the rate of heat transfer also increases.

The distance (D*=0.6) is a special case at Rayleigh number (Ra=105), in spite of the increase in the average number of Nusselt with increasing the number of Rayleigh for all cases, but the increase at the number of Rayleigh (Ra=105) is less than the rate of increase registered at a distance of (D*=0.4) because of the temperature gradient at the inner and outer edges of the heat sources which is less at the distance (D*=0.6) and (Ra=105). At this distance, the two heat sources are far from each other, and are far from the two cold vertical walls on the other side. Therefore, the temperature gradients at the inner edges are slightly compared to (D*=0.4), so the heat transfer rate decreases at (Ra=105), as shown in the Figures (8) and (9).

4.1 Effect of wave amplitude on (Nuav)

A numerical study was conducted showing that the number of Nusselt was affected by the wave amplitude according to the value of the distance between the thermal sources. Figure (10) at (D*=0.4) depicts the change in the average Nusselt number at Rayleigh number (Ra=105) and the change in the amplitude between (0.1-0.2) for the pair of discrete heat sources. The average Nusselt number increases as the wave amplitude increases, because the temperature stratification of isotherm lines in the middle of the space and the thermal separation area began to decrease with an increase of the wave amplitude and the heat transfer rate increased at Rayleigh number (Ra=105), this appears from the stream function lines and isotherm plots in Figures (4) to (7).

Finally, figure (11), which illustrates the results at (D*=0.6), exhibits the raise in the Nusselt number with an amplitude of the wave that increases for a Rayleigh numbers of (103-104) for the two heat sources due to temperature stratification in the isotherm lines over the heat sources and the increasing of temperature gradient on the outer edges of the thermal sources, as shown in Figures (4) and (6). However, this increase does not continue, it results in an oscillation of the magnitudes of the Nusselt number at a Rayleigh Number of 105.

Figure 4. Streamline and isotherm plots at (Ra = 103) and (A* = 0.1) for various distances between heat sources

Figure 5. Streamline and isotherm plots at (Ra = 105) and (A* = 0.1) for various distances between heat sources

Figure 6. Streamline and isotherm plots at (Ra = 103) and (A* = 0.2) for various distances between heat sources

Figure 7. Streamline and isotherm plots at (Ra = 105) and (A* = 0.2) for various distances between heat sources

Figure 8. Relation between (Nuav) and Rayleigh number for (source 1) and (D*=0.4, 0.6) at A*=0.1

Figure 9. Relation between (Nuav) and Rayleigh number for (source 1) and (D*=0.4, 0.6) at A*=0.2

5. Conclusions

The natural convection process was analyzed by means of numerical simulations applying the finite volume method in a bi-dimensional enclosure with a rectangular shape, in which the upper wall has a wavy form, the side walls are kept at a low temperature, and the lower wall is maintained at an elevated temperature. The principal aim of this investigation is to reach quantitative results that provide information about the cooling effects. It is necessary to indicate that the principal parameters that have a great interest are the Rayleigh number, the dimensionless separation between the sources of heat, and the amplitude of the wave. The results of the study helped in reaching the following conclusions:

An increase in the separation between the thermal sources results in a raise of the average Nusselt number for all values of Rayleigh number.

The average Nusselt number increases when the magnitude of the wave amplitude gets higher for all distances between the heat sources.



Cartesian Coordinates

Component of the velocity in the x axis (m/s)


Component of the velocity in the y axis (m/s)

u*, v* 

Dimensionless velocity component in the x and y direction


Temperature degree (K)


Rayleigh number

Number of waves


Average Nusselt number


Length of the heat source (m)

Length of the enclosure (m)

Fluid thermal conductivity (W/m.K)


Heat transfer coefficient (W/m2.K)

Enclosure height (m)


Grashof number


Acceleration due to gravity (ms−2 )


Wave amplitude (m)

Greek symbols


Density (kg/m3)


Kinematic viscosity (N.s/m2)


Volume expansion (K-1)


Thermal diffusivity (m2/s )



Dimensionless quantity










[1] Ataei-Dadavi, I., Rounaghi, N., Chakkingal, M., Kenjeres, S., Kleijn, C.R., Tummers, M.J. (2019). An experimental study of flow and heat transfer in a differentially side heated cavity filled with coarse porous media. International Journal of Heat and Mass Transfer, 143: 118591.

[2] Hasan, M.N., Saha, S., Saha, S.C. (2012). Effects of corrugation frequency and aspect ratio on natural convection within an enclosure having sinusoidal corrugation over a heated top surface. International Communications in Heat and Mass Transfer, 39(3): 368-377.

[3] Saha, S., Sultana, T., Saha, G., Rahman, M.M. (2008). Effects of discrete isoflux heat source size and angle of inclination on natural convection heat transfer flow inside a sinusoidal corrugated enclosure. International Communications in Heat and Mass Transfer, 35(10): 1288-1296.

[4] Das, P.K., Mahmud, S. (2003). Numerical investigation of natural convection inside a wavy enclosure. International Journal of Thermal Sciences, 42(4): 397-406.

[5] Oztop, H.F., Abu-Nada, E., Varol, Y., Chamkha, A. (2011). Natural convection in wavy enclosures with volumetric heat sources. International Journal of Thermal Sciences, 50(4): 502-514.

[6] Varol, Y., Oztop, H.F. (2006). Free convection in a shallow wavy enclosure. International Communications in Heat and Mass Transfer, 33(6): 764-771.

[7] Ahmed, G.R., Yovanovich, M. (1992). Numerical study of natural convection from discrete heat sources in a vertical square enclosure. Journal of Thermophysics and Heat Transfer, 6(1): 121-127.

[8] Saha, G., Saha, S., Islam, M.Q., Akhanda, M.R. (2007). Natural convection in enclosure with discrete isothermal heating from below. Journal of Naval Architecture and Marine Engineering, 4(1): 1-13.

[9] Musatfa, A.W., Khalif, H.J., Ali, H.H. (2011). Effect of sinusoidal wavy wall on heat transfer from discrete heat sources placed in two-dimensional channel. Al-Qadisiya Journal for Engineering Sciences, 4(4): 408-418.

[10] Al-Balushi, L.M., Uddin, M.J., Rahman, M.M. (2019). Natural convective heat transfer in a square enclosure utilizing magnetic nanoparticles. Propulsion and Power Research, 8(3): 194-209.

[11] Harish, R. (2018). Buoyancy driven turbulent plume induced by protruding heat source in vented enclosure. International Journal of Mechanical Sciences, 148: 209-222.

[12] Guestal, M., Kadja, M., Hoang, M.T. (2018). Study of heat transfer by natural convection of nanofluids in a partially heated cylindrical enclosure. Case Studies in Thermal Engineering, 11: 135-144.

[13] Ferziger, J.H., Perić, M. (2002). Computational Methods for Fluid Dynamics. 3: 196-200. Berlin: Springer.

[14] Hoffmann, K.A. (1989). Computational Fluid Dynamics for Engineers. Engineering Education System, Austin, Tex.

[15] Thompson, J.F., Warsi, Z.U., Mastin, C.W. (1985). Numerical grid generation: foundations and applications. 45. Amsterdam: North-holland.

[16] Rhie, C.M., Chow, W.L. (1983). Numerical study of the turbulent flow past an airfoil with trailing edge separation. AIAA journal, 21(11): 1525-1532.