Ind. Eng. Chem. Res. 2010, 49, 9771–9788
9771
Efficient Methodologies for Processing of Fluids by Thermal Convection within Porous Square Cavities Ram Satish Kaluri and Tanmay Basak* Department of Chemical Engineering, Indian Institute of Technology Madras, Chennai-600036, India
Effective use of thermal energy is the key for the energy-efficient processing of materials. A proper understanding of heat flow would be very useful in designing the systems of high energy efficiency with a minimal waste of precious energy resources. In the current study, a distributed heating methodology is proposed for the efficient thermal processing of materials. A detailed investigation on the processing of various fluids of industrial importance (with a Prandtl number of Pr ) 0.015, 0.7, 10, and 1000) in differentially and discretely heated porous square cavities is presented. Analysis of laminar convective heat flow within a range of Darcy number, Da ) 10-6-10-3 and Rayleigh number, Ra ) 103-106 has been carried out, based on a heatline visualization approach. The effect of Da and the role of distributed heating in enhancing the convection in the cavities is illustrated via heatline distributions, which represent the paths of the heat flow, the magnitude of heat flow, and zones of high heat transfer. It is observed that distributed heating plays an important role in enhancement of thermal mixing and temperature uniformity. Furthermore, the effect of Da for various Pr values on the variation of local Nusselt number (Nu) is analyzed, based on heatline distributions. 1. Introduction Porous materials continue to play an important role in various fields, such as chemical processing,1-7 food,8 and biomedical,9 environmental,4,10,11 reservoir,12 andelectrokineticapplications,13,14 as well as thermal management in electronics,15-18 fuel cells,19 electrochemical applications,20 combustion,21,22 building insulation,23 heat exchangers,24 etc. The subject of natural convection in porous media has significant importance, because many processes are governed by this phenomena. For example, Brooks et al.25 discussed the problem of spontaneous combustion of coal stockpiles and concluded that natural convection is the primary flow mechanism. Ejlali et al.26 investigated the key parameters such as approaching wind speed, porosity, and permeability of the porous medium and presented a criterion for the design of reactive coal stockpiles. Bejan27 presented a complete theory of the melting that occurs in a confined porous medium saturated with phase-change material and obtained the matched boundary layer solution for natural convectiondominated melting in the quasi-steady region. Gatica et al.28 investigated the interaction between chemical reaction and natural convection in porous media to determine the conditions under which secondary flows would develop, and showed that the presence of buoyancy-driven flows will drastically modify the topology of heterogeneous reaction system in an adiabatic packed bed. Benneker et al.29 experimentally illustrated the influence of natural convection on fluid mixing in packed beds at elevated pressure by varying flow rates, pressures, particle sizes, tube diameters, and flow directions, and they concluded that the usual correlations for mass and heat transfer may not be valid in high-pressure equipment. In light of numerous applications mentioned above, many studies have been conducted to gain a fundamental understanding of natural convection in porous media. A comprehensive experimental and numerical study on natural convection in a porous medium was carried out by Prasad et al.,30 and a model for effective thermal conductivity was proposed to account for the enhanced effect of fluid thermal conductivity during * To whom correspondence should be addressed. E-mail: tanmay@ iitm.ac.in.
convective flow. Chamkha31 numerically studied the transient hydromagnetic three-dimensional natural convection from an inclined porous stretching surface for metallurgical applications. Baytas and Pop32 numerically investigated the effects of natural convection on the heat transfer and temperature distribution within a tilted porous trapezoidal enclosure with cylindrical top and bottom walls at different temperatures and adiabatic sidewalls. Kathare et al.33 experimentally studied buoyant convection in superposed metal foam and water layers and concluded that the enhancement of heat transfer by natural convection could be obtained via varying the sublayer thickness of foam and its placement. In recent years, the approach of discrete/distributed heating is being considered as the potential methodology for enhanced thermal processing of materials. Distributed heating has advantages such as improved thermal mixing, greater temperature uniformity, greater control over flow field, and, importantly, high energy efficiency. A few interesting studies have been reported for enhanced processing of materials based on the distributed heating. Plumat34 reported the enhanced melting of glass by employing heated strips. Sarris et al.35,36 showed that the flow currents, temperature distribution, and thermal penetration could be controlled by heated strips for a homogeneous melting of glass which results in better quality of the final product. Enhancement of oxygen transport in liquid metal by natural convection during oxidation of liquid lead and lead-bismuth eutectic in a discretely heated enclosure was investigated by Ma et al.37 Studies on natural convection in porous media due to discrete heat sources are mainly focused on electronic cooling applications.38-41 However, thermal treatment of various fluids in porous media have several other important applications such as molten metal infiltration in porous media,42-45 drying and transport of gases in porous media,46-48 enhanced oil recovery by hot-water flooding in porous beds,49 combustion of heavy oils in porous reservoirs,50,51 etc. Thus, an understanding of the role of discrete heating for various fluids in porous media is important and such studies are yet to appear in the literature. The current work aims to address
10.1021/ie100569w 2010 American Chemical Society Published on Web 09/28/2010
9772
Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010
Figure 1. Schematic diagrams of the cavities filled with a fluid-saturated porous medium for different cases. The top wall is adiabatic. Thick lines represents the uniformly heated sections, while the remaining sections are maintained cold.
issues such as thermal mixing and temperature uniformity by distributed heating during natural convection in porous media, in the context of material processing applications. Numerical results of natural convection are mostly presented and analyzed by streamlines and isotherms. Although visual representation of fluid flow can be adequately illustrated by streamlines, isotherms are inadequate for convective heat flow visualization. The information of heat flow may be very useful, because that may be helpful with regard to devising strategies to control the temperature distribution in the cavity. A numerical tool to visualize heat flow was first proposed by Kimura and Bejan.52 They formulated a “heatfunction”, which is analogous to a streamfunction, by accounting for conductive and convective fluxes. The isolines of the heatfunction are called “heatlines”, and the trajectories of these heatlines represent the paths of heat flow. In addition, they indicate the magnitude of heat flow and zones of high heat transfer. The concept of heatlines has been employed and extended by several researchers to describe various physical phenomena.53-57 Few studies on the analysis of natural convection in porous media based on heatline concept have also been reported.58-60 In current work, the approach of heatlines is employed to study the heat flow distribution and thermal mixing within the square cavity in the presence of distributed heating of various walls. The main objective of the present study is to analyze the role of distributed heating and permeability of porous media (in terms of the Darcy number, Da) in enhancing the thermal mixing and temperature uniformity during natural convection in square cavities filled with fluid-saturated porous media, based on the heatline approach. Three different cases are considered: (1) hotisothermal bottom wall with cold-isothermal side walls, (2) discretely heated cavity with isothermal heat sources at the central regions of bottom and side walls, and (3) cavity with multiple isothermal heat sources, located at central as well as
lower corner regions. The porous medium is modeled based on a generalized non-Darcy model, neglecting the Forchheimer inertia term. This generalized model, based on volume averaging principles, was developed by Vafai and Tien.61 Cases 1-3 are studied for a range of parameters (Rayleigh number (Ra) of 103 e Ra e 106, Darcy number (Da) of 10-6 e Da e 10-3). Fluids of industrial importance, viz, molten metals (Pr ) 0.015), gases (Pr ) 0.7), aqueous chemical solutions (Pr ) 10), and olive/engine oils (Pr ) 1000) are used as working fluids. A detailed analysis on fluid and heat flow is presented, based on streamlines, heatlines, and isotherms. Thermal mixing and temperature uniformity in various cases is further analyzed using the cup-mixing temperature and the root-mean-square deviation (RMSD). Finally, the effect of Da on the heat transfer characteristics are studied using the local Nusselt number (Nu), and various qualitative and quantitative features of Nu are adequately explained, based on heatline distributions. 2. Mathematical Modeling and Simulation 2.1. Governing Equations: Temperature and Velocity. The physical domains for various cases are shown in Figures 1a-c. The top wall is maintained adiabatic in all of the cases. Case 1 involves an isothermal hot bottom wall with side walls being maintained isothermally cold (Figure 1a). Case 2 represents discrete heating of the cavity with heat sources applied along the walls of the cavity (Figure 1b), whereas Case 3 involves multiple heat sources, which are placed at the central portions of the walls and at the lower corners of the cavity. The location and the length of the heat sources may be noted from Figure 1c. The dimensionless length of the hot isothermal zone of the bottom wall is 0.5, while that along the side walls is 0.25. Note that the total dimensionless length of hot zones in Cases 2 and 3 is equal to that in case 1. All the physical properties are assumed to be constant,
Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010
except for the density in buoyancy term, and the change in density due to temperature variation is obtained using Boussinesq’s approximation. Also, it is assumed that the local thermal equilibrium (LTE) is valid, i.e., the temperature of the fluid phase is equal to the temperature of the solid phase within the porous medium. The momentum transfer in a porous medium is based on the generalized non-Darcy model of Vafai and Tien.61 However, the velocity squared term (or Forchheimer term), which models the inertia effect, is neglected in the present case, because this work involves only natural convection flow in a porous medium in an enclosed 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, the governing equations for steady twodimensional natural convection flow in a porous square cavity, using the conservation of mass, momentum, and energy, can be written with the following dimensionless variables or numbers: x X) , L P)
pL2 , FR2
y Y) , L Pr ) K)
ν , R
uL U) , R Da ) e3dp2
150(1 - e)
, 2
υL V) , R
T - Tc θ) Th - Tc
Ra )
gβ(Th - Tc)L3Pr 2
ν
∂U ∂V )0 + ∂X ∂Y
U
(1)
(2)
(
)
∂U ∂P ∂2U ∂U ∂2U Pr +V )+ Pr U + 2 ∂X ∂Y ∂X Da ∂X ∂Y2
(
continuity equation (eq 2) has been used as a constraint, because of mass conservation, and this constraint may be used to obtain the pressure distribution. 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 that is given by eq 2, which results in ∂V + ( ∂U ∂X ∂Y )
P ) -γ
(8)
The continuity equation (eq 2) is automatically satisfied for large values of γ. Typical value of γ that yield consistent solutions is 107. Using eq 8, the momentum balance equations (eqs 3 and 4) reduce to
U
(
)
∂U ∂ ∂U ∂V ∂2U ∂U ∂2U Pr +V )γ + U + Pr + 2 2 ∂X ∂Y ∂X ∂X ∂Y Da ∂X ∂Y (9)
(
)
and
K , L2
as
U
9773
(3)
)
∂V ∂P ∂2V ∂V ∂2V Pr +V )+ Pr V + RaPrθ + 2 2 ∂X ∂Y ∂Y Da ∂X ∂Y (4)
U
(
)
∂V ∂ ∂U ∂V ∂2V ∂V ∂2V +V )γ + + Pr + 2 2 ∂X ∂Y ∂Y ∂X ∂Y ∂X ∂Y Pr V + RaPrθ (10) Da
(
)
The system of equations (eqs 5, 9, and 10) with boundary conditions is solved using the Galerkin finite element method.62 Expanding the velocity components (U, V) and N temperature (θ), using basis set {Φk}k)1 , as N
N
∑ U Φ (X, Y),
U≈
k
∑ V Φ (X, Y),
V≈
k
k
k)1
and
k
k)1
N
∂θ ∂2θ ∂θ ∂2θ +V ) U + ∂X ∂Y ∂X2 ∂Y2
(5)
The boundary conditions for the velocities are U(X, 0) ) U(X, 1) ) U(0, Y) ) U(1, Y) ) 0 V(X, 0) ) V(X, 1) ) V(0, Y) ) V(1, Y) ) 0
k
k
(11)
k)1
for
(6) 0 e X,
and the boundary conditions for the temperature in Cases 1-3 are θ ) 1 (for the hot region) θ ) 0 (for the cold region) ∂θ ) 0 (for adiabatic topwall) ∂Y
∑ θ Φ (X, Y)
θ≈
The Galerkin finite element method yields the following nonlinear residual equations for eqs 9, 10, and 5, respectively, at nodes of the internal domain Ω:
(7)
Note that, in eqs 1-7, X and Y are dimensionless coordinates varying along the horizontal and vertical directions, respectively; U and V are dimensionless velocity components along the X- and Y-directions, respectively; θ is the dimensionless temperature; P is the dimensionless pressure; Ra is the Rayleigh number; Pr is the Prandl number; and Da is the Darcy number. The momentum and energy balance equations (eqs 3-5) are solved using the Galerkin finite element method. The
Ye1
N
R(1) i )
(
[( ) ]
N
∑ ∫ ∑U Φ Uk
k)1
N
Ω
k
k)1
)
k
∂Φk + ∂X
∂Φk VkΦk Φi dX dY + γ ∂Y k)1
∑
N
k)1
]
∂Φi ∂Φk dX dY + Ω ∂X ∂X
∑U ∫ k
k)1
∑ ∫ [ ∂X
N ∂Φi ∂Φk dX dY + Pr Uk Ω ∂X ∂Y k)1
∑V ∫ k
[
N
]
∂Φi ∂Φk Pr dX dY + ∂Y ∂Y Da
(
Ω
N
∫ ∑U Φ Ω
k
k)1
k
)
∂Φi ∂Φk + ∂X
Φi dX dY
(12)
9774
Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010 N
R(2) i )
(
[( ) ]
k)1
N
Ω
)
N
∑ ∫ ∑ Vk
UkΦk
k)1
[
∂Φk VkΦk Φi dX dY + γ ∂Y k)1
∑
N
Nus )
∂Φk + ∂X
]
N
∑U ∫
∂Φi ∂Φk dX dY + Ω ∂Y ∂X
k
k)1
∑ ∫ [ ∂X
N ∂Φi ∂Φk dX dY + Pr Vk Ω ∂Y ∂Y k)1
∑V ∫ k
k)1
]
∂Φi ∂Φk Pr dX dY + ∂Y ∂Y Da
(
N
∫ ∑V Φ Ω
RaPr
k
(
∂Φi ∂Φk + ∂X
Ω
)
k
k)1
N
Ω
k
k)1
∫
0.25
0
Nub,1 )
)
k
Nub dX
(20)
Nus dY
(for the left cold section)
0.25
(21)
Φi dX dY (13)
∫
0.75
0.25
Nub,2 )
Nub dX
(for the middle hot section)
0.5
and
(22) N
R(3) i )
1
0
In Cases 2 and 3, the cold and hot sections are present on the same wall, because of distributed heating. Therefore, the average Nusselt numbers for individual sections on the bottom wall in Case 2 are calculated as follows:
Φi dX dY -
∫ ∑θ Φ
∫
[(
N
∑θ ∫ ∑U Φ k
k)1
Ω
k
k)1 N
k
) ( ∂Φk + ∂X
∑θ ∫ [ k
k)1
Ω
N
∑V Φ k
k)1
k
) ]
∂Φk Φi dX dY + ∂Y
]
∂Φi ∂Φk ∂Φi ∂Φk + dX dY (14) ∂X ∂X ∂Y ∂Y
Biquadratic basis functions with three-point Gaussian quadrature are used to evaluate the integrals in the residual equations except the second term in eqs 12 and 13. In eqs 12 and 13, the second term, which contains the penalty parameter (γ), are evaluated with two-point Gaussian quadrature (reduced integration penalty formulation62). The nonlinear residual equations (eqs 12-14) are solved using the Newton-Raphson method to determine the coefficients of the expansions in eq 11. The detailed solution procedure is given in an earlier work.41 2.2. Nusselt Number and the Overall Heat Balance. The heat-transfer coefficient, in terms of the local Nusselt number (Nu) is defined by Nu ) -
∂θ ∂n
(15)
9
∂Φi
∑ θ ∂Y i
∫
1
0.75
Nub,3 )
Nub dX
(for the right cold section)
0.25
(23) In a similar manner, Nul, 1, Nul, 2, and Nul, 3 on the left vertical wall and Nur, 1, Nur, 2, and Nur, 3 on the right vertical wall are calculated. Note that subscripts “1” and “3” refer to cold sections, whereas the subscript “2” refers to the hot section on each wall. Following a similar manner, the average Nusselt number on hot and cold sections are evaluated in Case 3. The integrals are numerically evaluated using the adaptive trapezoidal rule. Furthermore, the convergence of integration results is checked with a heat balance, in terms of percentage error (discussed below). The heat balance was verified to ensure that the energy balance across the cavity is satisfied, i.e., the total heat lost by the hot sections is equal to the total heat gained by the cold sections. The heat balance equation for each case may be written as (i) Case 1: Nub ) Nul + Nur
where n denotes the normal direction on a plane. The local Nusselt numbers at the bottom wall (Nub), at the left wall (Nul), and at the right wall (Nur) are respectively defined as Nub )
and
(24)
(ii) Case 2: lb,2 ′ Nub,2 + ll,2 ′ Nul,2 + lr,2 ′ Nur,2 ) lb,1 ′ Nub,1 + lb,3 ′ Nub,3 + ll,1 ′ Nul,1 + ll,3 ′ Nul,3 + lr,1 ′ Nur,1 + lr,3 ′ Nur,3
(16)
i)1
(25)
(iii) Case 3: 9
∂Φi Nul ) θi ∂X i)1
∑
(17)
and
The error in heat balance is calculated as
9
∂Φi Nur ) θi ∂X i)1
∑
(18)
Nub )
∫
0
Nub dX
|∑ | |∑ | (| ∑ | | ∑ | ) Nh
Nc
-
lj,m ′ Nuj,m
Note that i is the local node number in an biquadratic element, which consists of nine local nodes and hence, the summation is written within the limits 1 to 9. In absence of any distributed heat sources, the average Nusselt number for the isothermal bottom and side walls (Nub and Nus, respectively) is given by 1
lb,1 ′ Nub,1 + lb,3 ′ Nub,3 + lb,5 ′ Nub,5 + ll,1 ′ Nul,1 + ll,3 ′ Nul,3 + lr,1 ′ Nur,1 + lr,3 ′ Nur,3 ) lb,2 ′ Nub,2 + lb,4 ′ Nub,4 + ll,2 ′ Nul,2 + ll,4 ′ Nul,4 + lr,2 ′ Nur,2 + lr,4 ′ Nur,4 (26)
(19)
)
lj,m ′ Nuj,m
hot
m)1
Nh
lj,m ′ Nuj,m
min
m)1
cold
m)1 Nc
lj,m ′ Nuj,m
,
hot
× 100
m)1
cold
(27)
where l′ refers to the length of cold or hot sections and the subscript j refers to the left (l), bottom (b), or right (r) walls. Subscripts m ) 1, 3 (cold sections) and 2 (hot section) on each
Ind. Eng. Chem. Res., Vol. 49, No. 20, 2010
wall in case 2 whereas in case 3, m ) 1, 3, 5 (hot sections) and 2, 4 (cold sections) and Nh and Nc refer to the total number of hot and cold sections in the cavity for given case. For each case, heat balance (eqs 24-26) is verified at all Ra values, and the error (), based on eq 27, is determined to be
