Analysis of Steady Convective Heating for Molten ... - ACS Publications

Material processing involving natural convection within a trapezoidal enclosure for uniformly and nonuniformly heated bottom wall, insulated top wall,...
1 downloads 0 Views 495KB Size
8652

Ind. Eng. Chem. Res. 2008, 47, 8652–8666

Analysis of Steady Convective Heating for Molten Materials Processing within Trapezoidal Enclosures Tanmay Basak,*,† S. Roy,‡,§ and E. Natarajan‡,§ Department of Chemical Engineering, Indian Institute of Technology Madras, Chennai 600036, India, and Department of Mathematics, Indian Institute of Technology Madras, Chennai 600036, India

Material processing involving natural convection within a trapezoidal enclosure for uniformly and nonuniformly heated bottom wall, insulated top wall, and isothermal sidewalls with inclination angle φ have been investigated. The penalty finite element method is used to obtain isotherm and streamline profiles for model liquids e.g. molten metal, salt water and olive oil. Parametric study for the wide range of Rayleigh number (Ra), 103 e Ra e 105 and Prandtl number (Pr) for model fluids with various tilt angles φ ) 45°, 30°, and 0° have been obtained. Secondary circulations were observed during molten metal processing. Streamlines show that the strength of convection is larger for φ ) 45° and flow intensities are also found to be larger for olive oil compared to molten metal and salt water. Heat transfer rates are shown via local and average Nusselt number plots. Local heat transfer rates are found to be relatively more for φ ) 0° than those with φ ) 45° and φ ) 30°. Average Nusselt number plots show higher heat transfer rates for φ ) 0° except for the nonuniform heating of the bottom wall with Pr ) 0.015 (molten metal). Overall, less heat transfer rates are observed for molten metal processing. 1. Introduction Convective heating processes are widely used in material processing industries (e.g, food, molten salt, molten metal etc).1-15 Studies on the various applications depend on the product specification, shape of the container, and heating characteristics. Numerical modeling can offer a way to reduce expensive experimental costs. A significant amount of earlier studies involve thermal processing of materials. Akterian1 developed a numerical model for the determination of the temperature field due to conductive heating of canned foods of various shapes under convective boundary conditions. Park et al.2 reported convective heat transfer in molten metals. Ulir´3 presented perspectives of molten salt reactors in nuclear industry. Basak and Ayappa4 and Ousegui et al.5 ensured that several attempts have been made to acquire a basic understanding of thawing processes and heat transfer characteristics in enclosures. Walle et al.6 reported pyrochemical techniques for nuclear systems based on molten salts. Moriyama et al.7 studied use of molten salts in fusion nuclear technology. Mishra and Olson8 reported various applications of molten salts in materials processing. Kumar et al.9 studied heating processes due to convection of different types of foods. Bergman and Hyun10 and Xu et al.11 reported convective heat transfer in molten metals. Pelekh et al.12 analyzed heating effects to determine the intrinsic kinetics of reactions involving nitrogen and molten titanium. Shang et al.13 analyzed the application of microwave heating on oil removal from the wastes. Various other studies involving sterilization and solidification of foods were reported by earlier researchers.14-26 Various applications on thermal processing further require a comprehensive understanding of heat transfer and flow circulations within cavities. Conduction, natural convection and forced convection are important means of heat transfer in thermal processing of * Corresponding author. E-mail address: [email protected]. † Department of Chemical Engineering. ‡ Department of Mathematics. § E-mail addresses: [email protected] (S.R.); [email protected] (E.N.).

food and applications of molten salts in nuclear industry. A comprehensive review by Bejan27 highlights that internal natural convection flow problems are more complex than external ones. During natural convection, the flow is solely due to the buoyancy force. The heat transfer rates are essential to obtain optimum processing conditions and to improve product quality. It is also essential to study the heat transfer characteristics in complex geometries in order to obtain the optimal design of the container. A few studies on steady natural convection within cavities have been carried out by earlier researchers, Gebhart28 and Hoogendoorn29 emphasized various aspects of natural convection flows in a square cavity. There were extensive studies using various numerical simulations reported by earlier researchers.30-37 Iyican and Bayazitoglu38 investigated the natural convective flow and heat transfer within a trapezoidal enclosure with parallel cylindrical top and bottom walls at different temperatures and plane adiabatic sidewalls. Karyakin39 reported two-dimensional laminar natural convection in isosceles trapezoidal cavity. The heat transfer rate is found to increase with the increase in base angle. Peric40 studied natural convection in a trapezoidal cavity using control volume method and observed the convergence of results for grid independent solutions. Kuyper and Hoogendoorn41 investigated laminar natural convection flow in trapezoidal enclosures to study the influence of the inclination angle on the flow and also the dependence of the average Nusselt number on the Rayleigh number. Thermosolutal heat transfer within trapezoidal cavity heated at the bottom and cooled at the inclined top part was investigated by Boussaid et al.42 Natural convection in cavities are also studied based on finite element method, and a few results are outlined by Reddy and Gartling.43 However, there is a conspicuous absence on studies of convective heating patterns for food preservation/ materials processing applications within trapezoidal containers which are commonly used for industrial processing. As a step toward the eventual development in the study of natural convection flows within enclosures for various potential practical applications, a comprehensive analysis on natural

10.1021/ie800263c CCC: $40.75  2008 American Chemical Society Published on Web 10/15/2008

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8653

convection flow within a trapezoidal enclosure with various thermal boundary conditions is important in materials processing research and that is yet to appear in the literature. The prime objective of this article is to investigate the circulation and temperature distribution with the detailed analysis of heat transfer rate in realistically shaped enclosures with uniformly and nonuniformly heated bottom plate and cooled side walls. The motivation of this study to analyze the situation where boundary conditions resemble a commonly used environment where the outside temperature is cold and inside temperature is sustained hot via any heating system. This study may also be useful to analyze the differential heating of materials with cold vertical walls. The analysis has been carried out for various materials with a range of Prandtl number (Pr) e.g. molten metals (Pr ) 0.026-0.004), gases (Pr ) 0.7-1), water (Pr ) 1.7-13.7), oils (Pr ) 50-105), etc., whereas earlier literatures primarily involve air and water only. The boundary conditions due to uniform heating correspond to jump discontinuities at corner points and similar boundary conditions were also used in earlier works on natural convection in square cavity.44,45 The boundary condition due to nonuniform heating has been represented by sinusoidal distribution of temperature and this type of boundary condition is particularly useful for processing molten glass.46 The effect of geometry has been illustrated for various angle of the sidewall varying within 0°-90°. The consistent penalty finite element method43-45,47 has been used to solve the nonlinear coupled partial differential equations for flow and temperature fields with both uniform and nonuniform temperature distributions prescribed at the bottom wall. Numerical results are presented in terms of isotherm and streamline contours along with the local and average heat transfer rates (Nusselt numbers). 2. Mathematical Formulation Let us consider a trapezoidal cavity of length of the bottom wall and height H with the left wall inclined at an angle φ ) 45°, 30°, and 0° with the y axis as seen in Figure 1a, b, and c, respectively. The velocity boundary conditions are considered as no-slip on solid boundaries. The liquid material is considered as incompressible, Newtonian, and the flow is assumed to be laminar. For the treatment of the buoyancy term in the momentum equation, the Boussinesq approximation is employed for the equation of the vertical component of velocity to account for the variations of density as a function of temperature and to couple in this way the temperature field to the flow field. The governing equations for steady natural convection flow using conservation of mass, momentum, and energy can be written as ∂u ∂V + )0 ∂x ∂y

(1)

(

∂u ∂u 1 ∂p ∂2u ∂2u u +V )+ν 2 + 2 ∂x ∂y F ∂x ∂x ∂y u

(

)

)

(2)

∂V 1 ∂p ∂2V ∂2V ∂V +V )+ ν 2 + 2 + gβ(T - Tc) (3) ∂x ∂y F ∂y ∂x ∂y

(

)

∂T ∂2T ∂2T ∂T +V )R 2 + 2 (4) ∂x ∂y ∂x ∂y The boundary conditions (see Figure 1) are as follows: At the bottom wall u

u ) 0, V ) 0, T ) Th,

πx or T ) (Th - Tc) sin + Tc L

( )

at the left and right vertical walls u ) 0, V ) 0, T ) Tc

(5)

and at the top wall ∂T )0 ∂y Here 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; T denotes the temperature; ν and R are kinematic viscosity and thermal diffusivity, respectively; p is the pressure, and F is the density; Th and Tc are the temperatures at hot bottom wall and cold vertical walls, respectively. Using following transformation of variables: u ) 0, V ) 0,

X)

T - Tc y uH VH pH2 x , Y) , U) , V) , P ) 2, θ ) H H R R Th - Tc FR (6)

The governing equations [eqs 1-4] reduce to nondimensional form: ∂U ∂V + )0 ∂X ∂Y U U

(7)

(

∂U ∂U ∂P ∂2U ∂2U +V )+ Pr + ∂X ∂Y ∂X ∂X2 ∂Y2

(

)

)

∂V ∂P ∂2V ∂2V ∂V +V )+ Pr + + RaPrθ ∂X ∂Y ∂Y ∂X2 ∂Y2

∂θ ∂2θ ∂2θ ∂θ +V ) + ∂X ∂Y ∂X2 ∂Y2 with the boundary conditions U

(8) (9) (10)

U ) 0, V ) 0, θ ) 1, or θ ) sin(πX), ∀Y ) 0, 0 e X e 1 U ) 0, V ) 0, θ ) 0, ∀X cos(φ) + Y sin(φ) ) 0, 0 e Y e 1 U ) 0, V ) 0, θ ) 0, ∀X cos(φ) - Y sin(φ) ) cos(φ), 0 e Y e 1 ∂θ ) 0, ∀Y ) 1, - tan(φ) e X e 1 + tan(φ) U ) 0, V ) 0, ∂Y (11)

Note that φ ) 45°, 30°, and 0° have been considered for simulation. Here X and Y are dimensionless coordinates varying along horizontal and vertical directions, respectively; U and V are dimensionless velocity components in the X and Y directions, respectively; θ is the dimensionless temperature; P is the dimensionless pressure; Ra and Pr are Rayleigh and Prandtl numbers, respectively. 3. Solution Procedure The momentum and energy balance equations [eqs 8-10] are solved using the Galerkin finite element method. The continuity equation [eq 7] is used as a constraint due to mass conservation andthisconstraintmaybeusedtoobtainthepressuredistribution.43-45 In order to solve eqs 8-9, we use the penalty finite element method where the pressure P is eliminated by a penalty parameter γ and the incompressibility criteria given by eq 7 (see the work of Reddy and Gartling43) which results in ∂V + ( ∂U ∂X ∂Y )

P ) -γ

(12)

The continuity equation [eq 7] is automatically satisfied for large values of γ. A typical value of γ that yields consistent solutions is 107.43-45

8654 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Using eq 12, the momentum balance equations [eqs 8 and 9] reduce to U

(

∂U ∂ ∂U ∂V ∂2U ∂2U ∂U +V )γ + + Pr + ∂X ∂Y ∂X ∂X ∂Y ∂X2 ∂Y2

(

)

)

(13)

and

(

)

∂V ∂ ∂U ∂V ∂2V ∂2V ∂V +V )γ + + Pr + + RaPrθ ∂X ∂Y ∂Y ∂X ∂Y ∂X2 ∂Y2 (14)

(

U

)

The system of equations [eqs 10, 13, and 14] with boundary conditions [eq 11] is solved using the Galerkin finite element method.43 Expanding the velocity components (U, V) and N temperature (θ) using basis set {Φk}k)1 as

∂ψ ∂ψ , V)∂Y ∂X which yield a single equation U)

∂2ψ ∂2ψ ∂U ∂V (20) + ) ∂X2 ∂Y2 ∂Y ∂X It may be noted that the positive sign of ψ denotes anticlockwise circulation, and the clockwise circulation is represented by the negative sign of ψ. The no-slip condition is valid at all boundaries as there is no cross-flow, hence ψ ) 0 is used for the boundaries. Expanding the stream function (ψ) using the basis set {Φ} as ψ ) ΣkN) 1 ψkΦk(X, Y) and the relation for U, V, the Galerkin finite element method yields the following linear residual equations for eq 20.

∑ ψ ∫ [ ∂X N

N

U≈

N

N

∑ U Φ (X, Y), V ≈ ∑ V Φ (X, Y), and θ ≈ ∑ θ Φ (X, Y) k

k

k

k)1

k

k

k)1

k

the Galerkin finite element method yields the following nonlinear residual equations for eqs 13, 14, and 10, respectively, at nodes of internal domain Ω:

[( ) ]

N

R(1) i )

(∑

)

N

∑U ∫ ∑U Φ k)1

k Ω

N

k

k

k)1

∂Φk + ∂X

[

∂Φk VkΦk Φi dX dY + γ ∂Y k)1 N

]

N

∑U ∫ k)1

k Ω

∂Φi ∂Φk dX dY + ∂X ∂X

∑ ∫[

N ∂Φi ∂Φk dX dY + Pr Uk Vk Ω ∂X ∂Y k)1 k)1

∑ ∫



∂Φi ∂Φk + ∂X ∂X

N

∑V ∫ ∑U Φ k)1

k Ω

N

VkΦk

k)1 N

k

k

k)1

N

Uk

k)1



N ∂Φi ∂Φk dX dY + Pr Vk Vk Ω ∂Y ∂Y k)1 k)1

∂Φi ∂Φk dX dY - RaPr ∂Y ∂Y

N

R(3) i )

[(

N

∑ ∫ ∑ θk

k)1



UkΦk

k)1 N

) ( ∂Φk + ∂X

∑ θ ∫ [ ∂X k)1

k Ω

k

k)1



]

∂Φi ∂Φk ∂Φi ∂Φk + dX dY ∂X ∂Y ∂Y

N

+

∑U ∫Φ k

k)1

i



N ∂Φk ∂Φk dX dY dX dY (21) Vk Φi ∂Y ∂X k)1 Ω

∑ ∫

The no-slip condition is valid at all boundaries as there is no cross-flow, hence ψ ) 0 is used as residual equations at the nodes for the boundaries. The biquadratic basis function is used to evaluate the integrals in eq 21 and ψ values are obtained by solving the N linear residual equations. 4.2. Nusselt Number. The heat transfer coefficient in terms of the local Nusselt number (Nu) is defined by ∂θ (22) ∂n where n denotes the normal direction on a plane. The normal derivative is evaluated by the biquadratics basis set in ξ-η domain. The local Nusselt numbers at bottom wall (Nub), left wall (Nul), and right wall (Nur) are defined as Nu ) -

∂Φk + ∂X

∂Φk Φi dX dY + γ ∂Y

and

]

∂Φi ∂Φk dXdY (16) ∂Y ∂Y

[( ) (∑ ) ] [∑ ∫ ∑ ∫ ] ∑ ∫[ ∫ (∑ ) ] N

R(2) i )

Rsi )

k)1

(15)

(19)

∂Φi ∂Φk dX dY + ∂Y ∂X Ω

∂Φi ∂Φk + ∂X ∂X

N



θkΦk Φi dX dY (17)

k)1

N

∑V Φ k

k)1

) ]

k

∂Φk Φi dXdY + ∂Y

]

∂Φi ∂Φk ∂Φi ∂Φk + dX dY (18) ∂X ∂Y ∂Y

The set of nonlinear algebraic equations (eqs 16-18) are solved using a reduced integration technique48 and the Newton-Raphson method as discussed in earlier work.44,45 4. Evaluation of Stream Function and Nusselt Number 4.1. Stream Function. The fluid motion is displayed using the stream function ψ obtained from velocity components U and V. The relationships between stream function, ψ (see the work of Batchelor49) and velocity components for twodimensional flows are

Figure 1. Schematic diagram of the physical system for (a) φ ) 45°, (b) φ ) 30°, (c) φ ) 0°.

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8655 9

Nub ) -

∑θ

i

i)1

(23)

∑ θ (cos φ

∂Φi ∂Φi + sin φ ∂X ∂Y

)

(24)

∑ θ (cos φ

∂Φi ∂Φi - sin φ ∂X ∂Y

)

(25)

9

Nul ) -

∂Φi ∂Y

i

i)1

and 9

Nur ) -

i

i)1

The average Nusselt numbers at the bottom, left, and right walls are



1

Nub )

0

Nub dX X|10

)



1

0

Nul ) cos φ



Nur ) cos φ



1/cosφ

0

Nub dX

(26)

Nul ds1

(27)

Nur ds2

(28)

and 1/cosφ

0

where ds1 and ds2 are the small elemental lengths along the left and right walls, respectively. Note that, Nul ) Nur ) Nus may be assumed due to symmetry in boundary conditions at side walls. 5. Results and Discussion 5.1. Numerical Tests. Equations 10, 13, and 14 with boundary conditions [eq 11] were solved numerically by penalty finite element method. The computational domain consists of 20 × 20 biquadratic elements which correspond to 41 × 41 grid points in the ξ-η domain as discussed in the Appendix and Figure 2. Numerical simulations were performed for Pr ) 0.015 (molten metals), Pr ) 7.2 (salt water with salinity of 35 g kg-1 at 20 °C) and Pr ) 988.24 (olive oil; see ref 50) with 103 e Ra e 105 for various inclination angles φ ) 45°, 30°, and 0° (see Figure 1a-c). The fluid motion and heating patterns were studied for heated bottom wall with cold isothermal side walls and adiabatic top wall. It may be noted that, the jump discontinuity in Dirichlet type of wall boundary conditions at the corner point (see Figure 1) corresponds to computational singularity. To ensure the convergence of the numerical solution to the exact solution, the grid sizes have been optimized and computed results are independent of grid sizes. In particular, the singularity at the corner nodes of the bottom wall needs special attention. The grid size dependent effect of the temperature discontinuity at the corner points upon the Nusselt numbers (local and average) tend to increase as the mesh spacing at the corner is reduced. One of the ways for handling the problem is assuming the average temperature of the two walls at the corner and keeping the adjacent grid-nodes at the respective wall temperatures. Alternatively, based on earlier work by Ganzarolli and Milanez,51 this procedure is still grid dependent unless a sufficiently refined mesh is implemented. Accordingly, once any corner formed by the intersection of two differently heated boundary walls is assumed at the average temperature of the adjacent walls, the optimal grid size obtained for each configuration corresponds to the mesh spacing over which further grid refinements lead to grid invariant results in both heat transfer rates and flow fields.

The stream functions and temperature contours have been compared with finite volume based solutions41 and the solutions are in well agreement. It may be noted that current solution is based on 20 × 20 biquadratic elements whereas the earlier work41 is based on 60 × 60 control volume grid. In the current investigation, Gaussian quadrature based finite element method provides the smooth solutions at the interior domain including the corner regions as evaluation of residual depends on interior Gauss points and thus the effect of corner nodes is less pronounced in the final solution. In general, the Nusselt numbers for finite difference/finite volume based methods are calculated at any surface using some interpolation functions40,41 which are now avoided in the current work. The present finite element approach offers special advantage on evaluation of local Nusselt number at the bottom and side walls as the element basis functions are used to evaluate the heat flux. Flow and temperature fields are shown in terms of streamlines and isotherms, respectively. On the basis of the symmetrical boundary conditions on the vertical walls, the flow and temperature fields are found symmetric about the midlength of the enclosure. The symmetrical boundary conditions in the vertical direction result in a pair of counter-rotating cells in the left and right halves of the enclosure for all the parametric values considered. 5.2. Uniform Heating at the Bottom Wall: Molten Metal vs Salt Water or Oil. Figures 3-8 illustrate the stream function and isotherm contours of numerical results for Ra ) 103-105 with Pr ) 0.015 (molten metal), Pr ) 7.2 (salt water), and Pr ) 988.24 (olive oil) in presence of uniformly heated bottom wall. As expected due to the cold vertical walls, fluids rise up from the middle portion of the bottom wall and flow down along the two vertical walls forming two symmetric rolls with clockwise and anticlockwise rotations inside the cavity. For Ra ) 103, the magnitudes of stream function are considerably smaller and the heat transfer is purely due to conduction as seen in Figure 3. For Ra ) 103 with φ ) 45°, the temperature contours with θ ) 0.1-0.4 occur symmetrically near the side walls of the enclosure (Figure 3a). The other temperature contours with θ g 0.5 are smooth curves symmetric with respect to vertical symmetrical line at the center. For φ ) 30°, the temperature contours with θ ) 0.1-0.3 occur symmetrically near the side walls of the enclosure (Figure 3b). The other temperature contours with θ g 0.4 are smooth curves symmetric with respect to vertical symmetric line. For φ ) 0° (square cavity), θ ) 0.1 is symmetric along the side walls and θ g 0.2 are smooth curves symmetric with respect to the vertical symmetric line (Figure 3c). It is interesting to note the stronger convection, however small, occurs for φ ) 45° and 30°. The central regime of top wall is at larger temperature for φ ) 45°. During conduction dominant heat transfer, the temperature profiles are almost invariant with respect to Ra and it is observed that the significant convection is initiated corresponding to a critical Ra, 3 × 103 for φ ) 45°, 3 × 103 for φ ) 30°, and 7 × 103 for φ ) 0° (Figure 4). At critical Ra, the distortion of the isotherm gradually increases and convection becomes dominant mode of heat transfer. At the onset of convection, the isotherms get distorted and move toward the sidewalls. As Ra increases, it is also observed that the circulations gradually become greater near the center. The isotherms with θ ) 0.5, 0.4, and 0.2 break into two symmetric contour lines for φ ) 45° (Figure 4a), 30° (Figure 4b), and 0° (Figure 4c), respectively. As Ra increases to 105, the effect of buoyancy force is stronger compared to viscous forces. Due to greater circulations near the central core at the top half of the enclosure, there are

8656 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

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

Figure 3. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 0.015 (molten metal) and Ra ) 103 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

small gradients in temperature at the central regime whereas a large stratification zone of temperature is observed at the vertical symmetry line near the bottom wall due to stagnation of flow [see Figure 5a-c]. It is also observed that the temperature contours are compressed strongly toward the side walls. For φ ) 45°, 30°, and 0°, the contour lines with θ e 0.4, θ e 0.3, and θ e 0.1, respectively, are compressed near the vertical walls of the enclosure. Isotherm lines for φ ) 45° are more distorted (Figure 5a) compared to 30° (Figure 5b) and 0° (Figure 5c). The shape of streamlines are found almost circular except near the side wall. It is interesting to observe that, for φ ) 45°,

symmetric secondary circulations with a maximum value of ψ ) 0.5 and tertiary circulations with ψ ) 0.2 near top corners were also observed (Figure 5a). For φ ) 30°, symmetric secondary circulations with a maximum value of ψ ) 0.1 were observed at all corners (Figure 5b) and for φ ) 0°, circulation patterns show neck formations and symmetric secondary circulation with the maximum value of ψ ) 0.2 and tertiary circulation with ψ ) 0.5 near the bottom corners are observed (Figure 5c). It may be remarked that these secondary circulations effects are the representatives for low prandtl number fluid Pr ) 0.015 (molten metal). The number

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8657

Figure 4. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 0.015 (molten metal) and Ra ) 3 × 103 for (a) φ ) 45°, Ra ) 3 × 103 for (b) φ ) 30°, and Ra ) 7 × 103 for (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

Figure 5. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 0.015 (molten metal) and Ra ) 105 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

of secondary/tertiary circulations is found to be dependent on tilt angle and therefore, the tilt angle may play significant role on thermal mixing within the fluid. For the salt water (Pr ) 7.2) at Ra ) 103, the isotherms are smooth and almost linear near the wall due to conduction dominant regime (figure not shown). The onset of convection takes place at critical Ra ) 3 × 103 for φ ) 45°, 3 × 103 for

φ ) 30°, and 5 × 103 for φ ) 0° as seen in Figure 6. Similar to molten metals, isotherms illustrate larger effect of convection for φ ) 45° (Figure 6a) than that for φ ) 30° (Figure 6b) and φ ) 0° (Figure 6c). For φ ) 45° and 30°, temperature contours with θ e 0.5 occur symmetrically near the sidewalls of the enclosure, whereas those for θ e 0.2 occur symmetrically near the sidewalls of the enclosure for φ ) 0°. The intensity of

8658 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Figure 6. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 7.2 (salt water) for (a) φ ) 45°, Ra ) 3 × 103, (b) φ ) 30°, Ra ) 3 × 103, and (c) φ ) 0°, Ra ) 5 × 103. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

Figure 7. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 7.2 (salt water) and Ra ) 105 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

circulation is stronger as seen from increased stream function (ψ) values for φ ) 45° and 30°. It is interesting to observe that the shape of streamlines near the wall start taking shape of the container (see Figure 6a-c). At Ra ) 105 (Figure 7), the isotherm contours get strongly compressed toward the sidewall as well as the bottom wall forming strong thermal boundary layers. The thickness of the boundary layer are reduced for salt water compared to molten metals (see Figures 5 and 7). At higher Ra for the salt water (Pr ) 7.2), the streamlines cover

the entire cavity and attain the shape of the container. For olive oil (Pr ) 988.24) at Ra ) 105 (Figure 8), temperature contours along the walls get further compressed and thickness of the thermal boundary layer is reduced. Similar to molten metals, streamlines indicate stronger circulations for φ ) 45° and 30°. The intensity of flow circulations for olive oil and salt water are represented with ψmax ) 17-18 for φ ) 45°, ψmax ) 17-18 for φ ) 30°, and ψmax ) 14-15 for φ ) 0° whereas for those for molten metals are ψmax ) 9 for φ ) 45°, ψmax ) 8 for φ )

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8659

Figure 8. Temperature and stream function contours for uniform bottom heating, θ(X,0) ) 1, with Pr ) 988.24 (olive oil) and Ra ) 105 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

Figure 9. Temperature and stream function contours for nonuniform bottom heating, θ(X,0) ) sin(πX), with Pr ) 0.015 (molten metal) and Ra ) 103 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

30°, and ψmax ) 5 for φ ) 0° (Figures 5, 7, and 8). These values further illustrate higher intensity of convection for olive oil and salt water compared to molten metals. It may also be remarked that the larger intensity of circulations for higher Pr fluid (salt water and olive oil) causes the shapes of stream functions almost trapezoidal near the walls and that signifies enhanced mixing effects.

5.3. Nonuniform Heating at the Bottom Wall: Molten Metal vs Salt Water or Oil. Figures 9-12 illustrate stream function and isotherm contours for Ra ) 103-105 during nonuniform heating of bottom wall. As seen in Figures 3-8, uniform heating of the bottom wall causes a finite discontinuity in Dirichlet type of boundary conditions for the temperature distribution at the edges of the bottom wall. In contrast, the

8660 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Figure 10. Temperature and stream function contours for nonuniform bottom heating, θ(X,0) ) sin(πX), with Pr ) 0.015 (molten metal) for (a) φ ) 45°, Ra ) 2 × 103, (b) φ ) 30°, Ra ) 2 × 103, and (c) φ ) 0°, Ra ) 8 × 103. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

nonuniform heating removes the singularities at the edges of the bottom wall and provides a smooth temperature distribution in the entire cavity. During conduction dominant heating (Ra ) 103) for molten metal (Pr ) 0.015) with φ ) 45°, the temperature contours with θ ) 0.1-0.2 occur symmetrically near the side walls of the enclosure as seen in Figure 9a, whereas those for θ e 0.4 occur symmetrically near the side walls for uniform heating case as seen in Figure 3a. Similar to φ ) 45°, temperature contours with θ ) 0.1-0.2 occur symmetrically near the side walls of the enclosure for φ ) 30° (Figure 9b). In contrast, for φ ) 0°, temperature contours with θ g 0.1 are smooth curves symmetric with respect to the vertical symmetric line (Figure 9c). It is observed that the significant convection is initiated at critical Ra ) 2 × 103 for φ ) 45°, Ra ) 2 × 103 for φ ) 30°, and Ra ) 8 × 103 for φ ) 0° (Figure 10a-c). At the onset of convection, θ e 0.3 occurs symmetrically for φ ) 45°, θ e 0.2 for φ ) 30°, and θ ) 0.1 for φ ) 0° near the side walls. Similar to uniform heating, streamlines show stronger convection for φ ) 45° than that with φ ) 30° and 0°. However, the strength of convection is less than that with uniform heating as seen in Figure 4a-c. At Ra ) 105, the temperature contours are compressed toward the side walls. For φ ) 0°, the contour line with θ ) 0.1 occur symmetrically near the side walls (Figure 11c). For φ ) 30°, the contour line θ ) 0.2 breaks into symmetric contour lines (Figure 11b). Further for φ ) 45°, the contour line θ ) 0.3 breaks into symmetric contour lines rising toward the middle of the enclosure due to stronger circulation at the center for nonuniform heating (Figure 11a). As seen from Figure 11a, for φ ) 45°, the isotherm lines are highly compressed near the bottom wall whereas the isotherms are less compressed near bottom wall for φ ) 30° and 0° due to less intense convective effects as seen in Figure 11b and c. Similar to uniform heating case, the stream lines are circular for Pr ) 0.015. It may also be noted that, for φ ) 45°, symmetric secondary circulation with ψ ) 0.2 and 0.5 are found near the side walls and weaker

tertiary circulations are also found for nonuniform heating case (see Figure 11a). For φ ) 30°, symmetric secondary circulation with ψ ) 0.2 and tertiary circulation with ψ ) 0.05 are found near the top portion of side walls whereas for φ ) 0°, secondary circulations are absent (Figure 11b and c). Similar to the uniform heating case, the stronger effect on convection occurs for salt water and olive oil. At Ra ) 105 (Figure 12), strong thermal boundary layers are formed along the sidewalls. The boundary layer thickness is much reduced for salt water compared to molten metals. However temperature gradient within the boundary layer is less compared with uniform heating see (Figures 7, 8, and 12). It may be noted that θ ) 0.1-0.4 are confined within the boundary layer for uniform heating case whereas θ ) 0.1-0.3 are confined within the boundary layer for nonuniform heating case. It is also interesting to observe that the strength of convection is lesser in nonuniform heating case (see Figures 7 and 12). It may be noted that ψmax ) 15 for φ ) 45° (Figure 12a) and 30° (Figure 12b) whereas ψmax ) 13 for φ ) 0° (Figure 12c). Similar to the uniform heating case, the circulations almost take the shape of the container similar to uniform heating for Pr ) 7.2. For olive oil (Pr ) 988.24), the thickness of thermal boundary layer is more compressed than that for salt water and molten metals and isotherm contours are highly dense near the side walls similar to uniform heating (figure not shown). The streamlines are almost trapezoidal near the walls for high Pr values (7.2 and 988.24) similar to uniform heating. 5.4. Heat Transfer Rates: Local and Average Nusselt Numbers. Figures 13 and 14 display the effects of Ra for model fluids (molten metals and salt water) on the local heat transfer rates or Nusselt numbers at the bottom and side walls (Nub, Nus) for various tilt angles φ. Figure 13a illustrates local Nusselt number distribution at the bottom wall (Nub) for Pr ) 0.015. As a result of symmetry in the temperature field, heat transfer at the bottom wall is symmetric with respect to the midlength

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8661

Figure 11. Temperature and stream function contours for nonuniform bottom heating, θ(X,0) ) sin(πX), with Pr ) 0.015 (molten metal) and Ra ) 105 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

Figure 12. Temperature and stream function contours for nonuniform bottom heating, θ(X,0) ) sin(πX), with Pr ) 7.2 (salt water) and Ra ) 105 for (a) φ ) 45°, (b) φ ) 30°, and (c) φ ) 0°. Clockwise and anticlockwise flows are shown via negative and positive signs of stream functions, respectively.

(X ) 1/2). Common to all cases (φ ) 0-45°) with uniform heating, the temperature contours are widely dispersed at the center of the bottom wall and therefore local Nusselt number has a minimum at X ) 1/2 for all Ra values. The local Nusselt number distribution has several interesting features as shown in Figure 13a. At Ra ) 103, the temperature contours are found to be most dense at the bottom wall for φ ) 0° whereas the temperature contours are widely dispersed at bottom for φ ) 45° (see Figure 3a-c). Therefore, local Nusselt

number for the bottom wall (Nub) corresponds to larger magnitude for φ ) 0° whereas it is smaller for φ ) 45°. The variation of local Nusselt number for φ ) 0-45° has been studied in detail for Ra ) 105. It is interesting to observe that the intensity of flow is significantly stronger for φ ) 45° as seen from the ψmax of primary vortex (see Figure 5a-c). The stronger circulation further enhances mixing within the core and that compresses the isotherms at the bottom. Therefore, overall Nub is larger for φ ) 45° except at bottom corner points.

8662 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008

Figure 13. Variation of local Nusselt number with distance for Pr ) 0.015 (molten metal) for at the (a) bottom wall and (b) side wall for uniform heating (s) and nonuniform heating (- - -). Inset plots show heat transfer rates for selected regimes with uniform heating (s) and nonuniform heating (- - -) for Ra ) 105.

The left inset plot illustrates Nub distribution within X ) 0.1-0.3 for φ ) 0-45° with Ra ) 105 for uniform heating. Due to stronger primary vortices, the temperature contours are compressed near X ) 0 for φ ) 0° and 45° leading to larger Nub values whereas due to secondary vortex formation at the corner, the temperature contours are less dense near X ) 0 for φ ) 30° leading to smaller magnitude of Nub. The local Nusselt number at X g 0.3 is larger for φ ) 45° due to dense temperature contours whereas Nub is smaller for φ ) 0° and 30° due to less temperature gradient resulting from less circulation. For nonuniform heating, due to the sinusoidal temperature field, the heat transfer at the bottom wall is also sinusoidal. At Ra ) 103, it is observed that Nub has a maxima at X ) 1/2 for φ ) 0° due to dense temperature contours. It is interesting to observe that Nub at X ) 1/2 for φ ) 0° is even larger than that for φ ) 45°. In contrast, the temperature contours are more dense for φ ) 45° than those for φ ) 30° and φ ) 0° with Ra ) 105. Therefore, the local Nusselt number is less at X ) 1/2 for φ ) 0°. The right inset plot for Ra ) 105 illustrates the Nub distribution within X ) 0.1-0.3 for φ ) 0-45°. Due to denser temperature contours toward the middle portion of bottom wall for φ ) 45°, the local Nusselt number for most of the regimes is larger for φ ) 45° than that for φ ) 30° and φ ) 0°. It is interesting to observe from inset plots that, Nub decreases with distance for φ ) 0-45° for uniform heating whereas Nub increases with distance for nonuniform heating. Due to singular temperature boundary conditions at bottom corners, the heat transfer rates are larger at the corner points for uniform heating. Heat transfer rates at the side wall for uniform and nonuniform heating cases are illustrated in Figure 13b. Due to

Figure 14. Variation of local Nusselt number with distance for Pr ) 7.2 (salt water) for at the (a) bottom wall and (b) side wall for uniform heating (s) and nonuniform heating (- - -). Inset plots show heat transfer rates for selected regimes with uniform heating (s) and nonuniform heating (- - -) for Ra ) 105.

symmetry in the boundary condition, the local Nusselt number is identical along both the side walls. Figure 13b illustrates that the local Nusselt number (Nus) decreases from the bottom edge with vertical distance at the side wall for Ra ) 103 for the uniform and nonuniform heating cases. The local Nusselt numbers at the side wall for both the heating case at Ra ) 103 are found to be identical for φ ) 45°, 30°, 0° and therefore Nus is shown for φ ) 45°. For Ra ) 103 (see Figures 3 and 9), it is observed that the temperature contours are widely spread toward the sidewalls for φ ) 45°, 30°, and 0°. It is interesting to note that, the heat transfer rate initially decreases and later it increases with distance for Ra ) 105. This is due to the fact that the temperature contours for Ra ) 105 are widely spread toward the top portion of the side walls forming thermal boundary layers. On the other hand, the temperature contours are compressed at the middle portion of the side wall due to larger intensity of circulation. For nonuniform heating case as seen in Figure 13b the heat transfer rates are less compared to uniform heating case. This may be due to the temperature contours to be less compressed toward the side walls compared to the uniform heating case (see Figures 5 and 11). It is interesting to note from the left inset plot for Ra ) 105, the local Nusselt numbers for φ ) 45° has increasing trend for the distance within 0.15-0.35, whereas for φ ) 0°, it shows decreasing trend and for φ ) 30°, both decreasing and increasing trends of Nus have been observed during uniform heating. In the region near to the bottom wall, the local Nusselt number is found to be higher for φ ) 0°. The right inset plot illustrates Nus distributions for nonuniform heating. The qualitative trends of Nus are similar to that of uniform heating. It may be noted

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8663

that, at higher Rayleigh numbers, the significant circulations as seen in Figures 5 and 11 result in highly dense contours at the top portion of the side walls and these dense temperature contours are in contrast with the conduction dominant cases as seen in Figures 3 and 9. Figure 14a illustrates the local heat transfer rates at the bottom wall for the salt water. For Ra ) 103, Nub exhibits similar qualitative distribution for molten metal as seen in Figure 13a. In addition, Nub is identical for φ ) 0°, 30°, and 45° for Ra ) 103. At Ra ) 105, the local Nusselt number for φ ) 0° is larger than that for φ ) 45° and φ ) 30° for both uniform and nonuniform heating cases. This is due to high temperature gradient for φ ) 0° as seen from the Figures 7 and 12. On the other hand, the local Nusselt number is found to be identical at a distance X ) 0.5 for φ ) 45° and 30°. The left inset plot within a distance 0.1-0.3 along the bottom wall shows higher heat transfer rate for φ ) 0° than that for φ ) 45° and φ ) 30° for uniform heating at Ra ) 105. At Ra ) 105, the nonuniform heating provides a sinusoidal type of local heat transfer rate symmetric with respect to midlength X ) 1/2. Similar to molten metals, the nonuniform heating shows larger Nub at X ) 1/2 for Ra ) 103. It may be noted that, the local Nusselt numbers for φ ) 45° and φ ) 30° are found to be identical whereas Nub is slightly larger for φ ) 0° at Ra ) 105. The right inset plot illustrates Nub distributions for nonuniform heating. The increasing trend of Nub is observed irrespective of φ and Nub at φ ) 0° is slightly larger within 0.1 e X e 0.3. It is also interesting to observe that, the local Nusselt number (Nub) shows the increase in the heat transfer rates for salt water compared to molten metals due to the increased convection effects of the salt water (see Figures 13a and 14a). Heat transfer rates for uniform and nonuniform heating at the side walls are shown in Figure 14b. For Ra ) 103, the heat transfer rates for uniform and nonuniform heating are almost identical for φ ) 45° except at corner points. At Ra ) 105, the heat transfer rates for φ ) 0° are larger except at bottom corner points as can be observed from the inset plot. It is interesting to note that Nus increases with distance except near the bottom corner point for all φ due to compression of temperature contours for larger intensity of circulation. In contrast, the compression is absent for Pr ) 0.015 and φ ) 0° and therefore Nus was observed to decrease. It may also be noted that, the local heat transfer rates are lesser at bottom points for nonuniform heating as the isotherms are dispersed. The qualitative trend of Nus distribution for nonuniform heating is similar to that with uniform heating for rest of the regimes. Similar to Pr ) 0.015, the nonuniform heating corresponds to lesser Nus than that with uniform heating. The overall effects upon the heat transfer rates are displayed in Figures 15 and 16, where the distribution of the average Nusselt number of the bottom and side walls, have been plotted vs the logarithmic Rayleigh number for Pr ) 0.015 of molten metal and Pr ) 7.2 of salt water, respectively. The average Nusselt numbers were obtained using eqs 26-28 where the integral is evaluated using Simpsons 1/3 rule. Note that Figure 15a and b illustrates heat transfer rates for uniform heating cases and Figure 15c and d represents nonuniform heating cases for molten metal. The values of the average Nusselt numbers along the side walls are less compared to the bottom wall irrespective of the tilt angle φ. This is due to the fact that the length of the bottom hot wall is lesser than the length of side or cold walls and also based on overall heat balance; Nub × lb ) 2Nusls, where lb ) length of the bottom wall and ls ) length of the side walls.

Figure 15. Variation of average Nusselt number with Rayleigh number for Pr ) 0.015 (molten metal) for uniform heating (a and b) and nonuniform heating (c and d).

Figure 16. Variation of average Nusselt number with Rayleigh number for Pr ) 7.2 (salt water) for uniform heating (a and b) and nonuniform heating (c and d).

For the uniform heating case, average Nusselt numbers decrease with increase in φ for specific Ra (see Figure 15a and b). This may be explained with the inset plot of Figure 13a and b

8664 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 Table 1. Correlations for the Average Nusselt Number of Molten Metals (Pr ) 0.015) Nub

Nus

φ

uniform heating

nonuniform heating

uniform heating

nonuniform heating

45° 30° 0°

1.7198Ra0.0896 2.6520Ra0.0528 4.3852Ra0.0231

0.4937Ra0.1439 0.7902Ra0.0934 1.0897Ra0.0705

0.6428Ra0.0915 1.1748Ra0.0541 2.2678Ra0.0192

0.1925Ra0.1413 0.3482Ra0.0952 0.5714Ra0.0652

Table 2. Correlations for Average Nusselt Number of Salt Water (Pr ) 7.2) Nub

Nus

φ

uniform heating

nonuniform heating

uniform heating

nonuniform heating

45° 30° 0°

0.6427Ra0.2204 0.8650Ra0.1980 1.2361Ra0.1765

0.1467Ra0.3168 0.1684Ra0.3057 0.1590Ra0.3148

0.2546Ra0.2145 0.3960Ra0.1948 0.6312Ra0.1737

0.0614Ra0.3033 0.0786Ra0.2990 0.0837Ra0.3082

especially for Ra ) 105. The inset plot illustrates that, for φ ) 0°, Nub and Nus are larger near the bottom corner points. Thus overall heat transfer rates are found to be higher for φ ) 0° for all Ra values. For nonuniform heating case at the bottom wall (Figure 15c), the average Nusselt numbers at Ra ) 103 for φ ) 45° are less compared to φ ) 30° and φ ) 0°. But as the Rayleigh number increases, Nub at φ ) 45° is found to be larger than that with φ ) 30° and φ ) 0° which may be due to larger Nub as explained in the right inset plot of Figure 13a. Figure 15d illustrates Nus for the nonuniform heating case, and it is observed that Nus is larger for φ ) 0° irrespective of Ra. This may be explained from right inset plot of Figure 13b. Figure 16a-d illustrates average heat transfer rates for salt water. For the uniform heating case (Figure 16a and b), variation of Nusselt number is qualitatively similar to molten metals but heat transfer rates are larger for Pr ) 7.2. For the uniform heating case, the average Nusselt numbers are found to be higher for φ ) 0° than that for φ ) 45° and 30° along the bottom and side walls of the enclosure. The average Nusselt numbers for the nonuniformly heated bottom wall (Figure 16c and d) are very close for φ ) 45°, 30°, and 0°, and the effect of the angle φ has less significance. This is also evident from Figure 14a and b and the right inset plot. The overall heat transfer rate for the side wall is larger for φ ) 0° than that with φ ) 45° and 30° for nonuniform heating cases and that is based on larger Nus values at φ ) 0° as seen from right inset plot of Figure 14b. The average Nusselt numbers show that the overall heat transfer rate decreases with an increase in angle for most of the cases. It may also be noted that the average Nusselt number increases with Ra for specific Ra in the presence of convection dominated heat transfer. The conduction dominated regime is shown as asymptotes in Figures 15 and 16. The correlations for the Nusselt number as a function of the Rayleigh number within the convection dominated regime for various angles φ ) 45°, 30°, and 0° are shown for Pr ) 0.015 and 7.2 in Tables 1 and 2, respectively.

performed for Pr ) 0.015 (molten metal), Pr ) 7.2 (salt water), and Pr ) 988.24 (olive oil) for various values of the Rayleigh number in the range 103 e Ra e 105 and side wall inclination angle φ ) 45°, 30°, and 0° (square). For Pr ) 0.015 (molten metal) and Ra ) 103, stronger convection is observed for φ ) 45° and 30° than for φ ) 0° for both uniform and nonuniform heating. At Ra ) 105, secondary circulations are obtained due to the significance of low viscous effects of molten metal. As the Prandtl number increases to Pr ) 7.2 and 988.24, the circulations become stronger and isotherms completely take the shape of the cavity with the absence of secondary circulation. For nonuniform heating of the bottom wall at Ra ) 103 for all model fluids (molten metal, salt water, and olive oil), the thermal boundary layer developed is less intense as compared to the uniform heating case within the cavity. During conduction dominant regime, variation of tilt angle from φ ) 45° to φ ) 30° has less significance. For Pr ) 0.015 and 7.2, the larger value of local Nusselt numbers were observed toward corner of the bottom wall for uniform heating especially at Ra ) 105. For nonuniform heating at Pr ) 0.015, the local Nusselt numbers at the bottom wall were found to be higher for φ ) 45° than that for φ ) 30° and 0°. As Pr increases to Pr ) 7.2, the local Nusselt number at the bottom wall is found to be higher for φ ) 0° than that for φ ) 45° and 30°. For Pr ) 0.015, the local Nusselt numbers at the side walls show decreasing trend with distance at φ ) 0° for both uniform and nonuniform heating, whereas for Pr ) 7.2, the local Nusselt numbers at φ ) 0° show increasing trend. The local Nusselt number values show increase in the heat transfer rates for salt water compared to molten metals which results due to larger convection. The average Nusselt numbers for Pr ) 0.015 are observed more for φ ) 0° than that for φ ) 45° and 30° at lower Ra whereas average nusselt number are larger for φ ) 45° at higher Ra. For all other cases the average Nusselt numbers were found to be more for φ ) 0° than 45° and 30°. Acknowledgment The authors would like to thank the anonymous reviewers for constructive comments which improved the quality of the manuscript. Appendix The name “isoparametric” derives from the fact that the same parametric function describing the geometry may be used for interpolating the spatial variable within an element. Figure 2 shows a trapezoidal domain mapping to a square domain. The transformation between (x, y) and (ξ, η) coordinates can be defined by X ) Σk9) 1 Φk(ξ, η)xk and Y ) Σk9) 1 Φk(ξ, η)yk where (xk, yk) are the X, Y coordinates of the k nodal points as seen in Figure 2a and b and Φk(ξ, η) is the basis function. The nine basis functions are the following: Φ1 ) (1 - 3ξ + 2ξ2)(1 - 3η + 2η2)

6. Conclusion

Φ2 ) (1 - 3ξ + 2ξ2)(4η - 4η2)

The role of uniform and nonuniform heating of the bottom wall and heat transfer characteristics due to natural convection heating of liquid food substances or molten materials in the trapezoidal enclosure has been studied in detail. The penalty finite element method has been used, and smooth solutions are obtained in terms of stream function and isotherm contours for wide ranges of Pr and Ra. Numerical simulations were

Φ3 ) (1 - 3ξ + 2ξ2)(-η + 2η2) Φ4 ) (4ξ - 4ξ2)(1 - 3η + 2η2) Φ5 ) (4ξ - 4ξ2)(4η - 4η2) Φ6 ) (4ξ - 4ξ2)(-η + 2η2)

Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 8665

Φ7 ) (-ξ + 2ξ )(1 - 3η + 2η ) 2

2

Φ8 ) (-ξ + 2ξ2)(4η - 4η2) Φ9 ) (-ξ + 2ξ2)(-η + 2η2) The above basis functions are used for mapping the trapezoidal domain into square domain and the evaluation of integrals of residuals [eqs 16, 18, and 21]. Nomenclature g ) acceleration due to gravity, m s-2 k ) thermal conductivity, W m-1 K-1 H ) height of the trapezoidal cavity Nu ) local Nusselt number Nu ) average Nusselt number p ) pressure, kg m-1 s-2 P ) dimensionless pressure Pr ) Prandtl number Ra ) Rayleigh number T ) temperature, K Th ) temperature of hot bottom wall, K Tc ) temperature of cold vertical wall, K 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 X ) dimensionless distance along x coordinate Y ) dimensionless distance along y coordinate Greek Symbols R ) thermal diffusivity, m2 s-1 β ) volume expansion coefficient, K-1 γ ) penalty parameter θ ) dimensionless temperature ν ) kinematic viscosity, m2 s-1 φ ) angle of inclination of the left wall F ) density, kg m-3 Φ ) basis functions ψ ) stream function ξ ) horizontal coordinate in a trapezoidal η ) vertical coordinate in a trapezoidal Subscripts b ) bottom wall s ) side wall

Literature Cited (1) Akterian, S. G. Numerical simulation of unsteady heat conduction in arbitrary shaped canned foods during sterilization processes. J. Food Eng. 1994, 21, 343–354. (2) Park, R. J.; Kim, S. B.; Kim, H. D.; Choi, S. M. Natural convection heat transfer with crust formation in the molten metal pool. Nuclear Technol. 1999, 127, 66–80. (3) Ulir´, J. Chemistry and technology of molten salts reactors - history and perspectives. J. Nucl. Mater. 2007, 360, 6–11. (4) Basak, T.; Ayappa, K. G. Influence of internal convection during microwave thawing of cylinders. AIChE J. 2001, 47, 835–850. (5) Ousegui, A.; Le Bail, A.; Havet, M. Numerical and experimental study of a natural convection thawing process. AIChE J. 2006, 52, 4240– 4247. (6) Walle, E.; Finne, J.; Picard, G.; Sanchez, S.; Boursier, J. M.; Noel, D. A computational tool for selective pyrochemical processes based on molten salts in nuclear industry. J. Nucl. Mater. 2005, 344, 158–164. (7) Moriyama, H.; Sagara, A.; Tanaka, S.; Moir, R. W.; Sze, D. K. Molten salts in fusion nuclear technology. Fusion Eng. Des. 1998, 3-940, 627–637. (8) Mishra, B.; Olson, D. L. Molten salt applications in materials processing. J. Phys. Chem. Sol. 2005, 66, 396–401.

(9) Kumar, A.; Bhattacharya, M.; Blaylock, J. Numerical simulation of natural convection heating of canned thick viscous liquid food products. J. Food Sci. 1990, 55, 1403–1411. (10) Bergman, T. L.; Hyun, M. T. Simulation of two-dimensional thermosolutal convection in liquid metals induced by horizontal temperature and species gradients. Int. J. Heat Mass Transfer 1996, 39, 2883–2894. (11) Xu, B.; Li, B. Q.; Stock, D. E. An experimental study of thermally induced convection of molten gallium in magnetic fields. Int. J. Heat Mass Transfer 2006, 49, 2009–2019. (12) Pelekh, A. E.; Mukasyan, A. S.; Varma, A. Kinetics of rapid hightemperature reactions: Titanium-nitrogen system. Ind. Eng. Chem. Res. 1999, 38, 793–798. (13) Shang, H.; Kingman, S. W.; Snape, C. E.; Robinson, J. P. Reactors effects on microwave decontamination of oily wastes in a multimode cavity. Ind. Eng. Chem. Res. 2007, 46, 4811–4818. (14) Jung, A.; Fryer, P. J. Optimising the quality of safe food: computational modelling of a continuous sterilisation Process. Chem. Eng. Sci. 1999, 54, 717–730. (15) Zhang, L.; Fryer, P. J. Models for the electrical heating of solid liquid food mixtures. Chem. Eng. Sci. 1993, 48, 633–642. (16) De Alwis, L.; Fryer, P. J. A Finite element analysis of heat generation and transfer during ohmic heating of food. Chem. Eng. Sci. 1990, 45, 1547–1559. (17) Farid, M. M.; Ghani, A. G. A. A new computational technique for the estimation of sterilization time in canned food. Chem. Eng. Process. 2004, 43, 523–531. (18) Laguerre, O.; Flick, D. Heat transfer by natural convection in domestic refrigerators. J. Food Eng. 2004, 62, 79–88. (19) Moraga, N. O.; Herna´n Barraza, G. Predicting heat conduction during solidification of food inside a freezer due to natural convection. J. Food Eng. 2003, 56, 17–26. (20) Moraga, N. O.; Salinas, C. H. Numerical study of unsteady 2D natural convection and solidification of a food inside a freezing chamber. Num. Heat Transfer Part A 2000, 37, 755–777. (21) Ghani, A. G. A.; Farid, M. M.; Chen, X. D.; Richards, P. Numerical simulation of natural convection heating of canned food by computational fluid dynamics. J. Food Eng. 1999, 41, 55–64. (22) Ghani, A. G. A.; Farid, M. M.; Zarrouk, S. J. The effect of can rotation on sterilization of liquid food using computational fluid dynamics. J. Food Eng. 2003, 57, 9–16. (23) Weerts, A. H.; Lian, G.; Martin, D. R. Modeling the hydration of foodstuffs: Temperature effects. AIChE J. 2003, 49, 1334–1339. (24) Chen, X. D.; Lin, S. X. Q. Air drying of milk droplet under constant and time-dependent conditions. AIChE J. 2005, 51, 1790–1799. (25) Sovova´, H.; Kuc´era, J; Jez´, J. Rate of the vegetable oil extraction with supercritical CO2sII. Extraction of grape oil. Chem. Eng. Sci. 1994, 49, 415–420. (26) Garrpeters, J. M. The neutral stability of surface-tension driven cavity flows subject to buoyant forces. 1. transverse and longitudinal disturbances. Chem. Eng. Sci. 1992, 47, 1247–1264. (27) Bejan A. ConVection Heat Transfer, 3rd ed.; Wiley: Hoboken, NJ, 2004. (28) Gebhart, B. Buoyancy induced fluid motions characteristics of applications technology: The 1978 Freeman Scholar Lecture. J. Fluids Eng. 1979, 101, 5–28. (29) Hoogendoorn C. J. Natural convection in enclosures. In Proceedings of the Eighth International Heat Transfer Conference, San Francisco; Hemisphere Publishing Corp.: Bristol, PA, 1986; Vol. I, pp 111-120. (30) Patterson, J.; Imberger, J. Unsteady natural convection in a rectangular cavity. J. Fluid Mech. 1980, 100, 65–86. (31) Nicolette, V. F.; Yang, K. T.; Lloyd, J. R. Transient cooling by natural convection in a two-dimensional square enclosure. Int. J. Heat Mass Transfer 1985, 28, 1721–1732. (32) Hall, J. D.; Bejan, A.; Chaddock, J. B. Transient natural convection in a rectangular enclosure with one heated side wall. Int. J. Heat Fluid Flow 1988, 9, 396–404. (33) Hyun, J. M.; Lee, J. W. Numerical solutions of transient natural convection in a square cavity with different sidewall temperature. Int. J. Heat Fluid Flow 1989, 10, 146–151. (34) Fusegi, T.; Hyun, J. M.; Kuwahara, K. Natural convection in a differentially heated square cavity with internal heat generation. Num. Heat Transfer Part A: Appl. 1992, 21, 215–229. (35) Lage, J. L.; Bejan, A. The Ra-Pr domain of laminar natural convection in an enclosure heated from the side. Num. Heat Transfer Part A: Appl. 1991, 19, 21–41. (36) Lage, J. L.; Bejan, A. The resonance of natural convection in an enclosure heated periodically from the side. Int. J. Heat Mass Transfer 1993, 36, 2027–2038.

8666 Ind. Eng. Chem. Res., Vol. 47, No. 22, 2008 (37) Xia, C.; Murthy, J. Y. Buoyancy-driven flow transitions in deep cavities heated from below. J. Heat Transfer 2002, 124, 650–659. (38) Iyican, L.; Bayazitoglu, Y. An analytical study of natural convective heat transfer within trapezoidal enclosure. J. Heat Transfer 1980, 102, 640–647. (39) Karyakin, Y. U. E. Transient natural convection in prismatic enclosures of arbitrary cross-section. Int. J. Heat Mass Transfer 1989, 32, 1095–1103. (40) Peric´, M. Natural convection in trapezoidal cavities. Num. Heat Transfer Part A: Appl. 1993, 24, 213–219. (41) Kuyper, R. A.; Hoogendoorn, C. J. Laminar Natural convection flow in trapezoidal enclosures. Num. Heat Transfer Part A: Appl. 1995, 28, 55–67. (42) Boussaid, M.; Djerrada, A.; Bouhadef, M. Thermosolutal transfer within trapezoidal cavity. Num. Heat Transfer Part A: Appl. 2003, 43, 431– 448. (43) Reddy, J. N.; Gartling, D. K. The Finite Element Method in Heat Transfer and Fluid Dynamics; CRC Press: Boca Raton, FL, 1993. (44) Roy, S.; Basak, T. Finite element analysis of natural convection flows in a square cavity with non-uniformly heated wall(s). Int. J. Eng. Sci. 2005, 43, 668–680. (45) Basak, T.; Roy, S.; Balakrishnan, A. R. Effects of thermal boundary conditions on natural convection flows within a square cavity. Int. J. Heat Mass Transfer 2006, 49, 4525–4535.

(46) Sarris, I. E.; Lekakis, I.; Vlachos, N. S. Natural convection in a 2D enclosure with sinusoidal upper wall temperature. Num. Heat Transfer Part A: Appl. 2002, 42 (5), 513–530. (47) Oden, J. T.; Reddy, J. N. An Introduction to the Mathematical Theory of Finite Elements; John Wiley & Sons: New York, 1976. (48) Zienkiewicz, O. C.; Taylor, R. L.; Too, J. M. Reduced integration technique in general analysis of plates and shells. Int. J. Num. Meth. Eng. 1971, 3, 275–290. (49) Batchelor, G. K. An Introduction to Fluid Dynamics; Cambridge University Press: New York, 1993. (50) Rahman, R. Food properties handbook; John Wiley & Sons: New York, 1995. (51) Ganzarolli, M. M.; Milanez, L. F. Natural convection in rectangular enclosures heated from below and symmetrically cooled from the sides. Int. J. Heat Mass Transfer 1995, 38, 1063–1073.

ReceiVed for reView February 15, 2008 ReVised manuscript receiVed July 14, 2008 Accepted August 13, 2008 IE800263C