Heatline Analysis for Natural Convection within Porous Rhombic

ARTICLE pubs.acs.org/IECR

Heatline Analysis for Natural Convection within Porous Rhombic Cavities with Isothermal/Nonisothermal Hot Bottom Wall R. Anandalakshmi† and Tanmay Basak* Department of Chemical Engineering, Indian Institute of Technology Madras, Chennai 600036, India †E-mail: [email protected]. ABSTRACT: Analysis has been carried out for energy distribution and thermal mixing in steady laminar natural convective flow through the porous rhombic cavities with inclination angle j for various applications such as geothermal, grain storage, electronic cooling, etc. A generalized non-Darcy model without the Forchheimer inertia term is used to predict the flow in porous media. The effect of the Darcy number (Da) and the role of j on the energy distribution and thermal mixing within porous rhombic cavities with isothermal (case 1) and nonisothermal (case 2) hot bottom walls are illustrated via “heatlines”. Heat transfer is found to be primarily conduction dominant at Da = 10
5 even at a higher Rayleigh number (Ra = 106). The onset of convection occurs at Da = 10
4, and the distorted heatlines from the hot bottom wall take a longer path to reach the cold side walls of the cavity. Larger heat transfer and thermal mixing occurs for Da = 10
3 at Ra = 106 irrespective of j and Pr. Multiple flow/convective circulations are observed at Pr = 0.015 for all j values at Da = 10
3. On the other hand, two asymmetric flow circulation cells are found to occupy the entire cavity at Pr = 0.7, 7.2, and 1000 for j = 75
at Da = 10
3. The cavity with inclination angle j = 30
enhances the convective heat transfer from the hot wall to the cold wall, and the heat transfer to the right cold wall is a maximum for j = 75
, as depicted by “heatlines” irrespective of Pr at Da = 10
3. Average Nusselt number studies based on heatfunction gradients also show that the cavity with j = 30
gives a maximum heat-transfer rate from the bottom to the left wall irrespective of Pr in case 1 at Da = 10
3. The cup mixing temperature (Θcup) is higher for case 1 compared to case 2, and it is almost invariant with j for higher Pr (Pr = 7.2 and 1000) in case 1 at Da = 10
3.

1. INTRODUCTION Natural-convection heat-transfer processes frequently occur in practical applications involving electronics cooling, solar heating, nuclear reactors, heat exchangers, trombe walls in the building industry, etc., because of its effective temperature control in many systems/processes with low cost and high reliability. Recent interest in convective heat transfer through fluid-saturated porous media has been mainly motivated by their importance for many natural and industrial problems, especially those associated with fluid flow in geothermal reservoirs, dispersion of chemical contaminants through water-saturated soil, petroleum extraction, migration of moisture in grain storage systems, energy-efficient drying processes, supply of reactant gases through porous membranes in fuel cells, heat and mass transfer in packed beds and fluidized beds, flow through turbo machines, separation processes in chemical industries, flow through tissues of the human body, building thermal insulation, solid matrix heat exchangers, ceramic processing, iron blast furnace, alloy solidification, energy-efficient drying, catalytic reactions, etc. A comprehensive review of these applications on the subject of porous media may be found in recent books.1
8 A significant amount of literature is available on the convection patterns in enclosures based on various approaches to problems such as the type of fluid, nature of fluid flow, geometry characteristics, orientation of the enclosure, etc. In view of the various industrially important applications,9
17 a fundamental understanding of transport phenomena in porous media is very much essential. Natural convection in fluid-saturated porous media has received considerable attention in recent past because many industrially important processes are governed by this phenomenon.18
23 Phenomena such as porous media, chemical r 2011 American Chemical Society

reaction, etc., may be considered in combination with natural convection to enhance energy transport in the heat-exchange systems. In view of the various applications of energy-efficient processes like microchannel flow, nuclear waste disposal, electronic cooling, material processing, and preservation of food materials and drugs, recent investigators24
29 have carried out extensive numerical investigations on natural convection with porous media in nonrectangular enclosures to analyze the fluid flow and heat transfer. Varol et al.24 performed numerical analysis to study the conjugate natural convection in a triangular enclosure filled with a fluid-saturated porous medium, and they presented flow patterns, a temperature distribution, and heat transfer for various control parameters such as the bottom wall thickness, Rayleigh number, and thermal conductivity ratio. Ahmed et al.25 investigated numerically the heat-transfer behavior inside a porous medium fixed within a complex geometry with the cone and cylinder merged together using a finite element method. Their results indicate that the cone angle of the conical cylinder has a strong influence on the average Nusselt number and the average Nusselt number increases with an increase in the cone angle. Double-diffusive natural convection in a trapezoidal porous cavity with a transverse magnetic field has been studied and validated numerically by Younsi.26 It is found that convective Special Issue: Nigam Issue Received: April 13, 2011 Accepted: July 28, 2011 Revised: July 27, 2011 Published: September 06, 2011 2113

dx.doi.org/10.1021/ie2007856 | Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research heat and mass transfer are strongly inhibited with increasing magnetic field. A stability analysis determining the onset of oscillatory convection in a bound cylindrical cavity containing a viscoelastic fluid-saturated porous medium was performed by Zhang et al.27 Baytas and Pop28 investigated the flow and heattransfer characteristics for a wide range of values of the Rayleigh number, inclined angle, and cavity aspect ratio for steady-state free convection within an inclined cavity filled with a fluidsaturated porous medium. Flow visualization and heat-transfer experiments are reported by Prasad et al.29 for buoyancy-induced flow in a cylindrical cavity filled with a fluid-superposed packed bed and heated from below. Thermal processing of various fluids in porous enclosures has several other important applications.30
39 A few of the numerical investigations in porous enclosures for various applications have been discussed as follows. A physical/mathematical model has been developed by Fedorov and Viskanta36 to simulate the conjugate heat transfer in an actively cooled electronic package. de Lemos and Fischer37 have carried out numerical simulation of a jet impinging against a flat plane covered with a layer of a porous material. A numerical study of the dynamics of grain storage in cylindrical silos was performed by Carrera-Rodriguez et al.38 Das and Lewis39 studied the hydrodynamics of coupled free and porous domains with various structures of porous media for problems encountered in many hydroenvironmental conditions. In view of numerous applications, additional efforts are still needed to understand convection patterns in various shapes of enclosures filled with porous media in order to ensure efficient thermal processing of materials. Therefore, the study of the fluid flow and temperature distribution during natural convection in various types of geometries is necessary, and such studies are yet to appear in the literature. The motivation for this work arises to address issues such as thermal mixing and maximum heattransfer rate during natural convection in rhombic cavities filled with porous media, in the context of material processing and electronics cooling applications. Natural convection in porous rhombic cavities is not studied extensively except in an earlier work by Moukalled and Darwish.40 Numerical solutions are presented by Moukalled and Darwish40 for laminar natural convection heat transfer in a fluid-saturated porous enclosure between two isothermal concentric cylinders of rhombic cross sections for various parameters such as Raleigh numbers, Darcy numbers, porosities, values of enclosure gaps, and Prandtl numbers. The current work attempts for the first time to analyze natural convection within rhombic porous cavities with various inclination angles of side walls based on a heatline approach to visualize heat flow and to analyze the optimal thermal mixing and temperature distribution. The heatline approach was first developed by Kimura and Bejan41 to visualize convective heat transfer. The information of heat flow may be very useful for controlling the temperature distribution in the cavity and also as heatfunctions, which are, in turn, related to the Nusselt number based on the proper dimensionless form. The concept of heatlines has been employed and extended by several researchers to describe various physical phenomena.42
50 Few studies on the analysis of natural convection in porous media based on the heatline concept have also been reported.51
55 The objective of this article is to analyze the thermal energy and fluid flow through the porous rhombic enclosures of the insulated top wall and isothermally cooled side walls with isothermally and nonisothermally heated bottom walls for various practical applications. There are numerous applications for

ARTICLE

Figure 1. Schematic diagram of the physical system for case 1 (isothermally heated bottom wall) and case 2 (nonisothermally heated bottom wall) for inclination angles j = 30
, 45
, and 75
.

isothermal heating. One such application is in the field of effective cooling of electronic components, which are treated as heat sources embedded in the bottom surfaces of the cavity.56 Other potential applications of isothermal heating are metal processing and the efficient design of energy systems.57,58 Nonisothermal boundary conditions could be seen in metallurgy to obtain regular melting for metal, and it can be obtained with a cylindrical heater for laboratory experiments.59,60 A generalized non-Darcy model without the Forchheimer inertia term is used to predict the flow in a porous medium. This generalized model based on volume-averaging principles was developed by Vafai and Tien.61 Numerical investigations on the heating characteristics within porous rhombic cavities have been carried out based on coupled partial differential equations of momentum and energy, which are solved using the Galerkin finite-element method with a penalty parameter to obtain numerical simulations in terms of streamfunctions, heatfunctions, and isotherms. The streamfunctions and heatfunctions of the Poisson equation are also solved using the Galerkin finiteelement method. The jump discontinuity in the Dirichlet type of wall boundary conditions for temperature at the corner points is tackled via implementation of exact boundary conditions at those singular points as mentioned earlier.62 A number of model fluids for almost all industrial application ranges (0.015 e Pr e 1000) have been considered for current analysis. The heat-transfer characteristics are studied with the help of local and average Nusselt numbers, and also thermal mixing for all j values are compared using a cup mixing temperature in both heating situations. Qualitative trends of local and average Nusselt numbers with the Darcy number are explained based on heatfunction gradients.

2. MATHEMATICAL MODELING AND SIMULATIONS 2.1. Governing Equations and Boundary Conditions. The physical domain is shown in Figure 1. The thermophysical properties of the fluid in the flow field are assumed to be constant except the density in the buoyancy term, and the change in density due to temperature variation is calculated using Boussinesq’s approximation. It may be noted that the local thermal equilibrium is valid; i.e., the temperature of the fluid phase is equal to the temperature of the solid phase within the porous medium, and similar approximations were also used by earlier researchers.2 The momentum transport in a porous medium is 2114

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

based on the generalized non-Darcy model. However, the velocity square term or Forchheimer term, which models the inertia effect, is neglected here in the present case because this work deals with natural convection flow within a porous cavity. This inertia effect is more important for non-Darcy effects on the convective boundary-layer flow over the surface of the body embedded in a high-porosity medium. Under these assumptions and following the methodology of earlier works61 with the Forchheimer inertia term neglected, the governing equations for steady two-dimensional natural convection flow in a porous rhombic cavity for the conservation of mass, momentum, and energy may be written with the following dimensionless variables or numbers: x X ¼ , L

y Y ¼ , L

U ¼

uL , R

vL , R

V ¼

θ¼

T
Tc Th
Tc

pL2 ν K , Pr ¼ , Da ¼ 2 , 2 R L FR gβðTh
Tc ÞL3 Pr Ra ¼ ν2 P¼

ð1Þ

The governing equations in dimensionless forms are ∂U ∂V þ ¼0 ∂X ∂Y U

ð2Þ

∂U ∂U ∂P þ V ¼
∂X ∂Y ∂X ∂2 U ∂2 U þ Pr þ ∂X 2 ∂Y 2

U

∂V ∂V ∂P þ V ¼
∂X ∂Y ∂Y ∂2 V ∂2 V þ Pr þ 2 ∂X ∂Y 2

!


Pr U Da

!

and

∂θ ∂θ ∂2 θ ∂2 θ þ V ¼ 2 þ 2 ∂X ∂Y ∂X ∂Y

U

V ¼ 0, V ¼ 0,

∂2 V ∂2 V þ Pr þ 2 ∂X ∂Y 2

ð5Þ

θ ¼ 1 or θ ¼ sinðπXÞ,

for Y ¼ 0, 0 e X e 1 on AB U ¼ 0,

  ∂V ∂V ∂ ∂U ∂V þ V ¼γ þ ∂X ∂Y ∂Y ∂X ∂Y

ð4Þ

and the governing equations (eqs 3
5) are subject to the following boundary conditions: U ¼ 0,

The continuity equation (eq 2) is automatically satisfied for large values of γ. A typical value of γ that yields consistent solutions is 107. Using eq 7, the momentum balance equations (eqs 3 and 4) reduce to   ∂U ∂U ∂ ∂U ∂V þ V ¼γ þ U ∂X ∂Y ∂X ∂X ∂Y ! ∂2 U ∂2 U Pr U ð8Þ þ Pr þ
∂X 2 ∂Y 2 Da

Pr V
Da

þ Ra  Prθ U

ð3Þ

T denotes the temperature, ν and R are kinematic viscosity and thermal diffusivity, respectively, K is the medium permeability, p is the pressure, F is the density, Th and Tc are the temperatures at hot and cold walls, respectively, g is the acceleration due to gravity, β is the volume expansion coefficient, and L is each side of the rhombic cavity. Note that Ra, Pr, and Da are the Rayleigh, Prandtl, and Darcy numbers, respectively. X and Y are the dimensionless distances measured along the horizontal and vertical directions, respectively, U and V are the dimensionless velocity components in the X and Y directions, respectively, P is the dimensionless pressure, θ is the dimensionless temperature, and j is the inclination angle with the positive direction of the X axis. 2.2. Solution Procedure. The momentum and energy balance equations (eqs 3
5) are solved using the Galerkin finite-element method. The continuity equation (eq 2) 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 3 and 4, we use the penalty finite-element method, where the pressure P is eliminated by a penalty parameter γ, and the incompressibility criteria given by eq 2, which results in   ∂U ∂V þ ð7Þ P ¼
γ ∂X ∂Y

V ¼ 0,

U ¼ 0,

V ¼ 0,

N

θ ¼ 0, for X sinðjÞ
Y cosðjÞ ¼ sinðjÞ, ∂θ ¼ 0, for Y ¼ sinðjÞ, ∂Y

cosðjÞ e X e 1 þ cosðjÞ on DC

ð6Þ

Note that, in eqs 1
6, x and y are the distances measured along the horizontal and vertical directions, respectively, u and v are the velocity components in the x and y directions, respectively,

Pr V Da ð9Þ

The system of equations (5), (8), and (9) with boundary conditions is solved by using the Galerkin finite-element method.63 Expanding the velocity components (U, V) and temperature (θ) using the basis set {Φk}N k=1 yields

θ ¼ 0, for X sinðjÞ
Y cosðjÞ ¼ 0,

0 e Y e sinðjÞ on BC




þ Ra  Prθ

U≈

∑ UkΦkðX, Y Þ, k¼1

θ≈

∑ θk Φk ðX, Y Þ k¼1

0 e Y e sinðjÞ on AD U ¼ 0,

!

N

V≈

N

∑ VkΦkðX, Y Þ, k¼1

and ð10Þ

The Galerkin finite-element method yields the nonlinear residual equations for eqs 8, 9, and 5 at nodes of internal domain Ω. The detailed solution procedure is given in earlier works. 62,64 2.3. Post-Processing: Streamfunction, Nusselt Number, and Heatfunction. 2.3.1. Streamfunction. The fluid motion is displayed using the streamfunction (ψ) obtained from velocity 2115

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

components (U and V). The relationships between streamfunction and velocity components for two-dimensional flows are65,66 U ¼

∂ψ ∂Y

and

V ¼


∂ψ ∂X

ð11Þ

asymmetric geometry with respect to the central symmetric line of the cavity. The percentage of error in heat balances within the cavity can be calculated using the average Nusselt numbers at the bottom, left, and right walls as follows:

which yield a single equation ∂ψ ∂ψ ∂U ∂V
þ ¼ ∂X 2 ∂Y 2 ∂Y ∂X 2

2

ε¼ ð12Þ

The no-slip condition is valid at all boundaries because there is no cross-flow. Hence, ψ = 0 is used as residual equations at the nodes for boundaries. Using the above definition of the streamfunction, the positive sign of ψ denotes anticlockwise circulation, and the clockwise circulation is represented by the negative sign of ψ. When the streamfunction (ψ) is expanded using the basis N set {Φk}N k=1 as ψ = ∑k=1ψkΦk(X,Y) and the relationships for U and V from eq 10, the Galerkin finite-element method yields the linear residual equations for eq 12, and the detailed solution procedure to obtain ψ values at each node point are given in earlier works.62,64 2.3.2. Nusselt Number. The heat-transfer coefficient in terms of the local Nusselt number (Nu) is defined by ∂θ Nu ¼
∂n

9

∂Φi ∂Y

Nub ¼

∑ θi i¼1

Nul ¼

  ∂Φi ∂Φi
cos j θi sin j ∂X ∂Y i¼1 9

∑

and Nur ¼


  ∂Φi ∂Φi
cos j θi sin j ∂X ∂Y i¼1 9

∑

ð14Þ

ð15Þ

∂Π ∂θ ¼ Uθ
∂Y ∂X ∂Π ∂θ ¼ Vθ

∂X ∂Y

∂2 Π ∂2 Π ∂ ∂ ðUθÞ
ðV θÞ þ ¼ 2 2 ∂X ∂Y ∂Y ∂X

ð20Þ

Using the above definition of the heatfunction, the positive sign of Π denotes anticlockwise heat flow, and the clockwise heat flow is represented by the negative sign of Π. When the heatfunction (Π) is expanded using the basis set {Φk}N k=1 as Π = ∑N k=1ΠkΦk(X,Y) and the relationship for U, V, and θ from eq 10, the Galerkin finite-element method yields the linear residual equations for eq 20, and the detailed solution procedure to obtain Π values is given in earlier works.62,64 The residual equation of eq 20 is further supplemented with various Dirichlet and Neumann boundary conditions in order to obtain a unique solution of eq 20. Boundary conditions are specified as follows for both cases as derived from eq 19. (a) bottom wall

ð16Þ

n 3 ∇Π ¼ 0

ðfor isothermal heating of the bottom wallÞ

ð21Þ n 3 ∇Π ¼ π cos πX

ðfor sinusoidal heating of the bottom wallÞ

ð22Þ (b) side wall n 3 ∇Π ¼ 0

1

Nul ¼

ð19Þ

which yield a single equation

The average Nusselt numbers at the bottom and inclined side walls are given by Z 1 Nub dX Z 1 Nub ¼ 0 ¼ Nub dX Xj10 0 Z

ð18Þ

2.3.3. Heatfunction. The heat flow within the enclosure is displayed using the heatfunction (Π) obtained from conductive heat fluxes [(
∂θ/∂X), (
∂θ/∂Y)] as well as convective heat fluxes (Uθ, Vθ). The heatfunction satisfies the steady energy balance equation (eq 5)41 such that

ð13Þ

where n denotes the normal direction on a plane. The normal derivative is evaluated by the biquadratic basis set in the ξ
η domain. The local Nusselt numbers at the bottom wall (Nub), at the left wall (Nul), and at the right wall (Nur) are defined as

jNub
ðNul þ Nur Þj  100 min½Nub , ðNul þ Nur Þ

ðfor the isothermal cold side wallsÞ ð23Þ

Nul dS1 0

and

Z

Nur ¼

1

Nur dS2

ð17Þ

0

Note that dS1 and dS2 are the small elemental lengths along the left and right walls, respectively. The average Nusselt numbers are also useful to benchmark the overall heat balance within the cavity. Note that Nul + Nur = Nub, whereas Nul 6¼ Nur because of

The top insulated wall may be represented by Dirichlet boundary conditions as obtained from eq 19, which is simplified into ∂Π/∂X = ∂θ/∂Y = 0 for an adiabatic wall, and a reference value of Π = 0 may be assumed at Y = sin j, which leads to Π = 0 "X at Y = sin j. The boundary conditions at hot
cold junctions are obtained by integrating eq 19 along boundaries from the reference point until the end point of each junction. The derivation of Dirichlet conditions for Π at hot
cold junctions is given below. 2116

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. (a) Mapping of the rhombus domain to a square domain in the ξ
η coordinate system (b) Mapping of an individual element to a single element in the ξ
η coordinate system.

(i) At (X = 0, Y = 0), Z Πð0, 0Þ ¼ Πðcos j, sin jÞ þ

Y ¼0

Y ¼ sin j



 ∂Π dS ¼ Nul ∂S

ð24Þ (ii) At (X = 1, Y = 0), Πð1, 0Þ ¼ Πð1 þ cos j, sin jÞ Z Y ¼0   ∂Π þ dS ¼ Nur ∂S Y ¼ sin j

convective flow exists. The cup mixing temperature (Θcup) and spatial average or area average temperature (Θavg) are given as Z Z V^ ðX, Y Þ θðX, Y Þ dX dY Z Z ð26Þ Θcup ¼ V^ ðX, Y Þ dX dY where V^ ðX, Y Þ ¼

ð25Þ

Because Π = 0 at Y = sin j, "X, Π(cos j, sin j) and Π(1 +cos j, sin j) are taken as zero in an adiabatic top wall. It may be noted that the unique solution of eq 20 is strongly dependent on the nonhomogeneous Dirichlet conditions. It may also be noted that most of the earlier works41,42,44,50 are limited within two adiabatic walls where Dirichlet boundary conditions are either 0 or Nu at the adiabatic walls for fluids confined within square cavities only. To assess the performance of isothermal and nonisothermal heating cases in terms of achieving a maximum heat-transfer rate or heat distribution, the average temperature across the cavity, i.e., the cup mixing temperature, is defined. The cup mixing temperature or the velocity-weighted average temperature (Θcup) is a suitable tool to analyze the intensity of mixing when

and Θavg ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U2 þ V 2

Z Z θðX, Y Þ dX dY Z Z dX dY

ð27Þ

3. RESULTS AND DISCUSSION 3.1. Numerical Procedure and Validation. The computational grid within the rhombus is generated via mapping of the rhombus into a square domain in the ξ
η coordinate system, as shown in Figure 2, and the procedure is outlined in the Appendix. The computational grid in ξ
η coordinates consists of 28  28 biquadratic elements, which correspond to 57  57 grid points. The biquadratic elements with a lesser number of nodes 2117

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Comparison of the Present Results with the Benchmark Resolutions of Earlier Works68
75 for Natural Convection in an Air-Filled (Pr = 0.71) Porous Square Cavity with 28  28 Biquadratic Elements (57  57 Grid Points) at Ram = Ra  Da = 1000 step no.

ref

Nu

1

Walker and Homsy68

12.960

2 3

Bejan69 Gross et al.70

15.800 13.448

4

Manole and Lage71

13.637

5

Goyeau et al.72

13.470

6

Baytas and Pop73

14.060

7

Saeid and Pop74

13.726

8

Saeid and Pop75

13.442

9

present work

14.182

Table 2. Comparison of the Average Nusselt Number (Nub) for Various Grid Systems in Cases 1 (a) and 2 (b) at Ra = 106, Da = 10
3, and Pr = 0.7 with Various Inclination angles (j)a j

24  24

26  26

28  28

Case 1 30


14.06 (0.28%)

14.22 (0.42%)

14.38 (0.42%)

45


12.72 (0.47%)

12.82 (0.31%)

12.93 (0.39%)

75


12.38 (0.81%)

12.43 (0.81%)

12.48 (0.73%)

30
45


7.52 (2.66%) 7.45 (0.95%)

7.37 (0.27%) 7.43 (0.81%)

7.36 (0.27%) 7.41 (0.40%)

75


8.01 (1.52%)

7.98 (1.66%)

7.96 (1.66%)

Case 2

Figure 3. Streamfunction (ψ) and temperature (θ) contours for the hot left wall and cold right wall with adiabatic horizontal walls at Pr = 5, Ra = 106, j = 90
(square cavity), and (a) Da = 10
5, (b) Da = 10
4, and (c) Da = 10
3.67 Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively.

smoothly capture the nonlinear variations of the field variables, which are in contrast with finite difference/finite volume solutions. In order to assess the accuracy of our numerical procedure, we have tested our algorithm based on the grid size (57  57) for a square enclosure with Pr = 5 fluid subject to the hot left wall and cold right wall in the presence of insulated horizontal walls at Da = 10
5, 10
4, and 10
3 as carried out in an earlier work,67 and the results are in good agreement. The results based on current simulation studies are shown in Figure 3. To validate our present numerical approach, benchmark studies have also been carried out for a porous square cavity filled with air (Pr = 0.71), in the presence of the isothermally heated left wall, while the right wall is maintained isothermally cold and the horizontal walls are maintained adiabatic. Table 1 shows a comparison of the average Nusselt number (Nu) based on the current numerical procedure, with the result the same as that obtained by Walker and Homsy,68 Bejan,69 Gross et al.,70 Manole and Lage,71 Goyeau et al.,72 Baytas and Pop73 and Saeid and Pop74,75 at a modified Rayleigh number (Ram = 1000), where Ram = Ra  Da. Note that the results of the Darcy model (employed by earlier works68
75) may be compared with those of the non-Darcy model (used in the present case) only at low Da values.

a

Percentage of error (ε) from the heat balance (eq 18) within the rhombic cavity is shown in brackets.

Detailed computations have been carried out for various values of Pr (Pr = 0.015, 0.7, 7.2, and 1000), inclination angles (j = 30
, 45
, and 75
), and Da (Da = 10
5
10
3). The isothermal situation along side walls and the bottom wall involve jump discontinuities at the hot
cold junctions, which correspond to mathematical singularities for isothermal heating situations. This problem has been resolved based on the average temperature of the hot and cold sections/walls at the junction and keeping the adjacent grid nodes at the respective section/wall temperatures, as discussed in an earlier work.62 In the current investigation, a Gaussian quadrature-based finite-element method provides smooth solutions in the computational domain including the singular points because evaluation of the residuals depends on interior Gauss points. A current solution scheme produces grid invariant results with reasonable percentage error, as shown in Table 2 for cases 1 and 2. The Galerkin finite-element approach as used in the current work also offers special advantage on evaluation of the local Nusselt number at cold side walls and the hot bottom wall because element basis functions have been used to obtain a heat flux whereas finite difference/finite volume based methods involve interpolation functions. In order to validate the heatfunction contours or heatlines, we have carried out simulations for a range of Darcy numbers (Da = 10
3
10
5) and Prandtl numbers (Pr = 0.015
1000). The results on the heatfunctions have been validated based on the 2118

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 1 with Pr = 0.015, Ra = 106, and Da = 10
5 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

grid size (57  57) for a nonporous square enclosure filled with air (Pr = 0.71) subject to the hot left wall and cold right wall in the presence of insulated horizontal walls at Ra = 103, 104, and 105 as carried out in an earlier work,50 and the results are in good agreement. The results based on current simulation studies are not shown here for the brevity of the manuscript. Simulations have been carried out at higher Ra (Ra = 106) in order to show the effect of Da (permeability of porous media) on convection in a pure fluid medium. At low Da (Da = 10
5), the fluid is almost stagnant and heat transfer is conduction-dominant even at high Ra (Ra = 106). It is observed that the heatlines emanate from the hot surface and end on the cold surface, which are perpendicular to isothermal surfaces, similar to heat flux lines, during the conduction-dominant regime. As they approach the adiabatic wall, they slowly bend and become parallel to that surface. Also, the heatlines and isotherms are found to be smooth curves, without any distortion in the presence of dominant conductive transport of heat (see Figure 4). At high Da (Da = 10
3), convection is enhanced in the cavity as the permeability of the porous medium increases or the hydraulic resistance decreases with Da. The magnitudes of the streamfunction increase, and heat flow during enhanced convection may be visualized by heatlines. A large amount of heat is drawn from the hot bottom wall, as indicated by dense heatlines with large magnitudes of Π (see Figure 6). These features are discussed in sections 3.2 and 3.3. The sign of the heatfunction needs special mention. The solution of the heatfunction, a Poisson equation (eq 20), is strongly dependent on the nonhomogeneous Dirichlet boundary conditions (eqs 24 and 25), and the sign of the heatfunction is governed by the sign of the “nonhomogeneous” Dirichlet conditions. In the current situation, a negative sign of heatlines represents the clockwise flow of heat, while a positive sign refers to the anticlockwise flow. This assumption is in accordance with the sign convention for the streamfunction. The streamfunction and heatfunction have identical signs for convective transport. The detailed discussion on heat

Figure 5. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 1 with Pr = 0.015, Ra = 106, and Da = 10
4 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

transport based on heatlines for various cases is presented in the following sections. 3.2. Case 1: Hot Isothermal Bottom Wall. Figures 4
8 show the effects of Da = 10
5
10
3 and Pr = 0.015
1000 on the fluid flow and heat transfer in porous rhombic cavities for various j values (30
, 45
, and 75
) induced by the hot isothermal bottom wall at Ra = 106. Because of buoyancy and an imposed temperature gradient between the hot isothermal bottom wall and cold isothermal side wall, fluid rises from the hot bottom wall and flows down along the cold side walls, forming two asymmetric rolls inside the cavity. Hence, the flow as well as temperature patterns within the cavity for Da = 10
5
10
3 with all Pr values (Figures 4
8) are not symmetric because of buoyancy-induced flow and geometrical asymmetry with respect to the vertical center line. This is further due to various inclinations of the side walls. 2119

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

Figure 6. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 1 with Pr = 0.015, Ra = 106, and Da = 10
3 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

Figure 7. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 1 with Pr = 0.7, Ra = 106, and Da = 10
3 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

Parts a
c of Figure 4 illustrate the distribution of the streamfunction (ψ), heatfunction (Π), and temperature (θ) for Da = 10
5, Ra = 106, and Pr = 0.015. The small values of Da (Da = 10
5) correspond to low permeability and high hydraulic resistance due to a porous matrix. Therefore, the magnitudes of the streamfunctions are small for all j values. Two asymmetric rolls of fluid circulation cells are formed, and weak fluid circulation cells are observed near the left portion of the cavity. As a result, the maximum streamfunction (|ψ|max) is found to be less for the left circulation cells (|ψ|max = 0.05) compared to the right circulation cells (|ψ|max = 0.19) for j = 30
(see Figure 4a). At j = 45
, fluid circulation cells increase longitudinally in size and the strength of the circulation cells is also increased slightly compared to j = 30
, as seen in Figure 4a,b. As j increases from 45
to 75
(Figure 4b,c), the size of the left circulation cells tends to become stronger. Also, both left and right circulation cells grow in size and magnitude, and two asymmetric circulations still occur within the cavity. However, the right circulation cell still dominates. It is interesting to observe that as j increases

ARTICLE

Figure 8. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 1 with Pr = 1000, Ra = 106, and Da = 10
3 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of the streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

the flow strength is enhanced for both the right and left circulation cells, as seen from the magnitudes of the streamfunctions (see Figure 4a
c). The heat-flow distribution for various j values in porous rhombic enclosures may further be explained with heatlines. The second column of Figure 4a
c illustrates the heatline distributions for j = 30
75
. As a consequence of the weak fluid flow at lower Da, the heat flow is also low and that implies smaller magnitudes of the heatfunction throughout the cavity except near the left bottom corner of the cavity. It is observed that, at the bottom portion of the side walls, the heatlines are denser and they have high magnitudes, which indicate large heat flow to the bottom portions, compared to the top portions of the side walls for all j values (see Figure 4a
c). As j increases, the large area of the top central region is maintained at uniform temperature because of increased heat flow near the top portion of the cavity. The stratification of isotherms near the bottom wall is based on almost vertical heatlines and a large gradient in temperature occurring near the bottom wall (see Figure 4a
c). It is interesting to note that the heat flow from the bottom wall to the left wall is larger near the left bottom corner because of the greater degree of inclination of the left wall toward the bottom wall, although the fluid-flow intensity is small near the lower portion of the left wall due to no-slip boundary conditions maintained at the walls. As j increases, magnitudes of the heatfunction (Π) along the left wall decrease because of a decrease in the heat flow from the bottom wall to the left wall (see Figure 4a
c). Therefore, the thickness of the boundary layer near the left wall increases with j. Heatlines in the right portion of the cavity are not parallel circular arcs, as seen in the left portion of the cavity, and they are found to emanate from a zone near a region to the right bottom corner (see Figure 4a
c). It is also interesting to note that dense heatlines emanating from a very small part of right portion of the hot bottom wall to the entire right cold wall indicate only a small amount of heat being received by the right cold wall. The heatfunction (Π) gradients are smaller along the right wall compared to the left wall because of lower thermal gradients along the right wall compared to the left wall (for example, Π varies within 0.05
9.37 and 0.0001
1.03 along the left and right walls, respectively, for j = 30
). Therefore, the boundary layer 2120

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research thickness along the right wall is larger than that along the left wall, as seen from the gray-scale plot of θ in Figure 4a
c. The zone near the top right corner is found to be cooler. It is also interesting to note that the boundary layer thickness near the right wall decreases as j increases because of an increase in the heat flow from the bottom wall to the right wall (Figure 4a
c). Overall, as j increases, the boundary layer thickness near the left wall is increased because of a decrease in the heat flow, whereas the boundary layer thickness near the right wall is decreased because of an increase in the heat flow. In summary, the left cold wall receives more heat from the large portion of the hot bottom wall for smaller inclination angles (j). As j gradually increases, the right cold wall gradually starts receiving a considerable amount of heat, as indicated by heatlines (see Figure 4a
c). The streamline patterns are qualitatively similar for the Da = 10
4 and Da = 10
5 cases for all j values (see Figures 5a
c and 4a
c). However, the intensity of the circulations is enhanced because of the onset of convection, as seen from the magnitudes of the streamfunction values. Distortion of heatlines indicate the onset of convection at Da = 10
4, and the heatlines take longer paths from the hot bottom wall to reach the cold surfaces as seen from the heatline contours (see Figure 5a
c). The heatlines emanate from the right corner of the cavity and end near the lower portion of the right wall. In addition, the heatlines in the region near the right corner of the bottom wall rise higher toward the top wall and thereafter bend toward the side walls because of the convective motion of the fluid. Also, there is a significant increase in the magnitudes of the heatfunctions compared to Da = 10
5 for all j values because of the enhanced heat flow (see Figures 4a
c and 5a
c). At j = 45
, enhanced heat flow occurs based on the heatline arcs throughout the left wall and the top portion of the right cold wall, as seen in Figure 5b. As the inclination angle (j) further increases to 75
, closed-loop heatline cells also appear in the left portion of the cavity. The heatline cells grow further in size and magnitude (|Π|max = 2.5 for the left heatline cell, whereas |Π|max = 2.3 for the right heatline cell) because of enhanced convective heat flow at j = 75
. It is interesting to note that dense heatlines emanating from the central zone of the bottom wall deliver heat to the top portions of the cavity. In addition, closed-loop heatline cells are found to occur near the lower portion of the cold walls, and those obstruct the heat flow to the lower portion of the cold walls. Also, a significant amount of heat is distributed near the top portion of the cavity. However, the heat is not distributed symmetrically along the cold walls, as illustrated from the heatfunction magnitudes. It is also interesting to note that the heat source of the bottom wall now distributes large heat to a larger regime of the cold side walls compared to the previous case with Da = 10
5 for all j values. It is interesting to observe that the top portion of the cavity is maintained at θ = 0.1
0.3, 0.3
0.5, and 0.5
0.6 for j = 30
, 45
, and 75
, respectively, as shown in Figure 5a
c. The top portion of the cavity receives more heat at j = 75
because of intense closed-loop heatline cells (|Π|max g 2) and that region is maintained at larger temperatures compared to 30
and 45
, as shown in the gray-scale plots of Figure 5a
c. On the other hand, the thickness of the boundary layer near the lower portion of the cold walls is larger for 75
compared to j = 30
because of less heat flow in that region and that is further due to intense closedloop heatline cells near those regions because the heatline cells obstruct the heat flow to the cold wall for 75
(see Figure 5a
c). Because of dense heatlines along the center line, isotherms are

ARTICLE

found to be distorted near the central portion of the cavity for all j values unlike in the previous case with Da = 10
5 (see Figures 4a
c and 5a
c). A larger temperature distribution, θ = 0.5
0.9, is also observed at the central regime of the cavity because of enhanced convection because more distorted heatlines are observed at the central regime for j = 75
(see Figure 5c). On the basis of the highly dense heatlines, isotherms are compressed near the left and right bottom corners of the cavity at j = 75
. Thus, the boundary layer thickness near the left and right corners of the bottom wall is lesser than that near the middle portion of the bottom wall at j = 75
(see Figure 5c). Figure 6 displays streamlines, heatlines, and isotherms for Pr = 0.015, Ra = 106, and Da = 10
3. Complex circulations are induced in the cavity at higher Da because hydraulic resistance of the porous medium is greatly reduced, leading to intense convection. Multiple circulations are formed as a result of the inclination angle j for low Pr fluid. A large primary circulation is observed in the central portion of the cavity. Secondary circulation cells are observed near the top right corner and bottom left corner of the cavity. Small tertiary circulation cells are also observed at the lower left corner of the cavity. At j = 75
, the central primary circulation cell grows bigger in size and magnitude. In addition, one or more fluid circulation cells with |ψ|max = 2 are found near the top portion of the cavity in contrast to other j values (see Figure 6a
c). The intensity of the fluid motion is greatly increased for all j values at Da = 10
3 compared to that of Da = 10
5, as indicated by |ψ|max. It may be noted that |ψ|max = 9.82, 12.38, and 13.87 for j = 30
, 45
, and 75
, respectively, occur because of the primary circulation cell at Da = 10
3, whereas |ψ|max = 0.19, 0.24, and 0.27 for j = 30
, 45
, and 75
, respectively, occur within the primary circulation cell at Da = 10
5 (Figures 6a
c and 4a
c). Heatlines are found to be distorted more because the heat distribution is dominated by convection. In contrast to Da = 10
5, much of the heat distribution in the cavity is due to heat sources from the right portion of the bottom wall and that is distributed to the upper portion of the side walls, as indicated by heatlines at higher Da (=10
3) for all j values. On the other hand, a lower amount of heat is drawn from heat sources of the left portion of the bottom wall and that is distributed to the lower portion of the left cold wall at higher Da for all j values (see Figure 6a
c). It is also interesting to observe that the trend on the heat-flow patterns via heatline cells and heat-flow distribution along the right portion of the cavity are qualitatively similar for all j values except that the primary heatline cell grows bigger in size and magnitude with j (see Figure 6a
c). As j increases, the convective heat transport in the central portion of the cavity gradually becomes stronger for Da = 10
3. Primary closed-loop heat circulation cells with |Π|max = 2.5, 3.8, and 5.8 are found near the central portion of the cavity for j = 30
, 45
, and 75
, respectively, and those heatline cells further obstruct the heat distribution near the lower portion of the right cold wall at Da = 10
3. These heatline cells also promote thermal mixing at the central zone for all j values. Localized heat circulation cells occur near the lower corners of the left cold wall because of the presence of a tiny amount of heat circulation in that regime, and they further resist the heat distribution near the lower portion of the left cold wall for all j values (see Figure 6a
c). It is also interesting to note that the boundary layer thickness is small near the corners of the bottom wall for all j values, indicating a large amount of heat being transferred from the 2121

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research bottom corners of the cavity, as seen from dense heatlines near those regions. On the other hand, the larger boundary layer thickness occurs near the middle portion of the bottom wall for all j values, indicating less heat transport from the middle portion of the bottom wall (see Figure 6a
c). The thermal boundary layer thickness near the lower portion of the left wall is larger because of the presence of secondary heatline cells for all j values at Da = 10
3. Because of the presence of sparse heatlines, isotherms are less compressed near the middle portion of the left cold wall for j = 75
, which further corresponds to comparatively larger boundary layer thicknesses than the j = 30
and 45
cases in a large regime along the middle and top portions of the left wall (Figure 6a
c). Stronger primary closed-loop convective heatline cells lead to enhanced thermal mixing especially for j = 75
, and thus a large central portion of the cavity is maintained at θ = 0.3
0.5 (Figure 6c). As a result of high heat transport to the right cold wall, the boundary layer thickness along the right wall is also reduced for all j values compared to Da = 10
5 (see Figures 6a
c and 4a
c). Even though the top right portion of the cavity receives more heat because of enhanced convection at Da = 10
3 compared to Da = 10
5, part of the fluid near the top right portion of the cavity remains cold because of the low intense heat flow, as seen from 0.002 e |Π| e 0.2 (see Figure 6a). Distributions of streamlines, heatlines, and isotherms for Pr = 0.7 and Da = 10
3 are displayed in Figure 7a
c. The intensity of fluid circulation is found to be stronger for Pr = 0.7 compared to Pr = 0.015, and the maximum values of the streamfunctions (|ψ|max) are increased for all j values (Figures 7a
c and 6a
c). Note that |ψ|max are 10.10, 11.57, and 14.61 for j = 30
, 45
, and 75
, respectively, at Pr = 0.7 and Da = 10
3, whereas |ψ|max are 9.82, 12.38, and 13.87 for j = 30
, 45
, and 75
, respectively, at Pr = 0.015 and Da = 10
3. Further, the cells expand, and they are diagonally skewed because of intense convection, unlike in the previous case with Pr = 0.015 for j = 30
(Figures 7a and 6a). It is also interesting to observe that the tertiary circulations are absent for Pr = 0.7 and Da = 10
3 because of intense convective motion. As j gradually increases, the left secondary circulation cells grow bigger in size and form two asymmetric counterrotating fluid circulation cells with identical |ψ|max for Pr = 0.7 (Figure 7b,c). The multiple fluid circulation cells near the top portion of the cavity, as observed for Pr = 0.015, are absent for Pr = 0.7 at j = 75
(Figure 7c). Heatlines take a longer path from the hot bottom wall to the cold side walls in contrast to the earlier case with Pr = 0.015 for j = 30
. Figure 7a illustrates that, at slightly higher Pr (Pr = 0.7), a larger portion of the right half of the bottom wall delivers heat to the top portion of the cold walls based on dense heatlines. Also, heat transfer to the top right portions of the cavity is enhanced, as seen from a large thermal gradient based on highly dense heatlines. In contrast to Pr = 0.015, secondary heatline cells near the left portion of the cavity grow bigger, with |Π|max almost identical with that of the primary heatline cell for j = 45
at Pr = 0.7 due to convection induced by flow circulation, as shown in Figure 7b. As j further increases to 75
, secondary heatline cells further increase in size and two asymmetric heatline cells occupy the entire part of the cavity. The magnitudes of the heatfunctions are the same for both heatline cells (|Π|max = 8), despite the fact that they are asymmetric with respect to the center symmetric line. Also, the intensity of closed-loop heatline cells in the cavity increases with j. The boundary layer thickness near the top portion of the right cold wall is greatly reduced because of high heat flow, as

ARTICLE

represented by heatlines of high magnitude of heatfunctions in that regime compared to Pr = 0.015. Also, the large zone of the top central portion of the cavity remains at uniform temperature (θ = 0.3
0.5 for j = 30
and θ = 0.5
0.6 for j = 45
and 75
) because of enhanced heat transfer, as seen from the heatfunction contours of Figure 7a
c. Because of geometrical asymmetry, the heat is not evenly distributed to the cold side walls, and the left cold wall gets more heat compared to the right cold wall, as seen from the heatlines. Therefore, isotherms are found to be distorted more near the left portion of the cavity for j = 30
because of the high thermal gradient near those regions (see Figure 7a). As j increases to j = 75
, the dense heatlines rise through the central symmetric line and distribute almost an equal amount of heat to the cold side walls. Thus, a larger temperature distribution, θ = 0.5
0.6, is observed at the central regime of the cavity for j = 75
(see Figure 7c). On the basis of the highly dense heatlines, isotherms are compressed near the left and right bottom corners of the cavity for all j values. Thus, the boundary layer thickness near the left and right corners of the bottom wall is less than that near the middle portion of the bottom wall because of the high thermal gradient near the corners of the bottom wall for all j values (see Figure 7a
c). As Pr increases further (Pr = 7.2 and 1000), the intensity of the flow circulations increases because of high momentum diffusivity compared to the case of fluids with Pr = 0.015 and 0.7 (see Figure 8). Because of higher momentum diffusivity, the left circulation cells at j = 30
for Pr = 0.7 (Figure 7a) tend to increase in size and form longitudinally expanded fluid circulation cells within the left portion of the cavity for Pr = 7.2 (figures not shown) and Pr = 1000 (see Figure 8a). The qualitative trends of the streamlines and heatlines for higher Pr (Pr = 7.2 and 1000) are similar to that of Pr = 0.7 for all j values (Figures 7a
c and 8a
c). The stratification zone of temperature is formed near the middle portion of the bottom wall because of highly dense heatlines, as seen from heatline contours (Figure 8a
c). A large part of the top portions of the cavity is maintained at uniform temperature with θ = 0.4
0.6 for j = 30
because of the high heat-transfer rate near the top portions of the cavity, as represented by highly dense heatlines (Figure 8a) compared to Pr = 0.7 (Figure 7a). There is no further reduction in the boundary layer thickness near the top portion of the cold walls at Pr = 7.2 and 1000 compared to Pr = 0.7 for j = 45
and 75
because of the presence of almost identical magnitudes of heatlines in that zone for higher Pr, as seen in Figure 8b,c. Similar to Pr = 0.7, the top portion of the cavity is maintained at uniform temperature (θ) of about 0.5
0.6 for j = 45
and 75
at Pr = 7.2 and 1000 (see Figures 7b,c and 8b,c). It may also be noted that the boundary layer thickness near the bottom corners of the cavity is smaller compared to the middle portion of the bottom wall because of the presence of high dense heatlines near the corners of the bottom wall for all j values at higher Pr (Pr = 7.2 and 1000) similar to Pr = 0.7 due to high thermal gradient near those regions. 3.3. Case 2: Nonisothermally Heated Bottom Wall. Figures 9a
c
11a
c illustrate the streamlines, heatlines, and isotherms for nonisothermal heating of the bottom wall with Pr = 0.015 and 1000. Although streamlines exhibit trends qualitatively similar to those of the isothermal heating case for Pr = 0.015 and Da = 10
5, as seen in Figures 4 and 9, both the magnitudes of the streamfunctions and the intensity of the circulations are found to be higher for the isothermal heating case (case 1) with all j values 2122

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

Figure 9. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 2 with Pr = 0.015, Ra = 106, and Da = 10
5 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

Figure 10. Streamfunction (ψ), heatfunction (Π) and temperature (θ) contours for case 2 with Pr = 0.015, Ra = 106 and Da = 10
3 for (a) j = 30
, (b) j = 45
and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

because of larger heating effects. In contrast to isothermal heating, heatlines are not perpendicular to the bottom and left walls of the cavity, and they are not smooth circular arcs because of significant conductive heat flux (
∂θ/∂X = π cos πX) in the case of the nonisothermal sinusoidal heating pattern (see Figures 9a
c and 4a
c). In contrast to the previous case with isothermal heating, the magnitudes of the heatfunctions are small near the top portion of the cavity because of less heat transfer from the bottom wall to the top portion of the cavity. Therefore, the boundary layer thickness is comparatively higher near the top portion of the right wall for all j values. It may also be noted that the top portion of the cavity is maintained at θ e 0.1 for all j values because of the poor thermal mixing based on conduction-dominant heat transfer at Da = 10
5. Therefore, fluid near the top portion of the cavity is cooler than that with the isothermal heating case (case 1) for all j values as displayed via isotherms (see Figures 9a
c and 4a
c). Because of nonisothermal heating, singularity in temperature is absent at the bottom corners of the cavity and that results in a smooth temperature distribution for the nonisothermal heating case for

ARTICLE

Figure 11. Streamfunction (ψ), heatfunction (Π), and temperature (θ) contours for case 2 with Pr = 1000, Ra = 106, and Da = 10
3 for (a) j = 30
, (b) j = 45
, and (c) j = 75
. Clockwise and anticlockwise flows are shown via negative and positive signs of streamfunction and heatfunction, respectively. White and black shades in temperature (θ) contours represent hot and cold zones, respectively.

all j values (see Figure 9a
c) in contrast to the isothermal heating case (Figure 4a
c). This is also due to lesser heating effects, as indicated by less dense heatlines near the right bottom corner of the cavity in the case of nonisothermal heating on the bottom wall compared to the isothermal heating case for all j values. In contrast to case 1, the boundary layer thickness near the cold surfaces is larger for the nonisothermal heating case for all j values because of lesser heating effects (sinusoidal pattern) near those regions. Case 2 exhibits qualitatively similar streamlines, heatlines and isotherms trends of case 1 for Pr = 0.015 and Da = 10
3 [see Figure 10(a-c) and Figure 6(a-c)]. As expected, both magnitudes of streamfunctions and heatfunctions are found to be comparatively lesser for the present case than the isothermal heating case for all js [see Figure 10(a-c) and Figure 6(a-c)]. It is found that local heat circulation cells near the left bottom corner of the cavity for case 2 are not as strong as that of the isothermal heating case (case 1) because of the low intensity of heating in the left corner zone for case 2. Also, heatlines are not distorted much near the top right and bottom left portions of the cavity for the nonisothermal heating case compared to the isothermal heating case irrespective of j at Pr = 0.015 and Da = 10
3. This is due to the low intensity of heat circulation due to the nonisothermal boundary conditions near those zones, as seen in Figure 10a
c. Consequently, the zone near the top right corner is found to be cooler in contrast to the previous case with isothermal heating (see Figures 10a
c and 6a
c). The central core of the cavity is maintained at lower values of temperature θ = 0.1
0.3 at j = 30
, θ = 0.1
0.4 at j = 45
, and θ = 0.2
0.4 at j = 75
because of the lower thermal mixing effect in contrast to the isothermal heating case, where the core is maintained at θ = 0.2
0.4 for j = 30
, θ = 0.2
0.5 for j = 45
, and θ = 0.3
0.6 for j = 75
(see Figures 10c and 6c). Parts a
c of Figure 11 illustrate the streamlines, heatlines, and isotherms for Pr = 1000 and Da = 10
3 in the nonisothermal heating case, which are qualitatively similar to that of the isothermal heating case. Although the intensity of the flow and heat circulations are found to be stronger for Pr = 1000 than those for Pr = 0.015, but the intensity of the fluid and thermal mixing near the core is lower compared to that of the isothermal heating case for all j values at Pr = 1000. 2123

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research Heatlines are not distorted much near the top right portion of the cavity for j = 30
and 45
. This shows that the intensity of the heat circulation in that regime is lower in contrast to the isothermal heating case for same Pr and Da. Therefore, the boundary layer thickness near the top portion of the cavity for j = 30
and 45
is larger for the present case compared to the isothermal heating case (see Figures 11a,b and 8a,b). It is found that isotherms (θ = 0.1
0.4) are more disperse near the cold side walls and the central core of the cavity is maintained at θ = 0.4
0.5 for j = 75
at Pr = 1000 and case 2 in contrast to Pr = 1000 and case 1 (see Figures 11c and 8c). Larger thermal mixing is observed for case 1 than for case 2, and that is illustrated via isotherms θ = 0.5
0.7 for all j values in case 1, whereas θ = 0.2
0.5 for j = 30
and θ = 0.4
0.6 for j = 45
and 75
occur in case 2 at Pr = 1000 (see Figures 11a
c and 8a
c). However, the zone of thermal mixing is less for j = 30
, the zone is bigger for j = 45
, and the zone is maximum for j = 75
irrespective of the cases (cases 1 and 2) considered. Overall, the heat-transfer rate in the nonisothermal heating case is less than that of the isothermal heating case and that is clearly indicated by small magnitudes of streamfunctions, heatfunctions, and shades of isotherm contours in the gray-scale plot. 3.4. Heat-Transfer Rates: Local and Average Nusselt Numbers. The variation of the local heat-transfer rates for the bottom wall (Nub), left wall (Nul), and right wall (Nur) versus the distance in cases 1 and 2 for Pr = 0.015 and 1000 at Da = 10
3 is presented in the panel plots of Figure 12a
f. The heat-transfer rates are shown only at larger Da (Da = 10
3) as a representative situation of convection dominance for the brevity of the manuscript. It may be noted from eq 19 that the magnitude of the gradient of the heatfunction, (∂Π/∂Y or ∂Π/∂X), is equal to the temperature gradient (∂θ/∂X or ∂θ/∂Y) at the wall because of no-slip boundary conditions maintained at the wall. Note that the gradient of the temperature is equal to the Nusselt number based on eq 13. In this way, the heatfunction (Π) is inherently related to the Nusselt number, Nu, and, therefore, the heatfunction gradients may also be used for analysis of the variation of the Nusselt number on the walls of the cavity, apart from visualizing the heat flow. It may also be noted that dense heatlines represent high heat-transfer rates and vice versa and the magnitude of the heatfunctions represents the total heat flux being transferred along the path of heatline. Variation of the Nusselt numbers versus distance with the heatfunction gradient is explained first for Pr = 0.015 (see Figure 12a
c). It is observed that Nub is found to be infinite at the corners of the bottom wall for all cases of Pr and j, in the isothermal heating case (see the left panel of Figure 12a). This may be explained based on the gradients of the heatfunctions, which are also related to the local Nusselt number at the wall. The infinite value of Nub near corners of the bottom wall (see Figure 12a) is due to a large heatfunction gradient and that is clearly illustrated via gradients of Π corresponding to heatlines at the intersection of the hot bottom wall with the cold walls (see Figure 6a
c). It is noteworthy to mention that a sharp fall of Nub over the distance X e 0.3 for j = 30
, X e 0.19 for j = 45
, and X e 0.08 for j = 75
from the left bottom corner is (see the left panel of Figure 12a) due to a sharp decrease in the heatfunction gradient from the left bottom corner to the respective distances (see Figure 6a
c). Further, Nub tends to increase because of an increase in the heatfunction gradient, as represented by highly dense heatlines, and that is further due to the presence of

ARTICLE

secondary heatline cells. It is observed that Nub shows local maxima at about 0.5, 0.3, and 0.15 for j = 30
, 45
, and 75
, respectively, because this region is occupied by highly dense heatlines that are further due to secondary heatline cells, as illustrated in Figure 6a
c. Thereafter, Nub decreases and exhibits minima for j = 30
, 45
, and 75
at about X = 0.7, 0.5, and 0.3, respectively, because of a decrease in the heatfunction gradient, as represented by sparse heatlines with a small magnitude heatfunctions (Figure 6a
c). Further, Nub tends to increase because of an increase in the heatfunction gradient near the right portion of the bottom wall, as seen from dense heatlines near that region for j = 45
and 75
. The local Nub is the maximum at X = 0.82 for j = 75
in contrast to other j values, and that is due to a positive heatfunction gradient (heat flow changes from a clockwise to an anticlockwise direction) in contrast to j = 30
and 45
in that region (heat flow is in the same direction). It is also interesting to observe that the variation of Nub for the nonisothermal heating case (see the right panel of Figure 12a) does not show any maximum near the bottom corner regime in contrast to the isothermal heating case for all j values. In addition, the heat-transfer rate is smaller for nonisothermal heating with identical parameters (j), as discussed in section 3.3 (Figure 10a
c). Therefore, the maximum value of Nub in the leading corners of the bottom wall is smaller for the nonisothermal heating case (see the right panel of Figure 12a). Thereafter, Nub shows an increasing trend for all j values because of an increase in the heatfunction gradient, and that is further due to highly dense heatlines resulting in the presence of secondary heat circulation cells (Figure 10a
c), which further show maxima in Nub near the left hot
cold junction. This is clearly illustrated in the Nub distribution, where local maxima are observed near X = 0.2, 0.18, and 0.3 for j = 30
, 45
, and j = 75
, respectively (see the right panel of Figure 12a). It may be noted that Nub shows one more additional maximum at about X = 0.4 for j = 30
because of a large heatfunction gradient based on highly dense heatlines, and that is further due to a strong secondary heatline cell near that zone compared to other j values, as seen from Figure 10a. Thereafter, Nub decreases, and that shows its minima at about X = 0.58, 0.45, and 0.35 for j = 30
, 45
, and 75
, respectively. The decrease in Nub is due to a sharp fall of the heatfunction gradient near the respective zones (Figure 10a
c). It is interesting to note that Nub shows a decreasing trend toward the right corner of the bottom wall for case 2 irrespective of j in contrast to case 1 (see Figure 12a). This is further due to a decrease in the heat-transfer rate based on a small heatfunction gradient for the nonisothermal heating pattern toward the right bottom corner region in contrast to the isothermal heating pattern, where a high heat-transfer rate is observed near that regime. It is observed that the heat-transfer rate (Nul) is infinite at the bottom corner of the left wall irrespective of j at case 1. The infinite maximum in Nul near the bottom corner of the left wall (see the left panel of Figure 12b) is due to the maximum heat flow in that regime for all j values, as seen from heatlines with a large heatfunction gradient, and that is further due to inclination of the cold left wall with the hot bottom wall (Figure 6a
c). Thereafter, Nul sharply decreases until the distance is 0.2 for j = 30
, 0.25 for j = 45
, and 0.1 for j = 75
because of a sharp decrease in the heatfunction gradient from the lower left corner region to the respective distances in the left cold wall (see Figure 6a
c). Further, Nul tends to increase based on highly dense heatlines because of secondary local convective heat 2124

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Figure 12. Variation of the local Nusselt numbers with distance for Pr = 0.015 [(a) Nub; (b) Nul; (c) Nur] and Pr = 1000 [(d) Nub; (e) Nul; (f) Nur] at Ra = 106 and Da = 10
3 in the presence of isothermal heating (case 1) and nonisothermal heating (case 2) on the bottom wall for inclination angles j = 30
( 3 3 3 ), 45
(- - -), and 75
(—). Plots for cases 1 and 2 are shown in the left and right panels, respectively.

circulation cells near the bottom portion of the left wall (see Figure 6a
c). The trend of Nul shows one more maximum and minimum at S = 0.25 and 0.4, respectively, for j = 30
because of the highly dense and sparse heatlines, respectively, for j = 30
.

Finally, Nul decreases and, thereafter, increases near the upper portion of the left wall because of sparse and dense heatlines in the respective regions (Figure 6a
c). The variation of Nul for case 2 is qualitatively similar to that of case 1 except that Nul does 2125

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research not show any infinite maximum near the left hot
cold junction for case 2 (see the right panel of Figure 12b). This is due to the nonisothermal variation in the heating pattern, as discussed earlier (Figure 10a
c). It is interesting to note that the maxima of Nul are not sharper for case 2 compared to case 1, which may be attributed to the sparse heatlines, and they are caused by less intense heat recirculation cells in the lower portion of the cavity. It is also found that Nur shows its maximum at the bottom corner of the right wall for case 1 at Pr = 0.015 because of the maximum heat flow at that zone (see the left panel plot of Figure 12c). As discussed earlier, heat transport to the right cold wall increases with an increase in j, as seen in Figure 6a
c. On the basis of heatfunction gradients, larger values of Nur occur at the bottom corner for j = 75
(see the left panel of Figure 12c). A sharp fall of Nur is observed from the bottom corner because of a decrease in the heatfunction gradient, as seen based on the magnitudes of the heatlines in Figure 6a
c. Thereafter, Nur starts to increase because of the high heat-transfer rate based on dense heatlines, which is further due to a large heatfunction gradient initiated by a primary heat circulation cell. Sparse heatlines with low magnitudes (|Π| e 0.1) in the outer region of the primary circulation cell result in a lesser value of Nur toward the top corner of the right cold wall. It is also interesting to note that Nur tends to almost zero at the top right corner (S g 0.8) for j = 30
and 45
, which is clearly represented by small magnitudes of heatfunctions (0.001 e Π e 0.005) at that corner for j = 30
and 45
(Figure 6a,b). However, Nur does not fall below zero for j = 75
near the top portion of the right wall based on |Π| values, as represented by 0.005 e Π e 0.5 in that zone (Figure 6c). The variation of Nur is qualitatively similar for case 2 (see the right panel of Figure 12c), except near the bottom right corner, where the heat-transfer rate is much lower than case 1 because of the nonisothermal heating effect (see Figure 10a
c). The variation of Nub is qualitatively similar for higher Pr (Pr = 1000; see the left panel of Figure 12d). Note that the variation of Nub near corners of the bottom wall for case 1 at Pr = 1000 is almost identical with the variation of Nub for case 1 at Pr = 0.015 for all j values. This may be explained based on gradients of heatfunctions in Figure 8a
c where the heatfunction gradient decreases from the left bottom corner and that increases in the region near the right bottom corner. It is also found that Nub shows only one local minimum for case 1 at 0.65 e X e 0.8 for j = 30
, 0.6 e X e 0.75 for j = 45
, and 0.45 e X e 0.6 for j = 75
(see the left panel of Figure 12d). This is due to a sharp fall of the heatfunction gradient near the intersection of two asymmetric heat circulation cells. It may be noted that local minima of Nub for Pr = 1000 are sharper than that for Pr = 0.015 in case 1 (see the left panel plots of Figure 12a,d) because of the relatively high sparse heatlines in the case of Pr = 1000 (see Figures 6a
c and 8a
c). It is also interesting to note that the variation of Nub for case 2 at Pr = 1000 is qualitatively similar to that of Pr = 0.015, and a similar explanation follows (see the right panel of Figure 12d). It is interesting to observe that almost symmetric distributions of Nub occur for all j values in case 2 (see the right panel of Figure 12d). In contrast to the Pr = 0.015 case, Nul for Pr = 1000 at j = 45
and 75
(see Figure 12e) shows an increasing trend after its local minima for both cases (cases 1 and 2) and that is due to an increase in the heatfunction gradient. This is further due to the presence of asymmetric heat circulation cells resulting in highly dense heatlines (see Figures 8a
c and 11a
c). In contrast to the Pr = 0.015 case, the variation of Nur for Pr = 1000 does not fall to zero near the upper portion of

ARTICLE

the right wall for j = 45
and 75
in both cases (cases 1 and 2; see Figure 12f). This is due to significant heat flow in the case of higher Pr fluids (Pr = 1000) near the top portion of the right wall for case 1 compared to case 2, as seen in Figures 8a,b and 11a,b. Similar to the low Pr case (Pr = 0.015), Nur for j = 75
also has significant values near the top portion of the right wall in both cases, as represented by 0.1 e Π e 4 in case 1 and 0.1 e Π e 2.5 in case 2 near the top portion of the right wall, as seen in Figures 8c and 11c. The overall effects on the heat-transfer rates via the average Nusselt numbers for the bottom, left, and right walls (Nub, Nul, and Nur) versus the logarithmic Darcy number for all of the cases are displayed in the panel plots of Figure 13a
c,d
f. Parts a
c of Figure 13show the variations for Pr = 0.015, while parts d
f of Figure 13 show the variations for Pr = 1000. The panel plots of Figure 13a
f illustrate that the average Nusselt numbers increase (Nub, Nul, and Nur) with the Darcy number irrespective of j and Pr. It may also be noted that the average heat-transfer rates (Nub, Nul, and Nur) are high at Da = 10
3 compared to those of Da = 10
5 (see the panel plots of Figure 13a
f) because of the high heatfunction gradient at higher Da. Variation of the average Nusselt numbers versus Darcy numbers via the heatfunction gradient is explained first for Pr = 0.015, as shown in Figure 13a
c. Because of the higher inclination of the left cold wall toward the hot bottom wall (j = 30
), the heatfunction gradient is larger for the j = 30
cavity than for the j = 45
and 75
cavities, and thus maximum heat flow occurs from the bottom wall to the left cold wall for j = 30
. Overall, the maximum heat-transfer rate (maximum Nub and Nul) is observed for j = 30
, followed by 45
and 75
for case 1 at Pr = 0.015 (see the upper panel of Figure 13a,b). Because of thermal equilibrium between the walls of the cavity (Nub = Nul + Nur), the maximum heat flow from the hot bottom wall to the cold left wall results in a minimum heat flow to the right cold wall. Therefore, a minimum heat-transfer rate (minimum Nur) is observed for j = 30
, whereas a maximum heat-transfer rate (maximum Nur) is observed for 75
at Pr = 0.015 (see the upper panel of Figure 13c). It is also interesting to note that the right heatline cells at Da = 10
4 grow in size and tend to form strong primary circulation cells for all j values at Da g 10
4 (figures not shown). Hence, the right cold wall receives more heat, and therefore Nur shows a local peak at Da = 1  10
4 for j = 30
and 45
, and that is more pronounced at j = 75
at Da = 2  10
4 (see the upper panel of Figure 13c). The variations of Nub, Nul, and Nur for case 2 are qualitatively similar to those of case 1 except Nub at higher Da for Pr = 0.015. There is no significant difference in the heat-transfer rates for the bottom wall with j (Nub) at Da = 10
3 (see the lower panel of Figures 13a and 6a
c). The gradient of the heatfunctions demonstrates that the heat-transfer rate is smallest for j = 75
. Even though the variations of Nub, Nul, and Nur for the nonisothermal heating case are qualitatively similar to those of the isothermal heating case (see the lower panel of Figure 13a
c), the extent of variations with j and magnitudes of Nub, Nul, and Nur are found to be lesser compared with the isothermal heating case for particular values of Da. This is due to lower heating effects in the nonisothermal heating case compared to the isothermal heating case, as seen from the heatfunction contours of Figures 4
6 and Figures 9 and 10. The average Nusselt number plots are qualitatively similar for Pr = 1000 except Nub and Nur at j = 45
and 75
(see Figure 13d
f). The variation of Nub at j = 45
and 75
shows 2126

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Figure 13. Variation of the average Nusselt numbers with the Darcy number for Pr = 0.015 [(a) Nub; (b) Nul; (c) Nur] and Pr = 1000 [(d) Nub; (e) Nul; (f) Nur] in the presence of isothermal (case 1) and nonisothermal heating (case 2) on the bottom wall for inclination angles j = 30
( 3 3 3 ), 45
(- - -), and j = 75
(—). Plots for cases 1 and 2 are shown in the upper and lower panels, respectively.

a relatively profound increase with Da for Pr = 1000 compared to that in Pr = 0.015 at Da g 10
4. This is due to a comparatively large increase in the heatfunction gradient for j = 45
and 75
at Pr = 1000 compared to Pr = 0.015 at Da g 10
4 in case 1 (see the upper panel of Figure 13d). It may be noted that the maximum average heat-transfer rate is observed at the bottom wall for j = 45
and 75
at higher Da (Da g 10
4) for Pr = 1000 with case 2, which is in contrast to Pr = 0.015 (see the lower panel of Figure 13d). This may be explained from the spatial Nub variations of case 2 at Pr = 1000 and Da = 10
3 (see the right panel of Figure 12d), where j = 45
and 75
correspond to higher Nub values along the bottom wall except near the left and middle portions of the bottom wall. This is also illustrated from the heatfunction contours of Figure 11, where large heatfunction gradients are observed for j = 45
and 75
compared to j = 30
at Pr = 1000 and Da = 10
3. The variation of Nul for higher Pr (Pr = 1000) is qualitatively similar to that of lower Pr (Pr = 0.015) for all j values except that the variation of Nul with Da is less pronounced beyond Da g 10
4 for Pr = 1000 (see Figure 13e). This is due to the presence

of closed-loop heatline cells in the case of Pr = 1000, and that further obstructs the heat flow near the lower portion of the left wall for all j values, as discussed in the earlier section (see Figure 8a
c). A sharp increase in Nur with Da at a high Darcy regime for j = 45
and 75
is due to a large increase in the heatfunction gradient, which is further due to an increase in the heat-transfer rate from the bottom wall to the right wall at Pr = 1000 for case 1 (see the upper panel of Figure 13f). 3.5. Cup Mixing Temperature. All of the test cases (j values) are quantitatively compared in terms of the “cup mixing temperature (Θcup)” and “spatial average temperature (Θavg)” to demonstrate the thermal mixing effects. The bulk mean temperature in the cavity is estimated via the cup mixing temperature Θcup, which is defined by eq 26. The cup mixing temperature is adequately represented by the velocity-weighted temperature at convection-dominant heat transfer. Moreover, the spatial or area averaged temperature Θavg (eq 27) is also presented in the same panel plot of Θcup (Figure 14), for the sake of comparison of all j values corresponding to various fluids (Pr = 0.015
1000) at Da = 10
5 and 10
3 at Ra = 106 (see Figure 14a
d). In order to 2127

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

Figure 14. Cup mixing temperature (O) and average temperature (+) for different inclination angles (j = 30
, 45
, 75
, and 90
) of the rhombic cavity for various Pr fluids and Da at Ra = 106. Plots for cases 1 and 2 are shown in the upper and lower panels, respectively.

illustrate the role of j values, Θcup and Θavg are also calculated for a cavity with j = 90
. The detailed distributions of streamfunctions and isotherms for j = 90
are not shown here for the brevity of the manuscript. The cup mixing temperatures for all j values (j = 30
90
) are low and vary within 0.3
0.32 at low Da (Da = 10
5) because of the conduction-dominant heat flow for all Pr (see the upper panel plots of Figure 14a for Pr = 0.015) except Pr = 1000 in case 1 (see the upper panel plots of Figure 14b). On the other hand, j = 90
shows higher Θcup compared to all other j values in case 2 for Pr = 0.015 (see the lower panel plots of Figure 14a). Lower values of Θcup for rhombic cavities in case 2 may be due to the presence of cold regions, and that is further due to asymmetric heat distributions in the case of rhombic cavities in the nonisothermal heating pattern (see Figure 9a
c). At Pr = 1000, Θcup increases linearly with j and maximum cup mixing temperatures (Θcup = 0.44) are observed for j = 75
and 90
in case 1. A similar trend is also observed for case 2 (see Figure 14b). It may be noted that the higher values of Θcup indicate higher thermal mixing during convective flow at higher Da (Da = 10
3;

see Figure 14c,d). Although the heat-transfer rate is high at j = 30
, the thermal mixing in the cavity for Pr = 0.015 (molten metals) is highest for j = 75
(Θcup = 0.43) followed by 45
(Θcup = 0.38) and 30
(Θcup = 0.37) (see the upper panel plot of Figure 14c). This is due to the fact that high intense closed-loop heat circulation cells at j = 75
enhance the thermal mixing, as seen from the uniform temperature of θ = 0.3
0.5 at j = 75
, Da = 10
3, and Pr = 0.015 (see Figure 6c). It may also be noted that Θcup is found to be 0.29, 0.34, and 0.31 for j = 30
, 45
, and 75
, respectively, for case 2. This clearly illustrates that nonisothermal heating at the bottom wall causes poor thermal mixing inside the cavity compared to isothermal heating (see the upper and lower panel plots of Figure 14c). It may also be noted that j = 45
results in a maximum cup mixing temperature in case 2. This may also be explained based on the heatfunction contours of the corresponding case. It is found that, even though the strength of the primary heat circulation cells as well as the overall heat distribution is high for j = 75
, the secondary flow circulation cells near the left wall reduce the heat distribution as represented by the larger thermal boundary layer thickness at Da = 10
3 and 2128

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research Pr = 0.015 in case 2. Thus, Θcup at j = 75
is smaller than that for j = 45
, where heat distribution near the left wall is significantly higher at Da = 10
3 and Pr = 0.015 in case 2 (see Figure 10). Also, Θcup = 0.26 for case 1 and Θcup = 0.22 for case 2 at j = 90
indicates poor thermal mixing in the porous square cavity with boundary conditions identical with those of porous rhombic cavities (see Figure 14c). Overall, j = 75
shows enhanced thermal mixing and heat distribution in case 1, whereas j = 45
shows enhanced thermal mixing in case 2 for Pr = 0.015 (Figure 14c). As Pr increases, the momentum diffusivity increases, and that further results in high thermal mixing at Pr = 0.7. Analysis of thermal mixing for Pr = 0.7 indicates that the cup mixing temperature is found to be almost identical for all j values in case 1 [figures not shown]. However, Θcup is higher than that observed in Pr = 0.015 for all j values except at j = 75
. It is worthwhile to mention that asymmetric heatline cells suppress the heat flow near the lower portion of the cavity, and therefore larger variations of temperature occur near the lower portion of the cavity at j = 75
compared to 45
. On the other hand, uneven heat distribution near the left and right portions of the cavity at 30
results in a large variation in the temperature near the central core regime (see Figure 7a
c). Hence, j = 45
corresponds to enhanced thermal mixing at Pr = 0.7, and they are comparatively higher than Pr = 0.015. Case 2 also shows enhanced thermal mixing at j = 45
but lower than that in case 1. However, j = 90
shows higher Θcup than that of rhombic cavities because of symmetrical heat distribution throughout the cavity. Overall, the rhombic cavities of the j = 45
configuration with the case 1 heating strategy may be preferred for fluids with Pr = 0.7 (e.g., gases, air) used for electronic cooling and solar heating applications [figures not shown]. Convective flows in the case of higher Prandtl numbers such as Pr = 7.2 (water) and Pr = 1000 (olive/engine oil) result in almost identical thermal mixing for all j valuess, which correspond to Θcup = 0.42 for case 1 because of an increase in the momentum diffusivity (see the upper panels of Figure 14d). However, j = 45
still shows enhanced thermal mixing compared to other j values of the rhombic configuration in the nonisothermal heating case because of lesser heating effects in spite of the higher momentum diffusivity at Pr = 1000 (see the lower panels of Figure 14d). Finally, it may be inferred that porous rhombic cavities of all configurations with isothermal heating on the bottom wall may be used to enhance optimum thermal mixing for chemical solutions and oils. It is interesting to note that the square cavity shows a significant improvement in thermal mixing compared to the rhombic cavities in case 2 especially at higher Pr (Pr = 0.7
1000). 3.6. Conclusion. The concepts of heatlines and streamlines have been used to analyze the natural convection heat flow inside porous rhombic cavities to find the role of j for the strategic design of rhombic cavities with efficient thermal management. The heatlines concept has been implemented here for visualization of the heat flow inside the porous cavity to understand the direction of heat flow. The local Nusselt numbers (heat-transfer rate) are found to be correlated with heatfunction gradients. Cup mixing and average temperatures have also been evaluated to quantify the overall thermal mixing for various rhombic cavities. The results of low (Da = 10
5
10
4) and high (Da g 10
4) Darcy regions are highlighted as follows: 3.6.1. Low Darcy Region (Da = 10
5
10
4) • At lower Darcy number (Da = 10
5), heat transfer is conduction-dominant for all j values irrespective of Pr

ARTICLE

because of the low permeability of porous media even at higher Ra (Ra = 106). • Isotherms, streamlines, and heatlines are monotonic and smooth curves for all j values at Da = 10
5 irrespective of heating strategies. The left cold wall receives more heat from the large portion of the hot bottom wall for j = 30
. As j increases, heat is almost evenly distributed, as indicated by heatlines at Da = 10
5. • Distortion of the heatlines indicates that the onset of convection occurs at Da = 10
4 because of an increase in the permeability of porous media, and heatlines from the hot bottom wall take a longer path to reach the cold side walls of the cavity for all j values. Less intense heatline circulation cells, which are indicative of convective transport, also indicate the onset of convection heat transfer and thermal mixing at Da = 10
4 for all j values. • The maximum heat-transfer rate from the bottom wall (Nub) to the left wall (Nul) occurs for j = 30
followed by 45
and 75
for both heating cases irrespective of Pr at Da = 10
5. • The cup mixing temperature (Θcup) is invariant with j for all Pr except at Pr = 1000 in case 1, whereas j = 90
shows high Θcup compared to rhombic cavities for all Pr in case 2 at Da = 10
5. 3.6.2. High Darcy Region (Da g 10
4) • At higher Darcy number (Da g 10
4), an increase in the permeability of porous media further enhances the convection at identical Ra (Ra = 106). Multiple flow circulations are observed at Pr = 0.015 for all j values at Da = 10
3. Multiple circulation cells are suppressed, and two asymmetric flow circulation cells are found to occupy the entire cavity for all other Pr values (0.7, 7.2, and 1000) at j = 45
and 75
in case 1. Larger magnitudes of the streamfunction are found to occur near the core of the cavity for all j values at Pr = 0.015, whereas that occurs only for j = 30
at Pr = 0.7 and 7.2 for Da = 10
3 in the case of isothermal heating (case 1). • It is found that the right portion of the bottom wall contributes primarily for overall heat distribution in the right wall and partial heat distribution in the left wall for all j values at Pr = 0.015 and j = 30
at Pr = 0.7
1000, whereas heat is almost evenly distributed from the bottom wall for j = 45
and 75
at Pr = 0.7
1000. • The magnitudes of the streamfunctions and heatfunctions are significantly low for the nonisothermal heating case because of less heating effects, and the temperature distribution is smooth along the bottom wall of the cavity for all j values irrespective of Da and Pr. • The maximum heat-transfer rate from the bottom wall to the left wall (maximum Nub and Nul) occurs for j = 30
followed by 45
and 75
for case 1 irrespective of Pr. In contrast to case 1, the maximum heat-transfer rate from the bottom wall occurs for j = 75
at Da = 10
3 in case 2 at higher Pr (Pr = 1000). • Rhombic porous cavities of j = 75
with case 1 are found to be useful in liquid metal processing applications (Pr = 0.015) because of a larger cup mixing temperature over square cavities, whereas j = 30
with case 1 gives larger thermal processing rates for gases (Pr = 0.7), chemical solutions (Pr = 7.2), and oils (Pr = 1000).

’ APPENDIX The name “isoparametric” derives from the fact that the same parametric function describing the geometry may be used for the 2129

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research interpolating spatial variable within an element. Figure 2 shows a rhombus domain mapping to a square domain. The transformation between (x, y) and (ξ, η) coordinates can be defined by X = ∑9k=1Φk(ξ,η) xk and Y = ∑9k=1Φk(ξ,η) yk, where (xk, yk) are the (X, Y) coordinates of the k nodal points, as seen in Figure 2, and Φk(ξ,η) is the basis function. The nine basis functions are as follows: Φ1 Φ2 Φ3 Φ4 Φ5 Φ6 Φ7 Φ8 Φ9

¼ ¼ ¼ ¼ ¼ ¼ ¼ ¼ ¼

ð1
3ξ þ 2ξ2 Þð1
3η þ 2η2 Þ ð1
3ξ þ 2ξ2 Þð4η
4η2 Þ ð1
3ξ þ 2ξ2 Þð
η þ 2η2 Þ ð4ξ
4ξ2 Þð1
3η þ 2η2 Þ ð4ξ
4ξ2 Þð4η
4η2 Þ ð4ξ
4ξ2 Þð
η þ 2η2 Þ ð
ξ þ 2ξ2 Þð1
3η þ 2η2 Þ ð
ξ þ 2ξ2 Þð4η
4η2 Þ ð
ξ þ 2ξ2 Þð
η þ 2η2 Þ

The above basis functions are used for the mapping of the rhombus domain into a square domain and the evaluation of integrals of nonlinear residual equations of momentum, heat, and streamfunctions and heatfunctions.62,64

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ NOMENCLATURE Da Darcy number g acceleration due to gravity, m s
2 K permeability, m2 k thermal conductivity, W m
1 K
1 L side of the rhombic cavity, m N total number of nodes n normal direction on a plane Nu local Nusselt number Nu average Nusselt number p pressure, Pa P dimensionless pressure Pr Prandtl number R residual of the weak form Ra Rayleigh number modified Rayleigh number Ram S dimensionless distance along inclined walls T temperature of the fluid, K temperature of the hot wall, K Th temperature of the cold wall, K Tc u x component of velocity, m s
1 U x component of dimensionless velocity v y component of velocity, m s
1 V y component of dimensionless velocity V^ dimensionless velocity X dimensionless distance along the x coordinate x distance along the x coordinate, m Y dimensionless distance along the y coordinate y distance along the y coordinate, m Greek Symbols

R β

thermal diffusivity, m2 s
1 volume expansion coefficient, K
1

ARTICLE

γ Γ θ ν F j Φ ψ Π Ω ε

penalty parameter boundary of the two-dimensional domain dimensionless temperature kinematic viscosity, m2 s
1 density, kg m
3 inclination angle with the positive direction of the X axis basis functions dimensionless streamfunction dimensionless heatfunction two-dimensional domain error in the heat balance within the cavity

Subscripts

i k b l r s cup avg

global node number local node number bottom wall left wall right wall side wall cup mixing spatial average

’ REFERENCES (1) Kaviany, S. Principles of Heat Transfer in Porous Media; Springer-Verlag: New York, 1995. (2) Nield, D. A.; Bejan, A. Convection in Porous Media; Springer-Verlag: New York, 2006. (3) Ingham, D. B.; Pop, I. Transport Phenomena in Porous Media; Pergamon: Oxford, U.K., 2002; Vol. II. (4) Ingham, D. B.; Pop, I. Transport Phenomena in Porous Media; Elsevier: Oxford, U.K., 2005; Vol. III. (5) Ingham, D. B.; Bejan, A.; Mamut, E.; Pop, I. Emerging Technologies and Techniques in Porous Media; Kluwer: Dordrecht, The Netherlands, 2004. (6) Vafai, K. Handbook of Porous Media; Taylor & Francis: New York, 2005. (7) Pop, I.; Ingham, D. B. Convective Heat Transfer, Mathematical and Computational Modeling of Viscous Fluids and Porous Media; Pergamon: Oxford, U.K., 2001. (8) Bejan, A.; Dincer, I.; Lorente, S.; Miguel, A. F.; Reis, A. H. Porous and Complex Flow Structures in Modern Technologies; Springer: New York, 2004. (9) Wang, T. J.; Baek, S. W.; Kim, S. J. Re-evaluation of the Nusselt number for determining the interfacial heat and mass transfer coefficients in a flow-through monolithic catalytic converter. Chem. Eng. Sci. 2008, 63 (12), 3152–3159. (10) Cariaga, E.; Concha, F.; Sepulveda, M. Flow through porous media with applications to heap leaching of copper ores. Chem. Eng. J. 2005, 111, 151–165. (11) Saguy, I. S.; Marabi, A.; Wallach, R. New approach to model rehydration of dry food particulates utilizing principles of liquid transport in porous media. Trends Food Sci. Technol. 2005, 16, 495–506. (12) Hanspal, N. S.; Waghode, A. N.; Nassehi, V.; Wakeman, R. J. Numerical analysis of coupled Stokes/Darcy flows in industrial filtrations. Transp. Porous Media 2006, 64, 73–101. (13) Brandon, N. P.; Brett, D. J. Engineering porous materials for fuel cell applications. Philos. Trans. R. Soc. London, A 2006, 364, 147–159. (14) Ejlali, A.; Ejlali, A.; Hooman, K.; Gurgenci, H. Application of high porosity metal foams as air-cooled heat exchangers to high heat load removal systems. Int. Commun. Heat Mass Transfer 2009, 36 (7), 674–679. (15) Kamal, M. M.; Mohamad, A. A. Combustion in porous media: A review. J. Power Energy 2006, 220, 487–508. 2130

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research (16) Alvarez, G.; Flick, D. Modelling turbulent flow and heat transfer using macroporous media approach used to predict cooling kinetics of stack of food products. J. Food Eng. 2007, 80, 391–401. (17) Li, Z. W.; Dong, M. Z. Experimental study of carbon dioxide diffusion in oil-saturated porous media under reservoir conditions. Ind. Eng. Chem. Res. 2009, 48, 9307–9317. (18) Gatica, J. E.; Viljoen, H. J.; Hlavacek, V. Interaction between chemical reaction and natural convection in porous media. Chem. Eng. Sci. 1989, 44, 1853–1870. (19) Ma, H.; Guo, P.; Zhang, J. S.; Li, N.; Fu, B. M. M. Enhancement of oxygen transfer in liquid lead and lead
bismuth eutectic by natural convection. Int. J. Heat Mass Transfer 2005, 48, 2601–2612. (20) Ejlali, A.; Aminossadati, S. M.; Hooman, K.; Beamish, B. B. A new criterion to design reactive coal stockpiles. Int. Commun. Heat Mass Transfer 2009, 36 (7), 669–673. (21) Chen, C.; Zhang, D. X. Pore-scale simulation of density-driven convection in fractured porous media during geological CO2 sequestration. Water Resour. Res. 2010, 46, W11527. (22) Javaheri, M.; Abedi, J.; Hassanzadeh, H. Linear stability analysis of double-diffusive convection in porous media, with application to geological storage of CO2. Transp. Porous Media 2010, 84, 441–456. (23) Chen, Z. Q.; Gu, M. W.; Peng, D. H. Heat transfer performance analysis of a solar flat-plate collector with an integrated metal foam porous structure filled with paraffin. Appl. Ther. Eng. 2010, 30, 1967–1973. (24) Varol, Y.; Oztop, H. F.; Pop, I. Conjugate heat transfer in porous triangular enclosures with thick bottom wall. Int. J. Numer. Methods Heat Fluid Flow 2009, 19, 650–664. (25) Ahmed, N. J. S.; Badruddin, I. A.; Zainal, Z. A.; Khaleed, H. M. T.; Kanesan, J. Heat transfer in a conical cylinder with porous medium. Int. J. Heat Mass Transfer 2009, 52, 3070–3078. (26) Younsi, R. Computational analysis of MHD flow, heat and mass transfer in trapezoidal porous cavity. Ther. Sci. 2009, 13, 13–22. (27) Zhang, Z. Y.; Fu, C. J.; Tan, W. C.; Wang, C. Y. Onset of oscillatory convection in a porous cylinder saturated with a viscoelastic fluid. Phys. Fluids 2007, 19, 098104. (28) Baytas, A. C.; Pop, I. Free convection in oblique enclosures filled with a porous medium. Int. J. Heat Mass Transfer 1999, 42, 1047–1057. (29) Prasad, V.; Brown, K.; Tian, Q. Flow visualization and heat transfer experiments in fluid-superposed packed beds heated from below. Exp. Ther. Fluid Sci. 1991, 4, 12–24. (30) Bhattacharyya, S.; Dhinakaran, S.; Khalili, A. Fluid motion around and through a porous cylinder. Chem. Eng. Sci. 2006, 61 (13), 4451–4461. (31) Li, Z. W.; Dong, M. Z.; Shirif, E. Transient natural convection induced by gas diffusion in liquid-saturated vertical porous columns. Ind. Eng. Chem. Res. 2006, 45 (9), 3311–3319. (32) Das, S.; Sahoo, R. K.; Morsi, Y. S. Natural convection in heat generating oval porous enclosures: A non-Darcian model. Can. J. Chem. Eng. 2003, 81 (2), 289–296. (33) Bejan, A. Designed porous media: maximal heat transfer density at decreasing length scales. Int. J. Heat Mass Transfer 2004, 47, 3073– 3083. (34) Shuja, S. Z.; Yilbas, B. S.; Khan, S. M. A. Flow subjected to porous blocks in the cavity: Consideration of block aspect ratio and porosity. Chem. Eng. J. 2008, 139, 84–92. (35) Lorente, S.; Bejan, A. Heterogeneous porous media as multiscale structures for maximum flow access. J. Appl. Phys. 2006, 100, 114909. (36) Fedorov, A. G.; Viskanta, R. A numerical simulation of conjugate heat transfer in an electronic package formed by embedded discrete heat sources in contact with a porous heat sink. J. Electron. Packag. 1997, 119, 8–16. (37) de Lemos, M. J. S.; Fischer, C. Use of foam-like materials to enhance heat transfer from surfaces subjected to impinging jets. ASME Power Conf. 2009, 89–95. (38) Carrera-Rodriguez, M.; Martinez-Gonzalez, G. M.; NavarreteBolanos, J. L.; Botello-Alvarez, J. E.; Rico-Martinez, R. Numerical study

ARTICLE

of the effect of the environmental temperature in the natural twodimensional convection of heat in grain stored at cylindrical silos. Rev. Mex. Ing. Quim. 2009, 8, 77–91. (39) Das, D. B.; Lewis, M. Dynamics of fluid circulation in coupled free and heterogeneous porous domains. Chem. Eng. Sci. 2007, 62, 3549–3573. (40) Moukalled, F.; Darwish, M. Natural convection heat transfer in a porous rhombic annulus. Numer. Heat Transfer, Part A 2010, 58, 101–124. (41) Kimura, S.; Bejan, A. The heatline visualization of convective heat transfer. J. Heat Transfer Trans. ASME 1983, 105, 916–919. (42) Bello-Ochende, F. L. A heatfunction formulation for thermal convection in a square cavity. Int. Commun. Heat Mass Transfer 1988, 15, 193–202. (43) Morega, A. M.; Bejan, A. Heatline visualization of forcedconvection laminar boundary-layers. Int. J. Heat Mass Transfer 1993, 36, 3957–3966. (44) Zhao, F. Y.; Liu, D.; Tang, G. F. Application issues of the streamline, heatline and massline for conjugate heat and mass transfer. Int. J. Heat Mass Transfer 2007, 50, 320–334. (45) Costa, V. A. F. Unification of the streamline, heatline and massline methods for the visualization of two-dimensional transport phenomena. Int. J. Heat Mass Transfer 1999, 42, 27–33. (46) Costa, V. A. F. Unified streamline, heatline and massline methods for the visualization of two-dimensional heat and mass transfer in anisotropic media. Int. J. Heat Mass Transfer 2003, 46, 1309–1320. (47) Costa, V. A. F. Bejan’s heatlines and masslines for convection visualization and analysis. Appl. Mech. Rev. 2006, 59, 126–145. (48) Dash, S. K. Heatline visualization in turbulent flow. Int. J. Numer. Methods Heat Fluid Flow 1996, 6, 37–46. (49) Zhao, F. Y.; Liu, D.; Tang, G. F. Natural convection in an enclosure with localized heating and salting from below. Int. J. Heat Mass Transfer 2008, 51, 2889–2904. (50) Deng, Q. H.; Tang, G. F. Numerical visualization of mass and heat transport for conjugate natural convection/heat conduction by streamline and heatline. Int. J. Heat Mass Transfer 2002, 45, 2373–2385. (51) Morega, A. M.; Bejan, A. Heatline visualization of forced convection in porous media. Int. J. Heat Fluid Flow 1994, 15, 42–47. (52) Zhao, F. Y.; Liu, D.; Tang, G. F. Natural convection in a porous enclosure with a partial heating and salting element. Int. J. Therm. Sci. 2008, 47, 569–583. (53) Hooman, K.; Gurgenci, H. Heatline visualization of natural convection in a porous cavity occupied by a fluid with temperaturedependent viscosity. J. Heat Transfer Trans. ASME 2008, 130 (1), 012501. (54) Hooman, K.; Gurgenci, H.; Dincer, I. Heatline and energy-fluxvector visualization of natural convection in a porous cavity occupied by a fluid with temperature-dependent viscosity. J. Porous Media 2009, 12 (3), 265–275. (55) Waheed, M. A. Heatfunction formulation of thermal convection in rectangular enclosures filled with porous media. Numer. Heat Transfer, Part A 2009, 55, 185–204. (56) Saha, G.; Saha, S.; Islam, M. Q.; Akhanda, M. A. R. Natural convection in enclosure with discrete isothermal heating from below. J. Naval Arch. Marine Eng. 2007, 4, 1–13. (57) Yamanaka, Y.; Kakimoto, K.; Ozoe, H.; Churchill, S. W. Rayleigh
Benard oscillatory natural convection of liquid gallium heated from below. Chem. Eng. J. 1998, 71, 201–205. (58) Omri, A.; Orfi, J.; Ben Nasrallah, S. Natural convection effects in solar stills. Desalination 2005, 183, 173–178. (59) Varol, Y.; Oztop, H. F.; Pop, I. Numerical analysis of natural convection for a porous rectangular enclosure with sinusoidally varying temperature profile on the bottom wall. Int. Commun. Heat Mass Transfer 2008, 35, 56–64. (60) Saeid, N. F. Natural convection in porous cavity with sinusoidal bottom wall temperature variation. Int. Commun. Heat Mass Transfer 2005, 32, 454–463. 2131

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132

Industrial & Engineering Chemistry Research

ARTICLE

(61) Vafai, K.; Tien, C. L. Boundary and inertia effects on flow and heat transfer in porous media. Int. J. Heat Mass Transfer 1981, 24, 195–203. (62) Basak, T.; Roy, S. Role of “Bejan’s heatlines” in heat flow visualization and optimal thermal mixing for differentially heated square enclosures. Int. J. Heat Mass Transfer 2008, 51, 3486–3503. (63) Reddy, J. N. An Introduction to the Finite Element Method; McGraw-Hill: New York, 1993. (64) Basak, T.; Roy, S.; Matta, A.; Pop, I. Analysis of heatlines for natural convection within porous trapezoidal enclosures: Effect of uniform and non-uniform heating of bottom wall. Int. J. Heat Mass Transfer 2010, 53, 5947–961. (65) White, F. M. Viscous Fluid Flow; McGraw-Hill: New York, 1977. (66) Batchelor, G. K. Introduction to Fluid Dynamics; Cambridge University Press: Cambridge, U.K., 1993. (67) Lauriat, G.; Prasad, V. Non-Darcian effects on natural convection in a vertical porous enclosure. Int. J. Heat Mass Transfer 1989, 32, 2135–2148. (68) Walker, K. L.; Homsy, G. M. Convection in a porous cavity. J. Fluid Mech. 1978, 87, 449–474. (69) Bejan, A. On the boundary-layer regime in a vertical enclosure filled with a porous medium. Lett. Heat Mass Transfer 1979, 6, 93–102. (70) Gross, R. J.; Bear, M. R.; Hickox, C. E. The application of fluxcorrected transport (FCT) to high Rayleigh number natural convection in a porous medium. Proceedings of 8th International Heat Transfer Conference, San Francisco, CA. (71) Manole, D. M.; Lage, J. L. Numerical benchmark results for natural convection in a porous medium cavity. Heat and mass transfer in porous media; ASME Conference HTD-216, 1992; pp 55
60. (72) Goyeau, B.; Songbe, J. P.; Gobin, D. Numerical study of doublediffusive natural convection in a porous cavity using the Darcy
Brinkman formulation. Int. J. Heat Mass Transfer 1996, 39, 1363–1378. (73) Baytas, A. C.; Pop, I. Free convection in a square porous cavity using a thermal nonequilibrium model. Int. J. Therm. Sci. 2002, 41, 861–870. (74) Saeid, N. H.; Pop, I. Transient free convection in a square cavity filled with a porous medium. Int. J. Heat Mass Transfer 2004, 47, 1917–1924. (75) Saeid, N. H.; Pop, I. Natural convection from a discrete heater in a square cavity filled with a porous medium. J. Porous Media 2005, 8, 55–63.

2132

dx.doi.org/10.1021/ie2007856 |Ind. Eng. Chem. Res. 2012, 51, 2113–2132