Numerical Study of Natural Convection Flow in Rectangular Cavity with Viscous Dissipation and Internal Heat Generation for Different Aspect Ratios

Abstract

Numerical simulations have been performed to investigate the influence of constant volumetric heat generation and viscous dissipation on the unsteady natural convection flow of an incompressible Newtonian fluid contained in a rectangular cavity. The left vertical wall of the cavity is cooled, while the right vertical wall is heated, and the bottom and top walls are adiabatic. A numerical technique based on the implicit finite difference method (IFDM), along with an upwind finite difference scheme and an iterative successive over relaxation (SOR) technique, is employed to solve the governing equations numerically. The effect of physical parameters, namely the modified Rayleigh number (${10}^3\le Ra\le {10}^7$), aspect ratio ($1\le A\le 4$), Prandtl number ($Pr=0.7$, $1.0$, $6.2$, 15), volumetric internal heat generation parameter ($Q_{\lambda }=0$, 1), and Eckert number ($0\le Ec\le {10}^{-6}$), on the streamlines and isotherms are discussed graphically. Variations of maximum stream function, as well as average and local Nusselt number, are also discussed. The results show that the increase in Eckert number from 0 to ${10}^{-4}$ causes the average heat transfer to decrease, while $Pr=0.71$, $Ra={10}^4$, and $Q_{\lambda }=0$. Additionally, the average heat transfer decreases as the cavity width increases from 1 to 4, while $Pr=0.71$, $Ra$ = $5\times {10}^4$, $Ec={10}^{-6}$ and $Q_{\lambda }=1$. The results of the numerical model used here are in excellent accord with earlier findings.

1. Introduction

The study of heat transfer is considered important as it attempts to optimize the efficiency of energy devices, and researchers try to find new solutions to enhance heat transfer. Natural convection in enclosures under different influencing parameters, such as internal heat generation and viscous dissipation, is one of the various techniques that are applied to strengthen the thermal performance of power systems. This phenomenon is the subject of intensive research in view of its variety of applications in several engineering processes and technological areas, like heat exchangers, nuclear reactors, thermal insulation engineering, cooling of microelectronic components, ventilation, heating control in building design, etc. This importance has led to an enormous amount of research, both theoretical and experimental. Recently, natural convection in rectangular cavities has become a subject of interest for different researchers [1,2,3,4,5,6] owing to the breadth of its applications. Many geometries, less or more complicated, have been investigated by experimental, numerical, or theoretical approximations. The difficulties often experienced in developing experimental research to determine the parameters associated with convection have led to various approaches for solving natural convection problems. The increase in computational power and the development of numerical algorithms implemented with systems of differential equations have made the numerical solution to these problems technically and computationally feasible and enables prediction of the results.

The natural convection in a cavity is caused either by the existence of internal heat sources or by differences in temperature across the cavity boundary walls. In certain technological applications, natural convection in enclosures with internal heat generation (IHG) is of major importance. The natural convection phenomenon caused by internal heating occurs in the problems of geothermal reservoir engineering, in the movement of fluids in the earth, melting of cores in the reactor sand, and many others. A review of literature shows that an immense amount of research has been conducted on the subject of heat transfer caused by IHG both in channels and in enclosed domains. Earlier studies related to IHG in differentially heated square cavities include the work of Acharya and Goldstein [7]. The governing equations were solved with a finite difference approach called the SIMPLER method. For a fluid with $R{a}_I$ (internal Rayleigh number) and $R{a}_E$ (external Rayleigh number) between ${10}^3$ and ${10}^7$, $Pr$ (Prandtl number) = 0.7, and inclination angle of cavity varied from ${30}^{\circ }$ to ${90}^{\circ }$. They found an important result that the flow was directed downward along both the cooled and heated vertical walls when $R{a}_I$ was remarkably greater than $R{a}_E$. Rahman and Sharif [8] examined the two-dimensional (2D) laminar natural convection in rectangular cavity of cooled bottom, heated top and insulated vertical walls both without and with heat generation. The finite-difference method (FDM) was used for numerical investigation with aspect ratios (A) ranging from 0.25 to 4, and $R{a}_I$ and $R{a}_E$ were $2\times {10}^5$ at different angles of inclination. It was observed that the impact of angle of inclination on the isothermal walls is especially strong in the upstream regions. Their results showed that when $R{a}_E/R{a}_I>1$, heat and flow transfer by convection is the same as that in an enclosure in the absence of internal heat generating fluid, and similar results were noticed as in Acharya and Goldstein [7]. They also found that the strength of convection increased as the aspect ratio increased. Churbanov et al. [9] showed that in a rectangular cavity, changing wall boundary conditions from isothermally cold to adiabatic reduces the oscillation’s amplitude at the beginning of instability. Additionally, the stability range increases by increasing the cavity width with the influence of IHG.

Hossain and Rees [10] investigated the laminar natural convection flow of water exposed to density inversion in a rectangular enclosure with IHG. The flow in the enclosure was generated by differential heating of the side walls and supported by IHG. The values of $Ra$ and $Pr$ were fixed at ${10}^5$ and 11.58, respectively. It was observed that the circulation of the flow was reversed when the parameter of IHG ($\lambda$) was sufficiently strong. Grosan et al. [11] examined the magnetohydrodynamics of free convection in a rectangular porous enclosure in the existence of IHG for various Rayleigh numbers and aspect ratio of the enclosure. They found that the examined parameters have a significant impact on the intensity of core convection. Additionally, the study reveals that the local Nusselt number decreases on the bottom wall when the inclined angle increases (where the magnetic field shifts from a horizontal to a vertical direction), while the opposite occurs on the top wall of the cavity. Joshi et al. [12] analyzed natural convection in a rectangular enclosure by using the semi-analytical method. The results were reported on temperature and velocity profiles for two various boundary conditions, namely adiabatic bottom and top walls, and isothermal side walls, as well as all isothermal walls. Molla et al. [13] and Molla et al. [14] examined natural convection flow along a uniformly heated perpendicular wavy surface and a horizontal cylinder in used for heat generation. They observed that in the heat generation scenario, the thermal state of the fluid intensifies, causing a decrease in the Nusselt number. In this situation, the buoyancy force increases, which results in an augmented flow rate within the boundary layer. As a result, both the velocity boundary layer and the thermal boundary layer become thicker. Conversely, the opposite phenomenon takes place in the case of heat absorption. Cheong et al. [15] conducted a detailed study on natural convection in a sinusoidally heated, curved, porous enclosure in the existence of IHG/absorption. They observed that the system’s heat transfer increases due to the wavy nature of the cavity; however, by increasing the parameter of IHG/absorption, rate of heat transfer in the enclosure decreases. Berrahil et al. [16] numerically examined the natural convective flow in a laterally heated enclosure that contains electrically conductive fluid in the existence of an external magnetic field and IHG. The results revealed that the average $Nu$ decreases with increasing heat generation and increases with increasing Grashof number ($Gr$). They found that at low $S_Q$ (volumetric heating), the average Nusselt number ($N{u}_{avR}$) near the hot plate decreases, while it increases near the cold plate with the increase in heat generation. Hdhiri et al. [17] studied natural convection flow in a differentially heated rectangular enclosure containing porous medium in the existence of IHG with an electrically conductive fluid $(0.015\le Pr\le 0.054)$. A comparison between the porous medium and homogeneous medium was presented. They demonstrated that the average Nusselt number ($N{u}_{avr}$) is greatly influenced by the increase in $Pr$, and the existence of homogeneous media overestimates the rates of heat transfer more than the existence of a porous medium. Furthermore, they established a correlation for heat transfer rates within the porous medium.

Viscous dissipation is the heat energy generated due to friction between adjacent layers of fluid. It is an irreversible process that involves the conversion of work performed by velocity against viscous forces into internal energy (heat), leading to an increase in the temperature of the fluid. Viscous dissipation (VD) heat becomes significant in natural convective flow when the flow field is of intense size or in a high gravitational field, or at a low temperature. Such effects are also important in some industrial operations and geophysical flows and are typically represented by a parameter called the Eckert number and is denoted by ($Ec$). All of these significant subjects have obtained considerable attention from several researchers. Mahajan [18] conducted a study on natural convection flows with viscous heating dissipation. They reported that heat transfer rates decrease by increasing the parameter of dissipation. Hossain [19] investigated the influence of viscous and joule heating on free convection flow of an incompressible viscous and electrically conducting fluid over a semi-infinite plate in the presence of a constant magnetic field. With the distance from the leading edge, the temperature of the semi-infinite plate varied linearly. Saeid and Pop [20] numerically studied the impact of VD in a porous cavity on natural convection. The results showed that by increasing the parameter of VD, the rate of heat transfer at the hot wall decreases. Al-Badawi and Duwairi [21] investigated numerically the magnetohydrodynamics of natural convection heat transfer in an iso-flux inclined rectangular porous enclosure under the impact of joule and viscous heating. An iso-heat flux was employed to cool and heat the two opposite sides of the cavity, while the other two sides remained adiabatic. Their study concluded that the effects of joule and viscous heating decreased the heat transfer rates. Hossain et al. [22] investigated the influence of VD in a thermally stratified medium on free convection from a perpendicular plate with a constant surface heat flux. They reported that the temperature of the fluid in the boundary layer region decreases by increasing the parameter of VD.

Israel-Cookey et al. [23] analysed the radiation and VD effects on magneto-hydrod-ynamic unsteady natural convection flow past a vertical plate within a porous medium. They showed that when viscous dissipation increases, the temperature profile increases. Considering the presence of nanofluids within the porous medium, RamReddy et al. [24] and Chamkha et al. [25] studied the influence of VD on convective heat transfer of nanofluids on the boundary layer. Badruddin et al. [26] analyzed the impact of VD and radiation on the heat transfer rate from a square enclosure containing porous media. They noticed that the $N{u}_{avR}$ at the hot wall along the height of the cavity decreases when the VD parameter increases. However, the local Nusselt number at the cold wall increases with the increase in the parameter of VD. Mughees et al. [27] examined the influence of the Soret and Dufour impacts on the bidirectional stretched flow of a Casson fluid with magneto-hydrodynamic effects. The analysis takes into account the effects of viscous dissipation, thermal radiation, and joule dissipation. The results indicate that the slip parameters have a significant impact on the flow field. Specifically, the Dufour effect leads to an increase in the fluid temperature, whereas the Soret effect results in a reduction in the temperature. Fahim et al. [28] conducted a numerical study of time-dependent blood flow in a $\omega$-shaped stenosed vessel, taking into account the influence of body acceleration and a pulsatile pressure gradient. The results revealed that as the amplitude of body acceleration increased, both the flow rate and blood velocity showed an upward trend. However, when the stenosis height, Weissenberg number, and power-law exponent were increased, the opposite effect was observed. Furthermore, the study found that the Brinkman number and Prandtl number had significant impacts on the temperature field. Sadiq et al. [29] conducted a flow analysis of nanofluid near a stagnation point over a lubricated stretching surface. The lubrication was achieved using a thin layer of power-law fluid with a variable thickness. They investigated the effects of thermophoresis and Brownian motion. The findings revealed that the addition of lubrication led to increased fluid velocity and decreased the temperature of the nanofluid at the stretchable surface.

The flow of fluid and heat transfer distributions in rectangular enclosure are highly influenced by the $Ra$ and aspect ratio (A) of the enclosure. For a 2D cavity, A is the ratio of cavity width to cavity height. Falahat [30] studied the impact of A and $Ra$ on natural convection in a partially heated enclosure for 0.5 $\le A\le$ 4 and ${10}^3\le Ra\le {10}^6$. They reported that the transfer of heat increases over the entire range of A with the increase in $Ra$. The results showed that $Nu$ increases for A from 0.5 to 1, and beyond that it decreases. Nogueira [2] investigated natural convection in rectangular enclosures with various A. They found that as the dimensionless ratio $L/D$ (Length/Height) increases, the $Nu$ increases as well.

Extensive research has been conducted on natural convection in various configurations, such as enclosures or cavities, with heat generation and viscous dissipation considered separately. These studies have provided valuable insights into heat transfer mechanisms and flow behavior. Generally, the physical phenomenon are not simple and many tangible factors on which physical processes depend. So, to better visualize this, we consider as many factors as possible, one of which is a viscous dissipation which is very important when it comes to measuring internal energy of the system. Natural convection with viscous dissipation is neglected in the literature. However, further investigation is required to explore the combined influence of all these factors, specifically focusing on the impact of cavity width. In this case, we consider that two vertical walls are kept at various isothermal temperatures and the horizontal walls are adiabatic. The reduced non-dimensional equations governing the flow are numerically studied by using the implicit finite difference method (IFDM), together with the iterative successive over-relaxation (SOR) technique [31,32]. The final results are represented in the form of streamlines, isotherms, as well as local and average heat transfer, which express the combined influence of heat due to constant volumetric heat generation and viscous dissipation for cavities with differential heating of the vertical walls.

2. Problem Formulation

Consider the 2D unsteady natural convection flow of an incompressible Newtonian fluid contained in a rectangular cavity. The cavity is developed by the regions between the two vertical walls at $\widehat{x}=0$ and $\widehat{x}=L$, and the two horizontal walls at $\widehat{y}=0$ and $\widehat{y}=H$. The temperature of the vertical left side is isothermally maintained at a cold temperature ${\widehat{T}}_C$ at $t=0$. The temperature of the right vertical wall is suddenly raised to ${\widehat{T}}_H$ and then is supposed to be isothermally maintained thereafter, where ${\widehat{T}}_H>{\widehat{T}}_C$. The lower and upper boundary walls of the rectangular cavity are held adiabatic. All solid boundary walls are rigid, no-slip walls. We also take into account the impact of heat in the flow field due to constant volumetric heat generation and viscous dissipation. Gravitational force direction is shown in the negative y-direction. Finally, with the assumption that fluid have constant characteristics with the exception of density, which varies with temperature according to the Boussinesq approximation. The flow pattern and the coordinate system of the problem are displayed in Figure 1.

Thus, equations governing the flow under the above assumption for conservation of mass, momentum, and energy in a 2D Cartesian coordinate system can be written as shown below [33]:

$$\begin{array}{c}\hfill \frac{\partial \widehat{u}}{\partial \widehat{x}}+\frac{\partial \widehat{v}}{\partial \widehat{y}}=0\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\partial \widehat{u}}{\partial \widehat{t}}+\widehat{u}\frac{\partial \widehat{u}}{\partial \widehat{x}}+\widehat{v}\frac{\partial \widehat{u}}{\partial \widehat{y}}=-\frac{1}{\rho }\frac{\partial \widehat{p}}{\partial \widehat{x}}+\nu \left(\frac{{\partial }^2\widehat{u}}{\partial {\widehat{x}}^2}+\frac{{\partial }^2\widehat{u}}{\partial {\widehat{y}}^2}\right)\end{array}$$

(2)

$$\begin{array}{c}\hfill \frac{\partial \widehat{v}}{\partial \widehat{t}}+\widehat{u}\frac{\partial \widehat{v}}{\partial \widehat{x}}+\widehat{v}\frac{\partial \widehat{v}}{\partial \widehat{y}}=-\frac{1}{\rho }\frac{\partial \widehat{p}}{\partial \widehat{y}}+\nu \left(\frac{{\partial }^2\widehat{v}}{\partial {\widehat{x}}^2}+\frac{{\partial }^2\widehat{v}}{\partial {\widehat{y}}^2}\right)+g\beta (\widehat{T}-{\widehat{T}}_C)\end{array}$$

(3)

$$\begin{array}{c}\hfill \frac{\partial \widehat{T}}{\partial \widehat{t}}+\widehat{u}\frac{\partial \widehat{T}}{\partial \widehat{x}}+\widehat{v}\frac{\partial \widehat{T}}{\partial \widehat{y}}=\alpha \left(\frac{{\partial }^2\widehat{T}}{\partial {\widehat{x}}^2}+\frac{{\partial }^2\widehat{T}}{\partial {\widehat{y}}^2}\right)+\frac{\mu }{\rho C_p}\widehat{\varphi }+\frac{S_{\lambda }}{\rho C_p}\end{array}$$

(4)

In the above equations, $\widehat{u}$ and $\widehat{v}$ represent the fluid velocity components along the $\widehat{x}$ and $\widehat{y}$-axis, respectively. The variable g denotes the magnitude of acceleration due to gravity, $\widehat{T}$ refers to the fluid temperature, $\rho$ represents the fluid’s density, $\widehat{p}$ denotes the pressure, $\widehat{t}$ is the time, $\nu =\frac{\mu }{\rho }$ denotes the kinematic viscosity of fluid, $\beta$ is the coefficient of thermal expansion, $\alpha =\frac{k}{\rho C_p}$ is the thermal diffusivity, in which k is the coefficient of thermal conductivity and $C_p$ is the molar specific heat at a constant pressure, $\mu$ is the dynamic viscosity, $S_{\lambda }$ denotes volumetric heat generation rate and viscous dissipation function $\widehat{\varphi }$ can be expressed as:

$$\begin{array}{c}\hfill \widehat{\varphi }=2\left[{\left(\frac{\partial \widehat{u}}{\partial \widehat{x}}\right)}^2+{\left(\frac{\partial \widehat{v}}{\partial \widehat{y}}\right)}^2\right]+{\left(\frac{\partial \widehat{u}}{\partial \widehat{y}}+\frac{\partial \widehat{v}}{\partial \widehat{x}}\right)}^2\end{array}$$

(5)

In accordance with the problem description, initial and boundary conditions are stated as [34]:

At $\widehat{t}<0$

$$\begin{array}{c}\hfill \widehat{u}=\widehat{v}=\widehat{T}=00\le \widehat{y}\le H0\le \widehat{x}\le L\end{array}$$

At $\widehat{t}\ge 0$

$$\begin{array}{c}\hfill \widehat{u}=\widehat{v}=0,\widehat{T}={\widehat{T}}_C\widehat{x}=00\le \widehat{y}\le H\end{array}$$

$$\begin{array}{c}\hfill \widehat{u}=\widehat{v}=0,\widehat{T}={\widehat{T}}_H\widehat{x}=L0\le \widehat{y}\le H\end{array}$$

$$\begin{array}{c}\hfill \widehat{u}=\widehat{v}=0,\frac{\partial \widehat{T}}{\partial \widehat{y}}=0\widehat{y}=00\le \widehat{x}\le L\end{array}$$

$$\begin{array}{c}\hfill \widehat{u}=\widehat{v}=0,\frac{\partial \widehat{T}}{\partial \widehat{y}}=0\widehat{y}=H0\le \widehat{x}\le L\end{array}$$

(6)

Now we introduce the stream function-vorticity formulation to remove the pressure terms from the above system of equations by using the relation given below

$$\begin{array}{c}\hfill \widehat{\omega }=\frac{\partial \widehat{v}}{\partial \widehat{x}}-\frac{\partial \widehat{u}}{\partial \widehat{y}}\end{array}$$

(7)

whereas in the form of stream function $\widehat{\psi }$, the components of velocity $\widehat{u}$ and $\widehat{v}$ are described by the expressions

$$\begin{array}{c}\hfill \widehat{u}=\frac{\partial \widehat{\psi }}{\partial \widehat{y}},\widehat{v}=-\frac{\partial \widehat{\psi }}{\partial \widehat{x}}\end{array}$$

(8)

Considering H as the cavity reference height and $U_o=\frac{\nu }{H}$ as reference velocity, we now apply the transformations given below to transform all the variables into non-dimensional form

$$\begin{array}{c}\hfill x=\frac{\widehat{x}}{H}y=\frac{\widehat{y}}{H}u=\frac{\widehat{u}H}{\nu }v=\frac{\widehat{v}H}{\nu }\psi =\frac{\widehat{\psi }}{\nu }\end{array}$$

$$\begin{array}{c}\hfill \omega =\frac{\widehat{\omega }{H}^2}{\nu }\theta =\frac{k}{S_{\lambda }{H}^2}(\widehat{T}-{\widehat{T}}_C)t=\frac{\widehat{t}\nu }{H^2}\end{array}$$

(9)

where x and y are non-dimensional coordinate axes, t, $\psi$, $\theta$, and $\omega$ are, respectively, non-dimensional time, stream function, temperature, and vorticity in the flow field, whereas u and v denote the components of velocity and $S_{\lambda }$ refers to the parameter of a constant volumetric heat generation rate. Equations (1)–(4) in dimensionless form can be written as

$$\begin{array}{c}\hfill \frac{{\partial }^2\psi }{\partial x^2}+\frac{{\partial }^2\psi }{\partial y^2}=-\omega \end{array}$$

(10)

$$\begin{array}{c}\hfill \frac{\partial \omega }{\partial t}+u\frac{\partial \omega }{\partial x}+v\frac{\partial \omega }{\partial y}=\left(\frac{{\partial }^2\omega }{\partial x^2}+\frac{{\partial }^2\omega }{\partial y^2}\right)+\frac{Ra}{Pr}\frac{\partial \theta }{\partial x}\end{array}$$

(11)

$$\begin{array}{c}\hfill \frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\frac{1}{Pr}\left(\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}\right)+Ec\left\{2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2\right\}+\frac{1}{Pr}\end{array}$$

(12)

where

$$\begin{array}{c}\hfill Ra=\frac{g\beta S_{\lambda }{H}^5}{\alpha \nu K},Ec=\frac{k{\nu }^2}{S_{\lambda }{C}_P{H}^4},Pr=\frac{\nu }{\alpha },\varphi =2\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2\right]+{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2\end{array}$$

are, respectively, the modified Rayleigh number and modified Eckert number because of the interaction of buoyancy and IHG, the Prandtl number, and the dimensionless form of the viscous dissipation function. $Ra$ and $Ec$ depend upon the rate of internal heat generation. The last term in Equation (12) is proposed to have the following form

$$\begin{array}{c}\hfill \lambda =\frac{Q_{\lambda }}{Pr}\end{array}$$

where $Q_{\lambda }$ is constant volumetric heat source term. Here, $Q_{\lambda }=0$ corresponds to the absence of IHG, whereas $Q_{\lambda }=1$ represents the presence of a constant volumetric heat generation rate. Further, by defining $A=\frac{L}{H}$, the initial and boundary conditions in dimensionless form now can be written as:

At $t<0$

$$\begin{array}{c}\hfill u=v=\psi =\omega =\theta =00\le y\le 10\le x\le A\end{array}$$

At $t\ge 0$

$$\begin{array}{c}\hfill u=v=\psi =0,\omega =\frac{\partial v}{\partial x}=\frac{-{\partial }^2\psi }{\partial x^2},\theta ={\theta }_{C}x=00\le y\le 1\end{array}$$

$$\begin{array}{c}\hfill u=v=\psi =0,\omega =\frac{\partial v}{\partial x}=\frac{-{\partial }^2\psi }{\partial x^2},\theta ={\theta }_{H}x=10\le y\le 1\end{array}$$

$$\begin{array}{c}\hfill u=v=\psi =0,\omega =-\frac{\partial u}{\partial y}=-\frac{{\partial }^2\psi }{\partial y^2},\frac{\partial \theta }{\partial y}=0y=00\le x\le A\end{array}$$

$$\begin{array}{c}\hfill u=v=\psi =0,\omega =-\frac{\partial u}{\partial y}=-\frac{{\partial }^2\psi }{\partial y^2},\frac{\partial \theta }{\partial y}=0y=10\le x\le A\end{array}$$

(13)

The motion is generated by the coupling of the term of buoyancy in Equation (11), the existence of constant volumetric heat generation source and viscous dissipation function in Equation (12), whose boundary conditions are represented in Equation (13), whereas the solution of Equation (10) is used to update the components of velocity. The local Nusselt number in non-dimensional form along the heated and cooled walls are, respectively, defined as

$$\begin{array}{c}\hfill N{u}_R=-{\left(\frac{\partial \theta }{\partial n}\right)}_{n=1}N{u}_L={\left(\frac{\partial \theta }{\partial n}\right)}_{n=0}\end{array}$$

(14)

where n refers to the direction perpendicular to the solid boundary wall. Average Nusselt number across the heated and cooled walls are, respectively, defined by the expression:

$$\begin{array}{c}\hfill N{u}_{avR}={\int }_0^{1}N{u}_{R}dyN{u}_{avL}={\int }_0^{1}N{u}_{L}dy\end{array}$$

(15)

The coupled partial differential Equations (10)–(12) are numerically solved in order to satisfy the boundary conditions defined by Equation (13) to study the flow behavior of the proposed model.

3. Method of Solution

The flow is generated by the coupling of the buoyancy term in Equation (11) with solution of Equation (12). The stream function is acquired by coupling of Equations (10) and (11). The dimensionless form of Equation (8) is used to update the components of velocity at each time step. The time-dependent vorticity transport Equation (11) and energy Equation (12) accompanied with non-dimensional boundary conditions shown in Equation (13) are solved by applying Implicit Finite Difference Method (IFDM) and upwind differencing scheme is used for convective terms. The aim is to obtain a steady state solution from the discretized equations of stream, vorticity, and energy equations. Successive over-relaxation technique is applied to obtain the solution of the steady state stream function Equation (10). We employ a tolerance level of ${10}^{-5}$ when solving the stream function equation. Consider H as the cavity reference height, we choose a constant mesh size $h=H/(jmax-1)$, in which $jmax$ represents the maximum number of mesh points that are equally spaced across H. Throughout the computation, the cavity width is assumed to be $X=AH$, where A stands for aspect ratio and $H=1$ for the reference height. For $\mathsf{\Delta }X=\mathsf{\Delta }Y=h$, the parameter of relaxation $\Phi$ is defined by the following Equation [31,32,35]

$$\Phi =2\left(\frac{1-\sqrt{\left(1-\zeta \right)}}{\zeta }\right)$$

(16)

where

$$\zeta ={\left[\frac{cos\left(\frac{\pi }{imax-1}\right)+cos\left(\frac{\pi }{jmax-1}\right)}{2}\right]}^2$$

(17)

Iterative technique on the discretized equation can be described by

$$\begin{array}{c}\hfill {\psi }_{\left(i,j\right)}^{k+1}={\psi }_{\left(i,j\right)}^k+\frac{\Phi }{4}\left[{\psi }_{\left(i+1,j\right)}^k+{\psi }_{\left(i-1,j\right)}^{k+1}+{\psi }_{\left(i,j+1\right)}^k+{\psi }_{\left(i,j-1\right)}^{k+1}-4{\psi }_{\left(i,j\right)}^k-h^2{\omega }_{\left(i,j\right)}^k\right]\end{array}$$

(18)

where $i,j$ represent grid points across x and y directions, respectively, $\omega$ is the vorticity function. Iteration is performed till the difference converges to the tolerance required among two successive k and $k+1$ iterations. Once the stream function is known components of velocity in the computational domain at each time step are updated by using Equation (8). For the time dependent vorticity transport and energy equations, Equations (11) and (12), given the variable values $\theta$ and $\omega$, at any time step in the flow field, we use the IFDM to compute new values of $\omega$ and $\theta$ at the next time step from Equations (11) and (12). Backward time central space discretization is employed for transient, diffusion, and source terms in the IFDM, while, for the convective terms, the IFDM is modified by applying the second order upwind differencing approach. The IFDM used here is first order accurate in time and second order accurate in terms of spatial discretization, and this method is unconditionally stable. Although, as with any other implicit technique, the stability condition for vorticity at implicit boundary walls needs a time step limitation, say $\mathsf{\Delta }t$, in a form like explicit scheme. For an equally spaced grid size, $\mathsf{\Delta }t$ is stated as [31,32,35]

$$\begin{array}{c}\hfill \mathsf{\Delta }t\le \frac{1}{2{\rm Y}\left(\frac{2}{h^2}\right)+\frac{\left|u\right|+\left|v\right|}{h}}\end{array}$$

(19)

where ${\rm Y}$ refers to the coefficient of diffusion terms in transport equations whose solutions are needed. At solid boundary walls, the time step limitation minimises to:

$$\begin{array}{c}\hfill \mathsf{\Delta }t\le \frac{h^2}{4{\rm Y}}\end{array}$$

(20)

For the entire computation in the proposed IFDM, time step size $\mathsf{\Delta }t=1\times {10}^{-7}$ is taken. It was found after the number of numerical runs that appropriate tolerance for temperature to obtain the steady state in this model for the entire computation would be ${10}^{-7}$. So the computations continue until

$$\begin{array}{c}\hfill \left|\frac{{\theta }_{\left(i,j\right)}^{n+1}-{\theta }_{\left(i,j\right)}^n}{{\theta }_{\left(i,j\right)}^n}\right|<{10}^{-7}\end{array}$$

(21)

is satisfied, where n stands for the number of time step and ($i,j$) denotes spatial grid locations across x and y directions, respectively. It is numerically proven that the steady state regime lies well within this range for all flow variables.

A study of grid dependence has been carried out at $Ra=2\times {10}^5$, $Pr=7$, $Ec=1\times {10}^{-6}$, A = 1, and $Q_{\lambda }=1$ to select the suitable grid size, for which results are represented in

Table 1. Grid dependence study of the maximum value of stream function ${\psi }_{max}$ and $N{u}_{avR}$ of the right wall at $Ra=2\times {10}^5$, $Pr=7$, $Ec=1\times {10}^{-6}$, $A=1$ and $Q_{\lambda }$ = 1.

Meshes	${\mathit{\psi }}_{\mathbf{max}}$	Error%	${\mathbf{Nu}}_{\mathbf{avR}}$	Error%	CPU Time
$31\times 31$	2.23302	5.27402	6.83177	5.8888	0:01:06:30
$41\times 41$	2.11525	2.62711	6.42946	3.55628	0:02:38:50
$51\times 51$	2.05968	1.50994	6.20081	2.13246	0:03:23:59
$61\times 61$	2.02858	0.99182	6.06858	1.33078	0:05:05:60
$71\times 71$	2.00846	0.6681	5.98782	0.87293	0:10:36:40
$81\times 81$	1.99504		5.93555		1:05:36:40

. For any X variable, the relative error percentage between the calculated values considering various mesh points can be defined by the following expression

$$\begin{array}{c}\hfill Error\%=\left|\frac{X_{\left(M+10,M+10\right)}-X_{\left(M,M\right)}}{X_{\left(M,M\right)}}\right|\end{array}$$

(22)

where $\left(X_{M,M}\right)$ is earlier calculated variable value for mesh size ($M\times M$). The percentage error has been calculated from the difference between the calculated values of maximum stream function, ${\psi }_{max}$ and average Nusselt number of the right wall, $N{u}_{avR}$ for various choices of mesh points at each time step. Table 1 indicates that the relative error drops to less than 1% between the grid size 71 × 71 and 81 × 81 in terms of computed values of both $N{u}_{avR}$ and ${\psi }_{max}$, whereas for the next grid size, CPU execution time is drastically increased. Results show that when mesh size decreases, the relative error reduces. It can be observed from Table 1 that very small difference between ${\psi }_{max}$ and $N{u}_{avR}$ occurs for grid size 71 × 71 and 81 × 81. Thus, grid size of 71 × 71 has been chosen throughout the computation. The decrease in relative error indicates the grid independence of the solution. The ${\psi }_{max}$ and $N{u}_{avR}$ against time for various choices of mesh size is also shown in Figure 2.

To verify the correctness of developed codes, we have reviewed the study of Churbanov et al. [9], who examined the natural convection of volumetric heat generating fluid in a rectangular enclosure whose four solid boundary walls were kept at a uniform isothermal or adiabatic temperature of $\theta =0$ or $\frac{\partial \theta }{\partial n}=0$ (where n represents the direction perpendicular to the corresponding wall) using FDM in the 2D vorticity-stream function formulation. The results thus obtained are represented in

Table 2. Maximum and minimum values of streamlines and isotherms compared with the literature [9] at $Ra=6.4\times {10}^5$, $Pr=7.0$, $A=1$, $Q_{\lambda }$ = 1, $Ec$ = 0 and walls at isothermal temperature $\theta$ = 0.

	${\mathit{\psi }}_{\mathbf{max}}$	${\mathit{\psi }}_{\mathbf{min}}$	${\mathit{\theta }}_{\mathbf{max}}$	${\mathit{\theta }}_{\mathbf{min}}$
Present	0.8918	−0.8918	0.0451	0
Churbanov et al. [9]	0.878	−0.878	0.0477	0

, as well as in Figure 3 for $Ra=6.4\times {10}^5$, $Pr=7.0$, $A=1$, $Q_{\lambda }=1$, $Ec$ = 0, for $41\times 41$ grid points and using isothermal conditions at $\theta$ = 0. Comparison shown in Table 2 represents a very good result’s agreement with previous work of Churbanov et al. [9].

We have also revisited the work of Orhan and Pop [35]. They investigated the natural convection flow of micropolar fluid in a square enclosure whose left boundary wall was kept at a temperature $\theta =0.5$ higher than the temperature of the right wall at $\theta =-0.5$, while horizontal walls were kept adiabatic by applying the alternate direction implicit method (ADI) along with the iterative SOR technique. Thus, at $Ra={10}^6$, $Pr=0.71$, $Ec=0$, and $A=1$, the value of heat transfer rate given by Orhan and Pop [35] was 8.945, while the value of heat transfer rate acquired by the current scheme is 8.921, as shown in Figure 4. It has been noticed that our results match well with those of Orhan and Pop [35]. On the basis of this successful validation, the proposed problem is numerically solved by applying the code (the effect of IHG and VD).

4. Results and Discussion

Numerical simulations for the natural convection of Newtonian fluid are performed in a differentially heated rectangular enclosure in the existence of IHG and VD. Results for streamlines, isotherms along with average and local heat transfer rate are represented graphically for various values of non-dimensional controlling parameters, e.g., modified Rayleigh number (${10}^3\le Ra\le {10}^7$), cavity aspect ratio ($1\le A\le 4$), modified Eckert number ($0\le Ec\le {10}^{-6}$), the Prandtl number ($Pr=0.7$, $1.0$, $6.2$, 15), and volumetric internal heat generation parameter ($Q_{\lambda }=0$, 1).

4.1. The Effect of Cavity Width

Figure 5 represents the influence of cavity width on streamlines for $Ra=5\times {10}^4$, $Pr=0.7$, $Ec={10}^{-6}$, and $Q_{\lambda }=1$. The aspect ratio effect is investigated by considering four different values, A = $1,2,3$, and 4. A single anticlockwise rotating cell covering the entire cavity is developed in each case. Recirculating vortices are generated, which are shown by closed streamlines because of the buoyant effects of the difference in temperature between the right hot wall and the left cold wall, and also accompanied by the impact of constant volumetric IHG and VD.

Figure 5b–d indicate a slight distortion of the streamline pattern in the core region as the aspect ratio increases. Flow intensity is analyzed through the maximum stream function value, ${\psi }_{max}$. Flow strength increases due to an increase in aspect ratio. Increasing A creates sufficient room for natural convection to take place, which develops due to buoyancy forces in general and internal heat generation in particular, enabling a greater number of particles to contribute due to internal heat generation along with viscous dissipation. Thus, with the increase in A, the maximum value of stream function, ${\psi }_{max}$, increases from $12.2571$ to $13.7879$. Due to the increasing contact area between the heated fluid and the surrounding cooled fluid, heat convection in the cavity is increased with an increased A for fixed $Ra$.

Figure 6 depicts the isotherms at $Ra=5\times {10}^4$, $Pr=0.7$, $Ec={10}^{-6}$, and $Q_{\lambda }=1$. For small aspect ratio, the isotherms cluster along the walls, particularly along the top portion of the cold wall and the bottom portion of the hot wall, suggesting steep temperature gradients exist in the normal direction of the vertical walls. The isotherms also represent the thermal stratification in the cavity center. For increasing values of A, the thermal stratification between layers in the center of the enclosure is closely packed, suggesting that a large temperature gradient exists in the middle of the enclosure in the vertical direction. For $A>1$, there is a wide distance between the cooled wall and the heated wall, which causes a lesser exchange of heat between cold and hot fluid in the middle of the enclosure before it reaches the opposing wall. However, a shorter distance between the cooled vertical wall and the heated cavity wall at a low aspect ratio ($A=1$) is more likely to result in greater interchange of heat at the thermally active side walls rather than between fluid molecules in the center of the enclosure. For low A, the temperature contours are, therefore, concentrated across the thermally active sides. Also, the thermal boundary layer rises with the increase in A.

The $N{u}_{avR}$ is measured by the Nusselt number along the heated vertical side. Consequently, Figure 7a depicts that average heat transfer $N{u}_{avR}$ decreases, and Figure 7b shows how local heat transfer $N{u}_R$ decreases along the heated vertical wall as aspect ratio increases.

4.2. Effect of Prandtl Number ($Pr$)

The steady state streamlines and isotherm pattern at $Ra=2\times {10}^4$, $Ec={10}^{-6}$, $A=2$ and $Q_{\lambda }$ = 1 for various values of $Pr=0.7$, $1.0$, $6.2$, and 15 are shown in Figure 8 and Figure 9. Figure 8a shows a very high vortex cell strength at $Pr=0.7$. The immense recirculation of vortex cell suggests that at this value of $Pr$, components of velocity are also very high. The flow strength decreases as the value of $Pr$ increases. The buoyancy effects decrease with the increase in viscous effects, which ultimately causes a decline in the maximum value of velocity components and stream function. The vortex cell strength drops from 9 to 0.55 as $Pr$ is increased from $Pr=0.7$ to 15.

Figure 9 shows the steady state isotherm pattern for $Pr=0.7$, $1.0$, $6.2$, and 15, respectively. It is observed that the isotherms look qualitatively similar, suggesting that the isotherms do not change significantly with a smaller $Pr$. However, for fluids with high $Pr$, the thermal boundary layer is small. It can be seen from a comparison of Figure 9a,d that the thermal boundary layer decreases as $Pr$ increases. Fluids with high $Pr$ have greater viscosity, due to which thermal boundary layer thickness decreases.

Figure 10a–c represent the impact of $Pr$ on average heat transfer rates, local heat transfer, and the maximum value of the stream function, respectively. The average and local heat transfer from the heated side wall increases with the rise in $Pr$ from $0.7$ to 15. This was expected since increasing the $Pr$ would increase the viscous effects, thereby causing a decrease in velocities and maximum value of stream function, and a significant decrease in the thermal boundary layer.

4.3. Effect of Rayleigh Number (When $Ec=0$ )

Figure 11 demonstrates the steady state streamlines pattern at $A=2$, $Pr=6.2$, $Ec=0$, and $Q_{\lambda }=1$ for various values of $Ra={10}^3$, ${10}^4$, ${10}^5$, ${10}^6$, and ${10}^7$, respectively. The flow fields are represented by streamlines. At $Ra={10}^3$, the streamlines form almost centrally located single cell. Recirculation in the cavity is very low, indicating that viscous effects are dominant over buoyancy effects. We note that the central vortices are distorted as the $Ra$ increases and the flow cell move towards the cold wall. The fluid molecules, when more heated near the hot vertical wall, are forced to move towards the cold wall because of the buoyancy forces. With the increase in $Ra$, streamlines become densely packed near the walls.

The maximum absolute value of the stream function ($|{\psi }_{max}|$) can be regarded as a measure of the strength of natural convection within the enclosure. It is obvious from Figure 11 that $|{\psi }_{max}|$ significantly rises from $0.18$ to $7.37$ by increasing the $Ra$ from ${10}^3$ to ${10}^7$ due to increase in buoyancy driven convection, which shows faster flow as natural convection intensifies. The strong recirculation of vortex cell in Figure 11e is, therefore, due to the dominant convective forces over the viscous forces.

The thermal fields are shown as isotherms in Figure 12 for various values of $Ra$. At $Ra={10}^3$, the isotherms display the features of low convection, forming a diagonally symmetrical structure. The isotherms are approximately parallel to the vertical boundary walls. This is because of the weak effect of convective flow. However, when $Ra$ increases, the isotherms change their structure, obviously due to the strong influence of buoyancy forces. The pattern of isotherms shows that as $Ra$ increases, the clustering of isotherms near the active walls becomes more pronounced and the temperature gradients near the vertical boundaries become more severe. At $Ra={10}^5$, a thermal boundary layer is established along the heated and cooled side walls. As $Ra$ increases from ${10}^6$ to ${10}^7$, the thermal boundary layer across the side boundary walls becomes denser. This is because convection is an eminent heat transfer mechanism at high $Ra$. As the Ra increases, the rate of heat transfer increases, the velocity increases, and the isotherms are not perpendicular any more and are deformed until they become normal to the vertical boundaries and parallel to the horizontal boundaries. Convection heat transfer changes the temperature distribution to such a level that the temperature gradients in the middle region approach zero. It is clear from the isotherm pattern that at $Ra={10}^7$, there is downward flow near the cold side and strong upward flow at the heated side wall.

Figure 13 depicts the influence of $Ra$ on average and local heat transfer with $A=2$, $Pr=6.2$, $Ec=0$, and $Q_{\lambda }=1$. The $N{u}_R$ and $N{u}_{avR}$ at the heated wall enhance significantly with the increase in $Ra$ from ${10}^3$ to ${10}^7$. The increase in $Ra$ means an increase in the fluid region’s temperature difference, which means an increase in the buoyancy forces that increase the natural convection, and, in turn, the average and local Nusselt numbers increase, as shown in Figure 13.

4.4. Effect of Rayleigh Number in the Presence of VD (When $Ec={10}^{-5}$ )

Figure 14 and Figure 15 demonstrate the impact of heat on streamlines and isotherm pattern due to viscous dissipation by taking the value of $Ec$ as ${10}^{-5}$ for $Ra={10}^3$, ${10}^4$, ${10}^5$, ${10}^6$, and ${10}^7$ at $A=2$, $Pr=6.2$, and $Q_{\lambda }=1$. In this case, the motion is generated by the coupling of buoyancy forces, viscous dissipation, and constant volumetric heat generation. The effect of the VD on the convective flow strength is more pronounced for higher $Ra$ values. The addition of heat because of VD, although small, will increase the temperature of the fluid internally, which will help to accelerate the fluid’s motion. With the increase in $Ra$ from ${10}^3$ to ${10}^6$, the convective flow strength rises from 0.2 to 3.88. The steady flow structure remains until $Ra={10}^6$, over which the flow becomes oscillatory for $Ra={10}^7$, as shown in Figure 14e and Figure 15e. The flow is oscillatory for $Ra>{10}^6$ and does not regain its steady state. This oscillatory effect is visible in Figure 16b, which shows a magnified curve for heat transfer on a smaller temperature and time scale. i.e., from $t=0.5$ to $t=1$, and $N{u}_{avR}$ from −0.5 to 1 is shown. This effect was obvious in Figure 16a due to the uniform scale. It can also be seen that at this value of $Ra$, heat transfer is negative due to the combined effects of internal heat generation, viscous dissipation, and high buoyancy; however, in the laminar flow regime, heat transfer increases as $Ra$ increases.

4.5. Effect of Volumetric Internal Heat Generation Parameter ($Q_{\lambda }$)

The influence of the IHG parameter ($Q_{\lambda }=0$, $Q_{\lambda }=1$) on the steady state pattern of flowfields while $Ra=5\times {10}^6$, $Ec=0$, $Pr=0.7$, and $A=2$ is shown in Figure 17. It is evident that in the presence of the parameter of IHG, ${\psi }_{max}=44.189$, and in the absence of the parameter of heat generation, ${\psi }_{max}=42.77$. This is because buoyancy forces increase in the former case, causing an increase in the flow rate. Consequently, in the case of IHG, the velocity distribution is higher. This influence of IHG on the streamlines is justified as it supports the forces of buoyancy by speeding up the fluid flow.

Due to the effect of heat generation, the fluid temperature rises significantly, as shown in Figure 18 for $Ra=5\times {10}^6$, $Ec=0$, $Pr=0.7$, and $A=2$. The thermal state of the fluid increases due to the existence of heat generation effect, which causes the thermal boundary layer to increase. At the core of the enclosure, the isotherms are horizontal, indicating a stratified flow.

The influence of the IHG parameter on $N{u}_{avR}$ and $N{u}_R$ from the heated wall while $Ra=5\times {10}^6$, $Ec=0$, $Pr=0.7$, and $A=2$ is shown in Figure 19. We observe that rates of heat transfer from the heated wall decrease as the $Q_{\lambda }$ increases. This is expected as the mechanism for heat generation raises the temperature of the fluid near the heated wall, resulting in lesser heat exchange due to the reduction in the temperature difference.

4.6. Effect of Eckert Number ($Ec$)

Figure 20 depicts the impact of the VD parameter on the streamlines and isotherms at $Ra={10}^4$, $Q_{\lambda }=0$, and $Pr=0.71$ for different values of $Ec=0$, ${10}^{-6}$, ${10}^{-5}$, and ${10}^{-4}$, respectively. In the case of $Ec=0$, Figure 20a clearly shows that the center is expanded by two circulations in the center. The influence of the $Ec$ on the streamlines appears to be insignificant for $Ec={10}^{-6}$ and $Ec={10}^{-5}$. The most significant influence of the existence of VD can be observed in the case of $Ec={10}^{-4}$, where there is high flow velocity. It can be seen that the existence of VD tends to the stream lines towards the cold vertical wall. The vortex cell’s strength rises from $6.539$ to $6.933$ as the $Ec$ number increases from $Ec=0$ to ${10}^{-4}$. The steady flow structure remains until $Ec={10}^{-4}$, beyond which the flow becomes oscillatory for $Ec={10}^{-3}$ and ${10}^{-2}$.

Figure 21 shows the impact of heat caused by VD on the isotherms. VD tends to elevate the temperature in the flow region. A comparison between Figure 21a,d reveals that the thermal boundary layer rises as $Ec$ rises from 0 to ${10}^{-4}$. The variations of $N{u}_{avR}$ and $N{u}_R$ over time at $Ra={10}^4$, $Q_{\lambda }=0$ and $Pr=0.71$ for different $Ec$ values are depicted in Figure 22. Due to the frictional or drag forces, the energy of heat is retained in the fluid as $Ec$ increases. Consequently, the fluid temperature rises, reducing the temperature variations close to the heated wall and resulting in decreased average and local heat transfer rates, as shown in Figure 22a,b. From Figure 22c,d, it can be observed that the oscillations in $N{u}_{avR}$ begin to appear for $Ec>{10}^{-4}$.

In order to have deeper understanding of this oscillatory effect, we have also plotted the heat transfer results for $Ra={10}^5$ in Figure 23. For greater values of Rayleigh number ($Ra={10}^5$), the oscillatory behavior of flow starts to appear earlier with increasing Eckert number compared to ($Ra={10}^4$). In this case, oscillatory flow starts to appear at $Ec={10}^{-4}$, as depicted in Figure 23c. The influence of $Ec$ on $N{u}_{avR}$ and $N{u}_R$ at $Ra={10}^5$, $Q_{\lambda }=0$, and $Pr=0.71$ is shown in Figure 23. It is now obvious that the average and local Nusselt number decrease with an increase in $Ec$.

The impact of VD on the rate of heat transfer across the heated wall is plotted in Figure 24 for three different $Ec$ values and for A values of 1 and 2 as a function of $Ra$ at $Pr=0.7$ and $Q_{\lambda }=0$. When there is no viscous dissipation effect, heat transfer rates rise with the increase in $Ra$. Since viscous dissipation is a measure of friction between the adjacent fluid layers, the rise in $Ec$ contributes to frictional heating, causing the fluid temperature to rise. Consequently, the heat transfer rate decreases with an increase in $Ec$. It can be observed from Figure 24b that for $Ec=5\times {10}^{-6}$, heat transfer rates become negative earlier compared to a smaller value of viscous dissipation parameter ($Ec={10}^{-6}$) due to increased frictional heating. For $Ec={10}^{-6}$, the average rates of heat transfer increase up to $Ra$ = ${10}^7$, then decrease further and remain negative with the increase in $Ra$. In contrast, for $Ec=0$, the $N{u}_{avR}$ continuously increases up to $Ra=3\times {10}^7$. The average heat transfer rates decrease for $Ec=5\times {10}^{-6}$ up to $Ra=2424430$, where $N{u}_{avR}=-1.93739$, then increases up to $Ra=4\times {10}^6$, where $N{u}_{avR}=2.68$. Beyond that $Ra$ value, the average heat transfer rates fall slowly and finally remain negative due to oscillatory nature of the flow. Also, as the aspect ratio increases, there are more molecules available to dissipate heat. Therefore, at $A=2$ heat transfer rates are also negative for $Ec={10}^{-6}$, whereas they are not negative at $A=1$.

Figure 25 focuses on the influence of $Ec$ on average rates of heat transfer in the existence of volumetric IHG. Upon careful observation of the graph, it is noticeable that in the presence of the IHG parameter, heat transfer rates are lower compared to the case when $Q_{\lambda }=0$. For instance, at $Ra=2\times {10}^7$ and A = 2, heat transfer rates of $Ec=0$ and $Ec={10}^{-6}$ are $24.20$ and $3.87$, when $Q_{\lambda }$ = 0. However, its value is $23.4$, $2.86$ for $Ec=0$ and $Ec={10}^{-6}$, respectively, in the presence of constant IHG.

5. Conclusions

We have investigated numerically the influence of constant volumetric IHG and VD on the natural convection flow in a differentially heated rectangular enclosure using the IFDM technique and with the buoyancy effects treated using the Boussinesq approximation. The effects of modified Rayleigh number, cavity aspect ratio, Eckert number, Prandtl number, and constant IHG parameter were investigated in terms of isotherms, streamlines, and heat transfer rates in the enclosure. The findings derived from the numerical examination result in the subsequent deductions.

• As the width of the cavity increases from 1 to 4, the intensity of fluid circulation increases. However, the average and local heat transfer decreases along the heated vertical wall as the aspect ratio increases.
• The flow strength and thermal boundary layer decrease as the value of Pr increases. The average and local heat transfer show an increase as Pr rises from 0.7 to 15.
• The thermal performance and flow field strengths within the enclosure amplify as the Rayleigh number increases. Also, the average heat transfer enhances significantly with the increase in Ra.
• The steady flow structure persists until $Ra={10}^6$, beyond which the flow transitions into an oscillatory state at $Ra={10}^7$ in the presence of viscous dissipation. In this scenario, the motion is driven by the combined effects of buoyancy forces, viscous dissipation, and constant volumetric heat generation.
• The average Nusselt number decreases with the increase in the volumetric internal heat generation parameter ($Q_{\lambda }$) and viscous dissipation.
• The strength of the vortex cell increases with the rise of the Eckert number from $Ec$ = 0 to ${10}^{-4}$. The stability range observed at $Pr=0.71$, $A=2$ and, $Ra={10}^4$ is $0\le Ec\le {10}^{-4}$. However, for $Ec>{10}^{-4}$, the flow exhibits an oscillatory pattern.
• A heat transfer comparison is also made for various values of VD parameter ($Ec=0$, ${10}^{-6}$, $5\times {10}^{-6}$) as a function of Rayleigh number, and it is observed that heat transfer rates decreases with the rise in viscous dissipation parameter.
• Our research extensively investigated the simultaneous impact of viscous dissipation, cavity width, and internal heat generation on fluid flow and temperature distribution. While numerous researchers have focused on studying natural convection in cavities either without or with only viscous dissipation [20,36,37].

This study can further be extended for non-zero mach number values such as subsonic flows.