Back to Advanced Heat and Mass Transfer Home
Table of Contents
6 NATURAL CONVECTION 6.1 Introduction Contrary to forced convection, which is driven by the external means, natural convection is induced by a density difference or density gradient in the flow field. This density difference may occur due to influence of temperature gradients on the density field, or due to differences in concentration or concentration gradients. Another condition that is necessary for natural convection to occur is body force that is proportional to the density. The most common body force is gravitational force, but the body force resulting from centrifugal or Coriolis forces is also possible. The combined body force and the density difference or gradient will produce buoyancy force within the fluid, which is the driving force behind natural convection. Natural convection can find its application in electronics cooling, the environmental science, manufacturing and materials processing, energy systems, and fire safety, etc. While the flow field for most forced convection problems can be solved independently of the solution for the temperature distribution, the flow field and temperature distribution in natural convection problems cannot be solved separately. A special technique must be developed to handle the coupling between the flow and heat transfer. Most natural convection problems that one encounters are due to density differences caused by temperature differences. The densities of most liquids and gases decrease with increasing temperature, i.e., ∂ρ / ∂T < 0 (with a few exceptions, such as water between 0°C and 4°C). So as to understand the mechanism of natural convection, let us consider the natural convection from a heated vertical flat plate placed in a quiescent fluid. The surface temperature of the flat plate, Tw, is higher than the temperature of the quiescent fluid, T∞. The fluid that is in contact with the heated vertical plate is heated and its density decreases ( ∂ρ / ∂T < 0 ). The heated lighter fluid flows up and the fluid from the neighbor moves in. The temperature in a thin layer of the fluid near the heated vertical plate is higher than the quiescent fluid temperature T∞ and a thermal boundary layer is formed. As the fluid near the heated vertical plate rises, a velocity boundary layer also forms. Depending on the Prandtl number of the fluid, the velocity boundary layer can be thicker (Pr > 1; see Fig. 6.1) or thinner (Pr < 1) than the thermal boundary layer. The distributions of temperature and velocity in the boundary layers are also sketched in Fig. 6.1. Depending on the Chapter 6 Natural Convection 515 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Thermal boundary x layer, δt Velocity boundary layer, δ Tw T u quiescent fluid, u∞=0 T∞, ρ∞ gravity y Figure 6.1 Natural convection over a vertical flat plate (Pr >1) temperature difference, the properties of the fluid, and the length of the flat plate, a transition from laminar to turbulent flow, similar to Fig. 5.3, can also occur. The local heat flux at the flat surface can be described using Newton’s law of cooling q′′ = hx (Tw − T∞ ) (6.1) where hx is the convective heat transfer coefficient, which depends on flow configurations and the properties of the fluid. The heat transfer coefficient usually also depends on the temperature difference, Tw − T∞ , which means that the heat flux is not a linear function of the temperature difference. At the surface of the heated plate, the velocity of the fluid is zero due to the non-slip condition. Therefore, conduction is the mechanism of heat transfer from the heated plate to the fluid immediately in contact with the heated plate so the heat transfer can be described by Fourier’s law of conduction:  ∂T  q′′ = −k    ∂y  y = 0 (6.2) Combining eqs. (6.1) and (6.2) yields hx = − k Tw − T∞  ∂T     ∂y  y = 0 (6.3) which can be used to determine the local heat transfer coefficient. It reveals that the temperature profile in the thermal boundary layer is required to determine the heat transfer coefficient. 516 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press While the above discussions are for the case wherein the temperature of the vertical flat plate is higher than that of the quiescent fluid, they are also valid for the case wherein the temperature of the vertical flat plate is lower than that of the quiescent fluid. The only difference in the latter case is that the flow near the vertical plate will be directed downward and the boundary layers will start to form from the top of the flat plate. Natural convection can also occur in an enclosure. In this case the development of the boundary layer is restricted by the wall of the enclosure. Natural convection in an enclosure is referred to as internal convection, which is different from the external natural convection from a heated or cooled vertical plate. The examples of internal natural convection include buoyancy-induced flow in a room caused by a heater or air-conditioner, electronics cooling in the computer chassis, melting or solidification of phase change materials in an enclosure, and propagation of fire in a room. Both internal and external natural convections are the subjects of discussion in this chapter. Section 6.2 presents the governing equations for natural convection, which are obtained by simplifying the generalized governing equations presented in Chapter 2, followed by scale analysis of natural convection in Section 6.3. Section 6.4 discusses external natural convection from vertical, inclined, and horizontal flat plates, cylinders, and spheres, as well as free flows including plumes, and wakes. Internal natural convection is the subject of Section 6.5. Natural convection in melting and solidification is discussed in Section 6.6. This chapter is closed by discussion of instability analysis of natural convection. 6.2 Governing Equations for Natural Convection 6.2.1 Generalized Governing Equations The governing equations for natural convection are special cases of the generalized governing equations that were discussed in Chapter 2. Consider a multicomponent system with N components, the governing equations for a stationary reference frame can be expressed as Dρ + ρ∇ ⋅ V = 0 (6.4) Dt DV ρ = ρ g − ∇p + ∇ ⋅ τ Dt (6.5) where the viscous stress tensor, τ, can be determined by using Newton’s law of viscosity [see eq. (1.53)]: 2 τ = 2μ D − μ (∇ ⋅ V )I 3 (6.6) where D is the rate of strain tensor: Chapter 6 Natural Convection 517 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press D= 1 T ∇ V + ( ∇V )   2 (6.7) For natural convection problems, it is often assumed that the fluid is incompressible, except in the first term on the right-hand side of eq. (6.5); this is referred to as the Boussinesq assumption. Under this assumption, the continuity equation (6.4) becomes: ∇⋅V = 0 (6.8) According to eq. (6.8), the second term on the right-hand side of eq. (6.6) will be zero. The momentum equation (6.5) then becomes: ρ DV = ρ g − ∇p + ∇ ⋅ ( μ∇V ) Dt (6.9) where the left-hand side is the inertial term; the three terms on the right-hand side represent body force per unit volume, pressure force per unit volume, and viscous force per unit volume, respectively. The density of a mixture is a function of its temperature and the mass fractions of its species. It can be expanded using a Taylor’s series near the vicinity of a reference point ( T∞ , ω1,∞ , ω2,∞ , ω N ,∞ ): ρ = ρ∞ + N ∂ρ ∂ρ (T − T∞ ) +  (ωi −ωi ,∞ ) +  ∂T ∂ωi i =1 where ρ∞ is density at the reference point. By defining the coefficient of thermal expansion, β, and composition coefficient of volume expansion, βm,i, as follows: 1  ∂ρ  β =−  (6.10) ρ∞  ∂T  p  β m,i = − 1  ∂ρ    ρ∞  ∂ωi  p (6.11) and neglecting the higher order terms in the Taylor’s series expansion, one obtains: ρ ≈ ρ∞ − ρ ∞ β (T − T∞ ) − ρ ∞  β m (ωi −ωi ,∞ ) i =1 N (6.12) which is valid only if β (T − T∞ )  1 and β m (ωi − ωi ,∞ )  1 . Substituting eq. (6.12) into eq. (6.9), the momentum equation for natural convection is obtained: ρ DV = ( −∇p + ρ∞ g ) − ρ∞ gβ (T − T∞ ) Dt − ρ ∞ g  β m ,i (ωi −ωi ,∞ ) + ∇ ⋅ ( μ∇V ) i =1 N (6.13) which is a generalized momentum equation because the effects of buoyancy forces due to both temperature and composition variations are considered. If it is assumed that the Dufour effect is negligible and the fluid is incompressible, the energy equation is: 518 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ρcp The conservation of species mass in terms of mass fraction for the ith species can be expressed as Dωi  ρ (6.15) = −∇ ⋅ J i + mi′′′ Dt DT = ∇ ⋅ (k∇T ) + q′′′ + ∇V : τ Dt (6.14) For a binary system of A and B, one can apply Fick’s law to eq. (6.15) to obtain: Dω A A ρ = − ρ∇ ⋅ ( D AB ∇ω A ) + m′′′ (6.16) Dt 6.2.2 External Natural Convection from Heated Vertical Plate For external natural convection near a vertical flat plate as shown in Fig. 6.1, the boundary layer assumption can be applied to simplify the above generalized governing equations. The boundary layer treatment for the case of natural convection is very similar to that for the case of forced convection that was discussed in Chapter 4. The difference between the natural convection problem shown in Fig. 6.1 and forced convection over a flat plate is that the free stream velocity in the outside of the velocity boundary layer in natural convection is zero for natural convection. In addition, the pressure outside the boundary layer is hydrostatic for the case of natural convection, instead of being externally imposed as in the case of forced convection. For 2-D external convection of an incompressible fluid as shown in Fig. 6.1, the continuity equation becomes ∂u ∂v + =0 ∂x ∂y (6.17) If one assumes that the fluid is single component so that the natural convection is driven by the density difference induced by the temperature gradient, eq. (6.13) becomes: ρ DV = ( −∇p + ρ∞ g ) − ρ∞ gβ (T − T∞ ) + ∇ ⋅ ( μ∇V ) Dt (6.18) Applying the boundary layer assumption and assuming steady flow with constant thermophysical properties, the momentum equation becomes: u ∂u ∂u 1 ∂p ∂2u +v =− − g + gβ (T − T∞ ) + ν 2 ρ∞ ∂x ∂x ∂y ∂y (6.19) The pressure in the boundary layer, p, is independent of y ( ∂p / ∂y = 0 ) and equals that outside the boundary layer at the same longitudinal position, p∞, i.e., ∂p dp dp∞ = = ∂x dx dx Chapter 6 Natural Convection 519 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The hydrostatic pressure, p∞, is dictated by the density and the longitudinal position: dp∞ = − ρ∞ g dx Substituting the above two equations into eq. (6.19), the momentum equation becomes: u ∂u ∂u ∂2u +v = ν 2 + gβ (T − T∞ ) ∂x ∂y ∂y (6.20) After applying the boundary layer assumption and assuming the viscous dissipation is negligible, the energy equation becomes: u ∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y (6.21) At the heated wall, the non-slip and impermeable conditions yield the following boundary condition for the momentum equation: u = v = 0, at y = 0 (6.22) The temperature at the heated wall is specified, i.e., T = Tw , at y = 0 (6.23) Since the quiescent fluid far away from the heated plate is not disturbed by the existence of the heated plate, the velocity at the locations away from the flat plate should be zero: u = v = 0, y → ∞ (6.24) Also, the temperature of the fluid outside the thermal boundary layer is not affected by the heated wall: T = T∞ , y → ∞ (6.25) 6.2.3 Dimensionless Parameters While the Reynolds number was used as a dimensionless parameter to characterize the flow for the case of forced convection, it cannot be used to characterize the natural convection because the characteristic velocity is not available. To identify the appropriate dimensionless parameters for description of natural convection, one defines the following dimensionless variables: X= T − T∞ x y u v ,Y= ,U= ,V= ,θ= L L u0 u0 Tw − T∞ (6.26) where u0 is a reference velocity that is unknown at this point. Equations (6.17), (6.20) and (6.21) can be respectively nondimensionalized as follows (Incropera et al., 2007): ∂U ∂V + =0 ∂X ∂Y 1 ∂ 2U gβ (Tw − T∞ ) L ∂U ∂U U +V = + θ 2 u0 ∂X ∂Y Re L ∂Y 2 (6.27) (6.28) 520 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U ∂θ ∂θ 1 ∂ 2θ +V = ∂X ∂Y Re L Pr ∂Y 2 Re L = u0 L (6.29) where (6.30) ν is the Reynolds number based on the reference velocity. To simplify the dimensionless governing equations, one can choose the reference velocity to be: u0 = gβ (Tw − T∞ ) L (6.31) so that the factor before dimensionless temperature, θ, in eq. (6.28) becomes unity. With this choice of reference velocity, the Reynolds number becomes: gβ (Tw − T∞ ) L3 Re L = (6.32) ν In natural convection problems, we define Grashof number, GrL, to be gβ (Tw − T∞ ) L3 GrL = Re 2 = (6.33) L ν2 as the dimensionless number which represents the ratio of the buoyancy force and the viscosity force acting on the fluid. Equations (6.28) and (6.29) then become: U U ∂U ∂U 1 ∂ 2U +V = 1/2 +θ ∂X ∂Y GrL ∂Y 2 ∂θ ∂θ 1 ∂ 2θ +V = 1/2 ∂X ∂Y GrL Pr ∂Y 2 (6.34) (6.35) The role of the Grashof number for a natural convection problem is similar to the role of Reynolds number for a forced convection problem. The Prandtl number, which reflects the ratio between momentum and thermal diffusions, is another parameter that affects natural convection. Therefore, it is expected that the average Nusselt number for natural convection, Nu L = hL / k , will be a function of the Grashof number and the Prandtl number: Nu L = f (GrL , Pr) (6.36) The objective of this chapter is thus to identify the appropriate forms of the above function for various geometric configurations and physical conditions. 6.3 Scale Analysis While scale analysis cannot provide the exact functions in eq. (6.36), the form of these functions can be provided (Bejan, 2004). Let us refer to Fig. 6.1 and consider the governing equations in the thermal boundary layer ( y  δ t ) for the entire flat plate ( x  L ). The thickness of the thermal boundary layer in which the effects of the heated wall are felt is much smaller than the length of the Chapter 6 Natural Convection 521 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press vertical plate, i.e. δ t  L . For the continuity equation (6.17) to be satisfied, the scales of the two terms must be the same: uv  L δt The scale of the velocity component in the y-direction is therefore: δt v L u (6.37) which indicates that v  u for flow in the thermal boundary layer. The scale of the velocity component in the x-direction, u, is still unknown at this point. While the left hand side of the energy equation (6.21) shows the effect of advection, the right-hand side shows the effect of diffusion. The scales of the two advective terms on the left hand side of eq. (6.21) are: u ∂T ΔT ∂T ΔT ΔT u ,v v u L δt L ∂x ∂y which indicates that the scale of the second term on the left hand side of eq. (6.21) is identical to the scale of the first term when the scale of the velocity component in the y-direction is given by eq. (6.37). The temperature difference ΔT = Tw − T∞ in the above scale analysis represents the scale of the excess temperature, T − T∞ . The scale of the right-hand side of eq. (6.21) is: α ∂ 2T ΔT α 2 2 ∂y δt ΔT ΔT α 2 L δt u α L The scales of the two sides of the energy equation must be the same: u The above equation can be used to estimate the scale of u as follows: (6.38) δ t2 The scale of v can be obtained by substituting eq. (6.38) into eq. (6.37): α v (6.39) δt where the scale of the thermal boundary layer thickness, δt, is still unknown at this point. The respective scales of the two inertial terms on the left-hand side of the momentum equation (6.20) are: ∂u u 2  ∂x L ∂u u u2 v v  L δt ∂y u The respective scales of the viscosity and buoyancy terms on the right-hand side of eq. (6.20) are: 522 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ν ∂2u u ν 2 2 δt ∂y gβ (T − T∞ )  gβ ΔT We can very well see that three forces are at play in the boundary layer region, which are inertia, viscosity and buoyancy forces. Considering the scale of u obtained from eq. (6.38), the scales of these three forces are α 2L να L g β ΔT , , (6.40) δ t4 δ t4 Inertia Viscous Buoyancy Among these three forces, the buoyancy force is never negligible because without it natural convection would not occur. Therefore, the scale of buoyancy force can be used to measure the importance of the inertial and viscous forces. Dividing eq. (6.40) by the scale of buoyancy force, gβΔT , one obtains the following: α2  L  να  L  ,  , 3 gβ ΔTL  δ t  gβΔTL3  δ t  Inertia 4 4 4 1 Buoyancy Viscous 4 which can be expressed in terms of dimensionless parameters as (Bejan, 2004): L L −1 −1 −1 1   Ra L Pr ,   Ra L , δt  δt    Inertia Viscous Buoyancy (6.41) where (6.42) να is the Rayleigh number that is related to the Grashof number by Ra L = GrL Pr . Therefore, the relative importance of inertia and viscous forces depends on the Prandtl number, which is a property of the fluid. Thus, if the Prandtl number is high ( Pr  1 ), the inertia term will be negligible and the viscosity term will balance the buoyancy term, whereas if the Prandtl number is low enough ( Pr  1 ) – as for liquid metals – then the inertia term is considerable and balances the buoyancy term in steady state. The scale analysis for high- and lowPrandtl number fluids is presented in detail below. Ra L = gβ (Tw − T∞ ) L3 6.3.1 High Prandtl Number Fluids (Pr  1) For the case of high Prandtl number fluids, the kinematic viscosity of the fluid is much greater than its thermal diffusivity, so viscous force has a significant effect on the momentum equation and the effect of inertial force is negligible. The viscous force in fact balances the buoyancy force, and so from eq. (6.41), one obtains the following: Chapter 6 Natural Convection 523 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ t  L Ra −1/ 4 L Substituting eq. (6.43) into eqs. the velocity components become: α 1/ 2 u v L Ra L (6.38) and (6.39), the (6.43) scales of (6.44) (6.45) α L Ra1/ 4 L The scale of the heat transfer coefficient can be obtained by analyzing eq. (6.3), as follows:  ∂T  −k    ∂y  y = 0 k (ΔT / δ t ) k  = hx = ΔT Tw − T∞ δt The scale of Nusselt number is therefore: Nu = hx L L   Ra1/ 4 L δt k (6.46) which will be confirmed by exact solution in the next section. For the case of high Prandtl number fluids the momentum boundary layer thickness, δ, is much greater than the thermal boundary layer thickness, δt (see Fig. 6.1). While the above scale analysis indicated that flow within the thermal boundary layer is dominated by balance between the buoyancy force and the viscous force, buoyancy force does not exist beyond the thermal boundary layer. Therefore, the velocity that is developed outside the thermal boundary layer but within the momentum boundary layer results from the viscous drag of the thermal boundary layer. Therefore, the flow between the thermal and momentum boundary layers is dominated by the balance between the viscous force and the inertial force. Since the thermal boundary layer is very thin for the case of Pr  1 , there exists a balance between inertia and viscosity in the momentum equation over the entire momentum boundary layer. Refer to Fig. 6.1 once again and consider the momentum equation of the momentum boundary layer ( y  δ ) for the entire flat plate ( x  L ). The thickness of the momentum boundary layer is much smaller than the length of the vertical plate, i.e. δ  L . The force balance in the momentum boundary layer between inertial and viscous forces gives us: u2 u ν 2 L δ where the vertical velocity component, u, is induced by the thin thermal boundary layer. Substituting eq. (6.44) into the above equation yields the scale of the momentum boundary layer: δ  LRa −1/ 4 Pr1/ 2 (6.47) L Comparing eqs. (6.47) and (6.43) gives us the following relationship between the thicknesses of the momentum and thermal boundary layers: 524 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ  Pr1/ 2 > 1 (6.48) δt It is evident that the greater the Prandtl number, the greater the ratio of δ/δt. This means that the region of unheated fluid which is being driven vertically by the heated layer due to viscous action is thicker for fluids with a higher Prandtl number. 6.3.2 Low Prandtl Number Fluids (Pr  1) In the case of low Prandtl number, eq. (6.41) indicates that inertial force has a significant effect on the momentum equation, while the effect of viscous force is negligible. The momentum equation requires that the inertia term to be balanced by the buoyancy term, i.e., δ t  L(Ra L Pr) −1/ 4 (6.49) Substituting eq. (6.49) into eqs. (6.38) and (6.39), the scales of the velocity components become: α u  (Ra L Pr)1/ 2 (6.50) L v α L (Ra L Pr)1/ 4 Momentum and thermal boundary layers, δ = δt (6.51) x shear layer, δu Tw T u quiescent fluid, u∞=0 T∞, ρ∞ gravity y Figure 6.2 Natural convection over a vertical flat plate (Pr < 1) Chapter 6 Natural Convection 525 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The scale of the Nusselt number is:  (Ra L Pr)1/ 4 (6.52) δt Since the Prandtl number is much less than 1, the thermal boundary layer is driven upwards by buoyancy and is restrained by inertia. Beyond the thermal boundary layer region, where there is no heating effect and consequently no buoyancy effect, the fluid is at rest and there is no motion whatsoever there. Therefore, the velocity boundary layer, δ, must end with the thermal boundary layer, δt. Due to the no slip boundary condition at the wall there should be a velocity peak somewhere between the wall and the thermal boundary layer. Consider a thickness of the velocity boundary layer from the wall up to the point of maximum velocity, δ u . The effect of inertia is negligible within this thin layer δ u – referred to as the shear layer (Bejan, 2004) – and the buoyancy force is balanced by the viscous force: Nu   g β ΔT δ u2 Substituting eq. (6.50) into the above equation, the scale of the shear layer can be determined as − δ u  LGrL 1/ 4 (6.53) The ratio between the thicknesses of the shear layer and the thermal boundary layer can be obtained by dividing eq. (6.53) by eq. (6.49): δu  Pr1/ 2 < 1 (6.54) δt It should be pointed out that the thickness of the shear layer, δu, is different from the momentum boundary layer, δ, which is approximately equal to the thermal boundary layer thickness for the case where the Prandtl number is less than 1. L ν u 6.4 External Natural Convection Laminar external natural convection on vertical and horizontal flat plates, over horizontal and vertical cylinders and spheres, as well as plumes, wakes and other types of free flow will be discussed in this section. 6.4.1 Similarity Solution for Natural Convection on a Vertical Surface The boundary layer-type governing equations for the external convection problem shown in Fig. 6.1 are eqs. (6.17), (6.20) and (6.21). Introducing the stream function ψ , and dimensionless temperature, θ : T − T∞ ∂ψ ∂ψ u= , v=− ,θ= (6.55) ∂y ∂x Tw − T∞ 526 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the continuity equation (6.17) is satisfied and the momentum and energy equations (6.20) and (6.21) become: ∂ψ ∂ 2ψ ∂ψ ∂ 2ψ ∂ 3ψ − = ν 3 + gβ (Tw − T∞ )θ (6.56) 2 ∂y ∂y∂x ∂x ∂y ∂y ∂ψ ∂θ ∂ψ ∂θ ∂ 2θ − =α 2 ∂y ∂x ∂x ∂y ∂y (6.57) with the following boundary conditions: ∂ψ ∂ψ = = 0, θ = 1 at y = 0 ∂y ∂x ∂ψ ∂ψ = = 0, θ = 0 ∂y ∂x at y → ∞ (6.58) (6.59) The idea behind the similarity solution is that the velocity and temperature profiles in different x in the boundary layers are geometrically similar, differing only by a stretching factor in the x-direction (Kays et al., 2005). Thus, the similarity variable should have the following form: η = y ⋅ H (x) (6.60) where H(x) is an unspecified stretching function. The objective now is to reduce eqs. (6.56) and (6.57) to ordinary differential equations. The stream function, ψ , which is a function of x and y, can be expressed as a function of x and η . If the similarity solution exists, one can express the stream function as ψ ( x ,η ) = ν F (η ) ⋅ G ( x ) (6.61) where F (η ) and G ( x ) are the similarity function and the stretching function, respectively. The dimensionless temperature can be assumed as a function of η only, i.e., θ ( x ,η ) = θ (η ) (6.62) The derivatives in eqs. (6.60)–(6.62) can be obtained as shown below ∂ψ ∂η = ν F ′(η ) G ( x ) = ν F ′(η ) H ( x )G ( x ) ∂y ∂y ∂ψ = ν F ′(η ) yH ′( x )G ( x ) + ν F (η ) ⋅ G ′( x ) ∂x ∂ 2ψ = ν F ′(η )[ H ′( x )G ( x ) + H ( x )G ′( x )] + ν F ′′(η )η H ′( x )G ( x ) ∂y∂x ∂ 2ψ = ν F ′′(η ) ⋅ H 2 ( x ) ⋅ G ( x ) ∂y 2 ∂ 3ψ = ν F ′′′(η ) ⋅ H 3 ( x ) ⋅ G ( x ) ∂y3 ∂θ ∂η = θ ′(η ) = θ ′(η ) yH ′( x ) ∂x ∂x ∂θ ∂η = θ ′(η ) = θ ′(η ) H ( x ) ∂y ∂y Chapter 6 Natural Convection 527 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where the primes for F and θ denote the derivatives with respect to η , while the primes for G and H denote the derivatives with respect to x. Substituting the above derivatives into eqs. (6.56) and (6.57) and considering eq. (6.60), one obtains: gβ (Tw − T∞ ) G′  H ′G G ′  F ′′′ + θ + FF ′′ −  2 +  ( F ′) 2 = 0 (6.63) 2 3 H H ν HG H θ ′′ + Pr G′ Fθ ′ = 0 H (6.64) In order to convert eqs. (6.63) and (6.64) to ordinary differential equations with η as the sole independent variable, all functions of x must be cancelled, which is possible only if the following combination of H and G and their derivatives are satisfied: H 3G = A = const (6.65) G′ = B = const H H ′G = C = const H2 (6.66) (6.67) Differentiation of eq. (6.65) and division of the resultant equation by H4 yields the following equation: 3 H ′G G ′ + =0 H2 H (6.68) which is satisfied if eqs. (6.66) and (6.67) are satisfied. Substituting eq. (6.68) into eq. (6.63), we have gβ (Tw − T∞ ) G′ 2 G′ F ′′′ + ( F ′) 2 = 0 (6.69) θ + FF ′′ − H 3H ν 2 H 3G Therefore, satisfaction of eqs. (6.65) and (6.66) is sufficient to ensure that eqs. (6.63) and (6.64) become ordinary differential equations. Although any constants A and B in eqs. (6.65) and (6.66) will transform eqs. (6.69) and (6.64) into ordinary differential equations, the proper choice of values of these two constants will yield ordinary differential equations with simple forms. Let us choose gβ (Tw − T∞ ) =1 ν 2 H 3G G′ =3 H The corresponding G and H functions then become: 1  G ( x ) = 4  Grx  4  H (x) = 1/ 4 (6.70) 1/ 4 11   Grx  x4  (6.71) where Grx is the local Grashof number defined as: 528 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (6.72) ν2 which is equivalent to the square of the Reynolds number based on the scale of the local velocity u0 = gβ (Tw − T∞ ) x . Grx = gβ (Tw − T∞ ) x 3 Substituting eqs. (6.70) and (6.71) into eqs. (6.69) and (6.64), the momentum and the energy equations become: F ′′′ + θ + 3FF ′′ − 2( F ′) 2 = 0 (6.73) θ ′′ + 3 Pr F θ ′ = 0 (6.74) The velocity components in the x- and y-directions can be expressed in terms of similarity variables by the following equations: ∂ψ ∂ψ ∂η 2ν = = u= Grx1/ 2 F ′(η ) (6.75) ∂y ∂η ∂y x v=− ∂ψ ∂ψ ∂η ν  Grx  [η F ′(η ) − 3F (η ) x1/ 4 ] =− = ∂x ∂η ∂x x  4    1/ 4 (6.76) which are obtained from eqs. (6.60), (6.61), (6.70) and (6.71). Substituting the above expressions into eqs. (6.58) and (6.59), the boundary conditions for eqs. (6.73) and (6.74) can be obtained as follows F (η ) = F ′(η ) = 0 and θ (η ) = 1 at η = 0 (6.77) F ′(η ) = θ (η ) = 0 at η → ∞ (6.78) Equations (6.73) and (6.74) are coupled nonlinear ordinary equations with boundary conditions specified at different η and they thus make external natural convection a boundary value problem. The original numerical solution was obtained by Ostrach (1953) for a wide range of Prandtl numbers from 0.01 to 1000. Figure 6.3 shows the dimensionless velocity and temperature for various Prandtl numbers obtained by using a shooting method with the Range-Kutta method. As the Prandtl number increases, the maximum velocity in the boundary layer decreases, and the location at which the peak velocity occurs shifts to smaller η. As shown in Fig. 6.3(a), the thickness of the momentum boundary layer decreases with increasing Prandtl number. Similarly, the thickness of the thermal boundary layer also decreases with increasing Prandtl number, as indicated by Fig. 6.3(b). The absolute value of the dimensionless temperature gradient also increases with increasing Prandtl number. After the dimensionless temperature in the boundary layer is obtained, the local heat transfer coefficient at the surface of the vertical plate can be obtained from eq. (6.3) as follows: hx = −  ∂T  1  Grx  k   = −kθ ′(0)  Tw − T∞  ∂y  y = 0 x 4   =− 1 1/ Grx / 4 = φ (Pr)Grx 4 1/ 4 (6.79) The local Nusselt number is hx x θ ′(0) Nu x = k 2 (6.80) Chapter 6 Natural Convection 529 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (a) Velocity profiles 1 Pr=(0.01,0.5,1,2,5,10,100,1000) 1 0 0 0 0 0 0 1 2 3 4 5 6 η=(Grx/4)1/4y/x θ=(T-T∞)/(Tw-T∞) 0.8 0.6 0.4 0.2 0 (b) Temperature profiles Figure 6.3 Velocity and temperature profile in the boundary layer for external natural convection over a vertical isothermal surface where φ (Pr) = −θ ′(0) / 2 is a function of Prandtl number. The dependence of φ on the Prandtl number is evidenced by eq. (6.74) and by Fig. 6.3(b). The values of φ for various Pr have been obtained numerically by Ostrach (1953). Ede (1964) proposed the following function that correlated the numerical results:  3 2 Pr 2 φ (Pr) =   1/ 2 4  5(1 + 2 Pr + 2 Pr)  1/ 4 The local Nusselt number thus becomes: 1/ Nu x = φ (Pr)Grx 4 =  3 2 Pr   4  5(1 + 2 Pr1/ 2 + 2 Pr)  1/ 4 1/ 4 (Grx Pr)1/4 (6.81) which can also be rewritten in terms of Rayleigh number Nu x =  3 2 Pr   1/ 2 4  5(1 + 2 Pr + 2 Pr)  Ra1/4 x (6.82) Equations (6.81) and (6.82) are valid for 0 < Pr < ∞ . As is demonstrated above, the similarity solution is obtained as a consequence of the geometrical similarity of the velocity and temperature profiles in the boundary layers, i.e., the velocity and temperature profiles vary with x according to the stretching functions, G(x) and H(x). If the geometrical 530 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1 1 Pr=(0.01,0.5,1,2,5, 10,100,1000) 0.8 0.6 0.4 0.2 0 0 1 2 η=Rax1/4 y/x 3 4 F'(η)=Rax -1/2 ux/α 0.8 0.6 0.4 0.2 0 (a) Velocity profiles 1 Pr=(0.01,0.5,1,2,5,10,100,1000) 1 0.8 0.6 0.4 0.2 0 0 1 2 η=Rax 1/4 y/x 3 4 0.8 0.6 0.4 0.2 0 θ=(T-T∞)/(Tw-T∞) (b) Temperature profiles Figure 6.4 Velocity and temperature profiles in the boundary layer based on modified scale. similarity exists, the selection of the stretching functions is not unique, and the resulting solutions based on different choices, as long as they all satisfy eqs. (6.65)–(6.67), are equivalent to the solution based on the stretching functions expressed in eqs. (6.70) and (6.71). However, a choice that better represents the physics of the problem will lead to the results presented in a physically more meaningful way. Figure 6.4 was obtained by modifying the numerical results shown in Fig. 6.3 by incorporating the modified stretching function for fluids of Pr > 1 (Bejan, 2004) expressed in eqs. (6.43) and (6.44), which is equivalent to a selection of the stretching functions in the following form G( x ) = 1 Ra1/ 4 x Pr 1 H ( x ) = Ra1/ 4 x x It is clear in Fig. 6.4 that, in the limit Pr→∞, the temperature profiles collapse onto a single curve, while the dimensionless velocity peak for fluids of Pr > 1 is consistently a number of order 1, showing that the velocity peak falls in the thermal boundary layer. Furthermore, as Pr increases, the velocity profile extends farther and farther into the isothermal fluid. All these features are Chapter 6 Natural Convection 531 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press anticipated by Fig. 6.1 and support the scale analysis, but cannot be seen from Fig. 6.3, in which the velocity and temperature profiles constantly shift as Pr changes, and the peak dimensionless velocity is not of order 1. Example 6.1 For the cases that Pr → 0 or Pr → ∞ , the similarity solution can be obtained by transforming eqs. (6.73) and (6.74) using the following new variables (Le Fevre, 1956; Kays et al., 2005): η1 = η ⋅ z1 (Pr) , F1 = F / z 2 (Pr) , and θ1 = θ , where z1 and z2 are unknown scaling coefficients. Obtain the similarity equation using the new similarity variables. Solution: To transform eqs. (6.73) and (6.74) to equations in terms of the new similarity variables, the relationship between the derivatives of original and new dimensionless stream function and temperature must be identified. The derivative of the original dimensionless temperature is: dθ dη1 = θ1′(η1 ) ⋅ z1 (Pr) θ ′(η ) = dη1 dη The second order derivative of original dimensionless temperature is: dη d d [θ ′(η )] 1 = [θ1′(η1 ) z1 (Pr)]z1 (Pr) = θ1′′(η1 ) z12 (Pr) θ ′′(η ) = dη1 dη dη1 Substituting the above derivatives and definition of F1 into eq. (6.74), the energy equation becomes θ1′′z12 + 3Pr F1θ1′z1 z 2 = 0 The choice of z1 and z2 should be such that the Prandtl number should not appear explicitly in the new energy equation, i.e., z 2 (Pr) = z1 (Pr) / (3Pr) . The derivatives of the stream function are dF z 2 (Pr) dF dη1 F ′(η ) = F1′(η1 ) = z 2 (Pr) 1 z1 (Pr) = 1 dη1 dη dη1 3Pr F ′′(η ) = F ′′′(η ) = dF ′ dη1 z13 (Pr) = F1′′ (η1 ) 3Pr dη1 dη dF ′′ dη1 z14 (Pr) F1′′′(η1 ) = dη1 dη 3Pr Substituting the above expressions into eq. (6.73), one obtains the following equation F1′′′(η1 ) 3F1 F1′′ (η1 ) − 2[ F1′(η1 )]2 θ1 3Pr + 9 Pr 2 + z14 (Pr) =0 where the three terms on the left-hand side represent the effects of viscosity, inertia, and buoyancy, respectively. The choice of the scaling factor must ensure that the equation is valid under limiting conditions. The results of the scaling analysis require that the viscous effect vanish 532 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press when Pr → 0 . On the other hand, the effect of inertia term should vanish when Pr → ∞ . Le Fevre (1956) suggested the following z14 (Pr) = 3Pr 2 1 + Pr Thus, the momentum and energy equations become: Pr( F1′′′+ θ1 ) + θ1 + F1 F1′′= 22 F1′ 3 (6.83) (6.84) θ1′′+ F1θ1′ = 0 Equations (6.83) and (6.84) allow simplification for the cases that Pr → 0 or Pr → ∞ (see Problem 6.4). Numerical solution of the simplified equation for these two cases yields the following results (Le Fevre, 1956): 0.600(Grx Pr 2 )1/ 4 = 0.600(Ra x Pr)1/ 4  Nu x =  1/ 4 1/ 4  0.503(Grx Pr) = 0.503Ra x  as Pr → 0 as Pr → ∞ (6.85) Equations (6.82) and (6.85) indicate that hx ∝ x −1/ 4 . The average Nusselt number over the entire vertical wall, Nu L , is related to the local Nusselt number at x = L by the following (see Problem 6.6): Nu L = 4 Nu L 3 (6.86) 6.4.2 Integral Solution for Laminar and Turbulent Natural Convection Laminar Flow Multiplying the continuity equation (6.17) by u and adding the resulting equation to the momentum equation (6.20) yields: ∂u 2 ∂ (uv) ∂2u + = ν 2 + gβ (T − T∞ ) ∂x ∂y ∂y Integrating the above equation with respect to y in the interval of (0, Y), where Y is greater than both δ and δ t , one obtains: dY2 ∂u 0 u dy = −ν ∂y dx y =0 + gβ  (T − T∞ )dy 0 Y (6.87) By following a similar procedure, the integral energy equation can be obtained as follows dY ∂T u(T − T∞ )dy = −α dx 0 ∂y (6.88) y =0 In order to obtain the solution of a natural convection problem, the appropriate velocity and temperature profiles must be utilized together with the integral equations (6.87) and (6.88). The velocity and temperature profiles depend on the thicknesses of the momentum and thermal boundary layers, which Chapter 6 Natural Convection 533 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press in turn depend on the Prandtl number. Assuming the velocity profile is the third degree polynomial function of y in the boundary layer and using the boundary conditions to determine the unspecified constant, the velocity profile becomes (see Problem 6.8): u y y = 1 −  U δ δ 2 (6.89) where U is a characteristic velocity that is a function of x. Similarly, the temperature profile can be obtained by assuming a second degree polynomial function and the result is (see Problem 6.9): T − T∞  y = 1 −  Tw − T∞  δ t  2 (6.90) The following analysis will be based on the assumption that the momentum and thermal boundary layers have the same thickness, i.e. δ t = δ . Substituting the velocity and temperature profiles into the integral form of the momentum equation (6.87) and energy equation (6.88) yields: 1d U1 (U 2δ ) = −ν + gβ (Tw − T∞ )δ 105 dx δ3 1d 2α (U δ ) = 30 dx δ (6.91) (6.92) At the leading edge of the vertical plate, the boundary layer thickness is zero and the characteristic velocity U is also zero: U = δ = 0, x = 0 (6.93) which are the initial conditions of eqs. (6.91) and (6.92). It is expected that as x increases, both U and δ should increase. Let us assume that they are functions of x such that U = C1 x m , δ = C 2 x n (6.94) Substituting eq. (6.94) into eqs. (6.91) and (6.92), one obtains the following: C 2m + n 2 1 C1 C 2 x 2 m + n −1 = − 1 ν x m − n + gβ (Tw − T∞ )C 2 x n C2 105 3 m+n 2α − n C1C 2 x m + n −1 = x C2 30 (6.95) (6.96) The above two relations can be true for all x only if the indices of x for all terms in the same equation are the same, i.e. when the following equations hold: 2m + n − 1 = m − n = n m + n − 1 = −n which are satisfied only if m = 1 / 2 and n = 1 / 4 . This suggests that U ∝ x1/ 2 and δ ∝ x1/ 4 ; this is in agreement with the result of the scaling analysis. Substituting back the values of m and n into eqs. (6.95) and (6.96), one obtains C12 C 2 C 1 = − 1 ν + gβ (Tw − T∞ )C 2 84 3 C2 534 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press C1C 2 2α = 40 C2 Solving for C1 and C2 from the above two equations yields: 1/ 2 −1/ 2 1/ 2 5  20 ν   gβ (Tw − T∞ )  C1 = 4   ν  +    ν2 3  21 α     15   20 ν  C2 = 4    +   16   21 α  = 3.93   1/ 4 1/ 4  gβ (Tw − T∞ )    ν2   Grx −1/ 4 ν   α  −1/ 2 The boundary layer thickness therefore becomes: 1/ 4 δ  0.952 + Pr  −1/ 4 x Pr 2   (6.97) The local heat transfer coefficient at the surface of the vertical plate can be obtained from eq. (6.3): hx = − k Tw − T∞  ∂T  2k  =  ∂y  y = 0 δ The local Nusselt number is: Nu x = 1/4 Substituting eq. (6.97) into the above expression yields:   Pr 2 Pr   1/ 4 1/ 4 Nu x = 0.508   Grx = 0.508   Ra x 0.952 + Pr  0.952 + Pr    1/4 hx x 2 x = k δ (6.98) Figure 6.5 shows the comparison between the integral solution, eq. (6.98), and the similarity solution, eq. (6.82). It can be seen that the integral solution under predicts the local Nusselt number for low Prandtl number but over predicts 0.6 0.5 0.4 Nux/Rax1/4 0.3 0.2 0.1 Integral Solution Similarity Solution 0.0 10-4 10-3 10-2 10-1 100 101 102 103 104 Pr Figure 6.5 Comparison between integral solution and similarity solutions for natural convection over a heated vertical wall. Chapter 6 Natural Convection 535 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the local Nusselt number for high Prandtl number. At Pr = 10−4 , the integral solution yields Nu x / Ra1x/ 4 = 0.051 , which is 13% lower than the value of 0.059 obtained from the similarity solution. At Pr = 0.72, which is the Prandtl number for air, the result obtained by the integral solution ( Nu x / Ra1/ 4 = 0.412 ) is 6.8% x higher than the similarity solution. It can also be seen that the agreement between the integral and similarity solutions is the best at high Prandtl number. When Pr = 104, the difference between the integral and similarity solution is only 1.5%. Example 6.2 The thickness of the momentum and thermal boundary layers are assumed to be the same in the above analysis. For the fluid with Prandtl number greater than 1, the following temperature and velocity profiles can be used (Bejan, 2004): T − T∞ = e − y /δt Tw − T∞ u = e − y / δ 1 − e − y / δt U (6.99) ( ) (6.100) where δ , δ t and U are unknown functions of x. Find the Nusselt number using integral solution based on the above profiles. Solution: Substituting eqs. (6.99) and (6.100) into eqs. (6.87) and (6.88) and setting Y → ∞ , the following equations are obtained:  d  U 2δ q 2 Uq δ + gβ (Tw − T∞ )   = −ν dx  2(2 + q)(1 + q)  q δ α d Uδ  = + q)(1 + 2q)  δ dx  (1 where δ δT is the ratio between the thicknesses of momentum and thermal boundary layers. Since q is a function of Prandtl number only, the above differential equations can be rewritten as: q2 d Uq δ (U 2δ ) = −ν + gβ (Tw − T∞ ) (6.101) q 2(2 + q)(1 + q) dx δ 1 d α (6.102) (U δ ) = (1 + q)(1 + 2q) dx δ An additional equation is needed since there are three unknowns in eqs. (6.101) and (6.102). Due to the non-slip and non-permeable boundary condition at the vertical wall (u = v = 0 at y = 0), eq. (6.20) in the region near the wall becomes: q(Pr) = 0 =ν ∂ 2u ∂y 2 + gβ (Tw − T∞ ) y =0 536 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Substituting eqs. (6.99) and (6.100) into the above momentum equation, one obtains: ν Uq(2 + q) (6.103) = gβ (Tw − T∞ ) δ2 Equations (6.101) and (6.102) are similar to eqs. (6.91) and (6.92); thus, one can expect that U ∝ x1/ 2 and δ ∝ x1/ 4 . By following the similar procedure that we used to solve eqs. (6.91) and (6.92), the local Nusselt number can be obtained: 3  q3 1/ 4 Nu x =   Ra x  8 (q + 1)(q + 1/ 2)(q + 2)  1/4 (6.104) where the ratio of the boundary layer thickness can be obtained from the following equation: 5 q +1/ 2 Pr = q 2 6 q+2 (6.105) The derivations of eqs. (6.104) and (6.105) are left as a homework problem (Problem 6.10). The above discussion is limited to the case where the surface temperature of the vertical plate is uniform. For the case where the surface heat flux of the flat plate is uniform, an integral solution can be performed to obtain the following results (see Problem 6.11): Nu x = 2  Pr  1/5   Ra *x 3601/5  0.8 + Pr  Ra *x = gβ q′′x 4 αν k 1/5 (6.106) where (6.107) is the modified Rayleigh number based on the heat flux. Turbulent Flow For small temperature difference between the heated wall and bulk fluid and short vertical plate, the natural convection is laminar and the above similarity and integral solutions are valid. Once the Rayleigh number exceeds a critical value, the natural convection will become turbulent and the above results will be invalid. It was believed that the transition from laminar to turbulent occurs at Ra x  109 until Bejan and Lage (1990) showed that for a wide range of Prandtl number ( 0.001 < Pr < 1000 ) the criterion for transition from laminar to turbulent is actually Grx  109 . Alternatively, one can say that the transition takes place at Ra x  109 Pr . Thus, the critical Rayleigh number for low-Prandtl number fluid is less than 109. For high-Prandtl number fluid, on the other hand, the critical Rayleigh number is higher than 109. The advantage of the integral solution is that it also works for turbulent flow. The integral momentum and energy equations (6.87) and (6.88) are applicable Chapter 6 Natural Convection 537 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press x Turbulent Tw quiescent fluid, u∞=0 Laminar gravity y Figure 6.6 Transition from laminar to turbulent if all terms are time-averaged as should be used for turbulent flows. Using the definition of shear stress and heat flux, eqs. (6.87) and (6.88) can be respectively modified as Y τw dY2 (6.108) 0 u dy = gβ 0 (T − T∞ )dy − ρ dx ′′ qw dY 0 u(T − T∞ )dy = ρ c p dx (6.109) The velocity and temperature profiles in the boundary layer must be correctly determined so as to reflect the behavior of the turbulent boundary layer. A velocity profile constructed as follows gives a good description of the velocity distribution in natural convection for turbulent flow over a vertical flat plate (Eckert and Jackson, 1951): u y =  U δ  1/ 7  y 1 −   δ 4 (6.110) where U is a characteristic velocity for the near wall. Equation (6.110) satisfies all of the velocity boundary conditions. It is further assumed that the velocity and temperature boundary layer thicknesses are the same ( δ t = δ ) and the temperature profile is: T − T∞ y = 1−   Tw − T∞ δ  1/ 7 (6.111) which yields T = Tw at y = 0 and T = T∞ at y = δ . It should be pointed out that eqs. (6.110) and (6.111) are valid only for y < δ . For large y, one has u = 0 and T = T∞ . 538 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Substituting eqs. (6.110) and (6.111) into eqs. (6.108) and (6.109), the integral equations become: 1 τ d  2 1 2/ 7 8 U δ  η (1 − η ) dη  = gβ (Tw − T∞ )δ  (1 − η 1/ 7 )dη − w   0 0   dx ρ 1 q′′ d 4 U (Tw − T∞ )δ  η 1/ 7 (1 − η ) (1 − η 1/ 7 ) dη  = w   ρc 0  dx  p where η = y / δ . Evaluating the integrals in the above two equations yields: τ d (6.112) 0.0523 (U 2δ ) = 0.125gβ (Tw − T∞ )δ − w dx ρ 0.0366(Tw − T∞ ) q′′ d (U δ ) = w ρcp dx (6.113) ′′ At this point, we have two equations and four unknowns (U, δ , qw , τ w ). Thus, assumptions have to be made regarding the forms of the expressions for ′′ τ w and qw . It is generally assumed that the flow near the wall in a turbulent natural convective boundary layer is similar to that of a turbulent forced ′′ convection so that the expressions for τ w and qw derived for forced convection can be applied: τw 0.0225 = (6.114) ρU 2 ( ρU δ / μ )0.25 ′′ qw ν  = 0.0225   ρ c pU (Tw − T∞ )  Uδ  1/ 4 Pr −2/3 (6.115) Substituting the above expressions into the integral momentum and energy equations, one obtains: 0.0523 1.75 d (U 2δ ) = 0.125gβ (Tw − T∞ )δ − 0.0225ν 0.25 U0.25 dx δ d U 0.75 0.0366 (U δ ) = 0.0225 Pr −0.67 ν 0.25 0.25 dx δ (6.116) (6.117) The boundary layer can be assumed to be turbulent from the leading edge of the surface, so the solutions to the above equations should be of the following form: U = C1 x m , δ = C 2 x n (6.118) Substituting the above equations into eqs. (6.116) and (6.117), one finds that the values which satisfy these two equations are m = 0.5 and n = 0.7 (Problem 6.13). By following the procedure similar to the case of laminar natural convection (Problem 6.14), the local Nusselt number can then be obtained: Nu x = 0.0295 Pr1/15 Ra 2/5 x (1 + 0.494 Pr 2/3 ) 2/5 (6.119) which suggests that the local heat transfer coefficient is proportional to x0.2. The average Nusselt number for the entire vertical plate, Nu L , is related to the local Nusselt number at x = L as follows (see Problem 6.15). Chapter 6 Natural Convection 539 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Nu L = 0.834Nu L (6.120) Empirical Correlation The above analyses for laminar and turbulent natural convection over a vertical flat plate suggest that the average Nusselt number can be expressed in the following format: Nu L = C Ra n (6.121) L which is confirmed by experimental studies. For laminar flow ( GrL ≤ 109 ), one can utilize the values of C = 0.59 and n = 1 / 4 . For turbulent flow ( GrL > 109 ), C = 0.1 and n = 1 / 3 can be used (McAdams, 1954; Warner and Arpaci, 1968). Churchill and Chu (1975) studied numerous sets of experimental data and recommended the following correlation:   0.387Ra1/ 6 L Nu L = 0.825 + 9/16 8/ 27  [1 + (0.492 / Pr) ]   2 (6.122) which covers all Prandtl number and Grashof number between 0.1 and 1012. For the case of laminar convection ( Grx < 109 ), the following correlation yields better results: Nu L = 0.68 + 0.67Ra1/ 4 L [1 + (0.492 / Pr)9/16 ]4/9 (6.123) Practically, one can use eq. (6.123) for laminar natural convection and eq. (6.122) for turbulent natural convection. These two correlations can provide better accuracy than the earlier correlation in simpler form. For the case of constant heat flux where the surface temperature, Tw, increases with increasing x, Churchill and Chu (1975) suggested that eq. (6.122) is still valid provided the constant 0.492 is changed to 0.437, and Nu L and RaL must be defined using the averaged wall temperature of the vertical plate. 6.4.3 Natural Convection over Inclined and Horizontal Surfaces Inclined Flat Surface Depending on whether the wall is heated or cooled and tilted upward or downward, there are four different configurations for convection over an inclined flat plate (see Fig. 6.7). While the effect of inclination thickens the boundary layers for upper surface of case (a) and bottom surface of case (b), the boundary layer becomes thinner for bottom surface of case (a) and upper surface of case (b). When the inclination angle is moderate ( −60 < φ < 60 ), the momentum equation for laminar natural convection over an inclined wall is identical to eq. (6.20), except that the gravitational acceleration g must be replaced by the component of the gravitational acceleration along the direction that is parallel to the inclined plate, g cos φ . For laminar natural convection, the empirical 540 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press T w < T∞ g g Tw > T∞ (a) heated surface Figure 6.7 Natural convection over inclined surface (b) cooled surface correlations that were developed in the previous subsection, e.g. eq. (6.123), are still valid if the Rayleigh number is defined as: g cos φβ (Tw − T∞ ) x 3 Ra x = (6.124) αν For turbulent flow, Vliet (1969) suggested that the experimental results correlated better using g, instead of g cos φ , i.e. heat transfer in turbulent natural convection is not sensitive to the inclination angle and eq. (6.122) can be used. Horizontal Flat Surface Natural convection on a horizontal surface can find its application in nature such as natural convection on ground or water surfaces, as well as on the top of a building. It can also be applied to electronic cooling. Due to its complicated nature, analytical solutions are only available for limited cases. Rotem and Claasen (1969) studied natural convection over a semi-infinite isothermal flat plate with a single leading edge (see Fig. 6.8). The fluid near the flat plate is Quiescent fluid, u∞=0 y gravity T∞, ρ∞ Velocity boundary layer, δ T u Tw Thermal boundary layer, δt x Figure 6.8 Natural convection over a horizontal flat plate Chapter 6 Natural Convection 541 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press heated and becomes lighter, which results in a pressure difference that causes flow over the surface near the leading edge. The continuity and energy equations are the same as the case of the vertical plate, eqs. (6.17) and (6.21). The momentum equation is, however, different from that of the vertical plate because the direction of gravity is normal to the horizontal plate. The momentum equations are: u ∂u ∂u ∂ 2 u 1 ∂pd +v =ν 2 − ∂x ∂y ρ ∂x ∂y gβ (T − T∞ ) = 1 ∂pd ρ ∂y (6.125) (6.126) where pd is the dynamic pressure whose gradient is the driving force for the boundary layer flow. Defining the following similarity variables: 1/5 1/5 4/5 5ν 2 ρ y  Gr   Gr   Gr  F (η )  x  G (η ) η =  x  , ψ = 5ν F (η )  x  , pd = (6.127) x 5  5 x 5 eqs. (6.17), (6.125), (6.126) and (6.21) become (Pera and Gebhart, 1973): 2 2 F ′′′ + 3F ′′F − 2 f ′2 + G ′η − G = 0 5 5 G′ = θ θ ′′ + 3 Pr F θ ′ = 0 (6.128) (6.129) (6.130) which can be numerically solved to obtain the following local Nusselt number: 1 Nu x = 0.394Grx /5 Pr1/ 4 (6.131) Equation (6.131) is valid for an isothermal surface. For a horizontal surface with constant heat flux, the local Nusselt number is: 1 Nu x = 0.501Grx /5 Pr1/ 4 (6.132) While the above solution is for the case of a heated surface facing upward, it can also be applied to the case of a cold surface facing downward. However, the similarity solution cannot be applied for the case of a hot surface facing downward or a cold surface facing upward. For most applications, the horizontal surface cannot be treated as a semiinfinite surface with only one leading edge. Two extreme types of flow situation arise for the horizontal surface with finite size. If the hot surface of the plate faces upwards or cold surface of the plate faces downwards, a plume is formed whereby the hotter fluid layers rise up or colder fluid mass goes down, respectively. For the cases that heated surface facing downward [bottom surface of Fig. 6.9 (a)] or cold surface facing upward [upper surface of Fig. 6.9. (b)], the flow leaves the boundary layer as a vertical plume at the central region. If the temperature difference happens to be large enough, then the flow occurs from all over the horizontal surface. On the contrary, when the hot surface faces downward [bottom surface of Fig. 6.9(a)] or the cold one faces upward [upper surface of Fig. 6.9(b)], the boundary layer covers the entire surface and the flow is spilled over the edges. 542 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Tw < T∞ Tw > T∞ (a) hot surface Figure 6.9 Natural convection over horizontal surface (b) cold surface For the case of upward facing hot isothermal surfaces or downward facing cold isothermal surfaces, the following correlation can be applied (McAdams, 1954): 0.54 Ra1/ 4 104 < Ra L < 107  L Nu L =  1/3 7 9  0.15 RaL 10 < Ra L < 10  (6.133) The characteristic length in the average Nusselt number and Rayleigh number is L = A / P , where A refers to the area of the plate and P is the perimeter. It should be pointed out that the average Nusselt numbers in eqs. (6.133) and 1 (6.134) are proportional to Ra1L/ 4 or Ra1/3 , not GrL/5 Pr1/4 as was suggested by the L above similarity solution. Thus, the applicability of the similarity solution is limited to natural convection over a horizontal surface due to the complicated nature of the problem. The corresponding correlation for downward facing hot surfaces or upward facing cold surfaces is: Nu L = 0.27 Ra1/ 4 105 < Ra L < 1010 (6.134) L While eqs. (6.133) and (6.134) are for isothermal surfaces, they can also be used for the case of uniform heat flux provided the average Nusselt and Rayleigh numbers are defined based on the average surface temperature of the flat surface. 6.4.4 Natural Convection over Cylinders and Spheres Vertical Cylinder Natural convection on a vertical cylinder can find its applications in the cooling of nuclear reactors, large-current bus in power plants, as well as in the human body. A vertical cylinder has three surfaces: the top surface, the bottom surface and the vertical surface. The natural convection from the top and bottom surfaces can be calculated using the approaches discussed in the previous subsection. For the vertical surface of the cylinder, the development of the boundary layer is shown in Fig. 6.10. If the ratio of diameter to height, D/L, is Chapter 6 Natural Convection 543 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δt L Tw x, u T∞ r, v D = 2R Figure 6.10 Boundary layer for natural convection over vertical cylinder large enough, the thermal boundary layer thickness will be much smaller than the radius of the cylinder and the correlations for natural convection over a vertical flat plate can be used to calculate the natural convection over the vertical cylinder. For fluid with Prandtl number greater than 1, the condition under which the correlation for flat plate is applicable is (Bejan, 2004): D > Ra −1/ 4 L L (6.135) This condition is referred to as the “thick cylinder limit.” If D/L is less than Ra −1/ 4 , the boundary layer thickness will be comparable to the radius of the L cylinder and the effect of surface curvature cannot be neglected. For this case – referred to as “thin cylinder limit” – the governing equations must be given in a cylindrical coordinate system: ∂u 1 ∂ (vr ) + =0 ∂x r ∂y (6.136) (6.137) (6.138) u ∂u v ∂ (ur ) 1 ∂  ∂u  + =ν  r  + gβ (T − T∞ ) r ∂r  ∂r  ∂x r ∂r ∂T v ∂ (Tr ) 1 ∂  ∂T  + =α u r  ∂x r ∂y r ∂r  ∂r  The boundary conditions for eqs. (6.136) – (6.138) are (6.139) u = v = 0, T = T∞ , at r → ∞ (6.140) The integral momentum and energy equations can be obtained by integrating eqs. (6.136) – (6.138) with respect to r in the interval of (R, R + δ): u = v = 0, T = Tw , at r = R R +δ d R +δ 2 ∂u R u rdr = −ν R ∂r r = R + gβ R (T − T∞ )rdr dx (6.141) (6.142) d R +δ ∂T R u(T − T∞ )rdr = −α R ∂r dx r=R 544 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press In order to simplify the solution procedure, it is assumed that the velocity and temperature profiles that were used in the integral solution of natural convection over a vertical flat plate, eqs. (6.89) and (6.90), are still valid for natural convection over a thin cylinder – except to replace y with r − R. Substituting eqs. (6.89) and (6.90) into eqs. (6.141) and (6.142), one obtains the following integral equations: d U (8R + 3δ )δ U 2  = −840ν R + 70 gβ (Tw − T∞ )δ (4 R + δ )  δ dx  d 420α R [(7 R + 2δ )U δ ] = dx δ (6.143) (6.144) Unlike the case of natural convection over a vertical flat plate, eqs. (6.143) and (6.144) do not suggest that U and δ can be expressed as a simple function of x similar to eq. (6.94). Le Fevre and Ede (1956) assumed that: UU δδ U = U 0 + 1 + 2 +  and δ = δ 0 + 1 + 22 +  (6.145) 2 R R R R Substituting the above expression into eqs. (6.143) and (6.144) and grouping the terms for the same order of R, equations for U0, δ0, U1, δ1, and so on, can be obtained: d 2 (U 0 δ 0 ) = −840ν U 0 + 280gβ (Tw − T∞ )δ 02 dx d 7δ 0 (U 0δ 0 ) = 420α dx d d 2 2 2 δ0 (8δ 02U 0 + 8δ1U 0 + 16δ 0U 0U1 ) + 8δ1 (δ 0U 0 ) dx dx = −840ν U1 + 70 gβ (Tw − T∞ )(8δ 0δ1 + δ 03 ) 8δ 0 (6.146) (6.147) (6.148) (6.149) δ0 d d (2δ 02U 0 + 7δ1U 0 + 7δ 0U1 ) + 7δ1 (δ 0U 0 ) = 0 dx dx It can be seen that eqs. (6.146) and (6.147) have a form similar to eqs. (6.91) and (6.92). Thus, one can assume that: U 0 = A0 x m , δ 0 = B0 x n and one can use eqs. (6.146) and (6.147) to obtain m = 1 / 2 and n = 1/ 4 and: 80α 80(20 + 21Pr)α 2 4 A0 = 2 , B0 = (6.150) B0 7 gβ (Tw − T∞ ) By following a similar procedure, one can use eqs. (6.148) and (6.149) to obtain: U1 = A1 x 3/ 4 , δ1 = B1 x1/ 2 where 2 (272 + 315 Pr) B0 4(656 + 315 Pr)α , B1 = − A1 = − (6.151) 7(64 + 63Pr) B0 35(64 + 63Pr) If one only considers the first two terms in eq. (6.145), the boundary layer thickness becomes: Chapter 6 Natural Convection 545 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ ≈ B0 x1/ 4 + hx = − B1 x1/ 2 R (6.152) The local heat transfer coefficient at the surface of the vertical cylinder is: k 2k  ∂T  =   Tw − T∞  ∂r r = R δ and the local Nusselt number is: Nu x = hx x 2 x = k δ Substituting eq. (6.152) into the above expression, one obtains  7Grx Pr  2   80(20 + 21Pr)  Nu x = 3/ 4  80(20 + 21Pr)  272 + 315 Pr x 1−   7Grx Pr 2  35(64 + 63Pr) R  1/ 4 Since the second term in the denominator of the above expression is much less than 1, Le Fevre and Ede (1956) suggested that the above expression for local Nusselt number can be simplified to:  7Grx Pr 2  Nu x =    5(20 + 21Pr)  Nu L = 1/ 4 + 2(272 + 315 Pr) x 35(64 + 63Pr) R 4(272 + 315 Pr) L 35(64 + 63Pr) D (6.153) The average Nusselt number can then be obtained as: 4  7Ra L Pr    3  5(20+21Pr)  1/ 4 + (6.154) where the characteristic length in Nu L and RaL is the height of the cylinder. It can be seen that the first term on the right-hand side of eq. (6.154) has a form similar to eq. (6.98), and that the second term represents the effect of aspect ratio of the vertical cylinder. The second term does not depend on the Rayleigh number. Example 6.3 A cone with the surface temperature of Tw is immersed into a quiescent fluid at a temperature of T∞ ( Tw > T∞ ). A boundary layer develops from the tip of the cone and thickens as x increases (see Fig. 6.11). Assuming that the thicknesses of the momentum and thermal boundary layers are the same, obtain the average Nusselt number. Solution: The dimensionless continuity, momentum and energy equations are (Alamgir, 1979): ∂ (UR) ∂ (VR) + =0 ∂X ∂Y ∂U ∂U 1 ∂  ∂U  +V = U R  + GrLθ ∂X ∂Y R ∂Y  ∂Y  (6.155) (6.156) 546 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ Tw L x,u φ T∞ 0 y, v Figure 6.11 Boundary layer for natural convection over vertical cone u ∂θ ∂θ Pr ∂  ∂θ  +V = R  ∂x ∂y R ∂Y  ∂r  (6.157) where the dimensionless variables are defined as: T − T∞ x y r uL vL δ X = ,Y = ,R= ,Δ= ,U = ,V = ,θ = ν ν L L L L Tw − T∞ and R is related to X and Y by: R = X sin φ + Y cos φ The boundary conditions of eqs. (6.155) – (6.157) are: U = V = 0, θ = 1, at Y = 0 (6.158) U = 0, θ = 0, at Y → Δ (6.159) The integral momentum and energy equations can be obtained by integrating eqs. (6.155) – (6.157) with respect to Y in the interval of (0, Δ): A 2  BΔ  U1 X Δ  1 + A X cot φ       DΔ  X ∂U = C GrL X Δ 1 + cot φ  −  CX  Δ ∂Y d dX (6.160) Y =0 E Pr d X ∂θ  FΔ  cot φ   = − U1 X Δ 1 + dX  EX Δ ∂Y    (6.161) Y =0 where U1 is a dimensionless characteristic velocity. The other variables are defined as 1 0 A =  (U / U1 ) 2 dη , 0 1 1 0 0 1 B =  η (U / U1 ) 2 dη , 0 1 C =  θ dη , 0 1 D =  ηθ dη , E =  (U / U1 )θ dη , F =  η (U / U1 )θ dη where η =Y /Δ . The effect of curvature is reflected by terms that contain (Δ / X ) cot φ . For the cases that the effect of curvature is negligible, eqs. (6.160) and (6.161) respectively become Chapter 6 Natural Convection 547 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press A d X ∂U (U12 X Δ ) = C GrL X Δ − dX Δ ∂Y E Pr d X ∂θ (U1 X δ ) = − dX Δ ∂Y Y =0 (6.162) Y =0 (6.163) Assuming the velocity and temperature in the boundary layer are described by the following equations: 3 U = U1η (1 − η ) θ = (1 − η ) 2 which satisfy eqs. (6.158) and (6.159), eqs. (6.162) and (6.163) become: A UX d (U12 X Δ) = 84GrL X Δ − 252 1 dX Δ d 84 X (U 1 X δ ) = dX Pr Δ (6.164) (6.165) Since eqs. (6.164) and (6.165) have a form similar to eqs. (6.91) and (6.92), one can use a similar approach to obtain: X Δ = a0    GrL  336 1 U1 = 2 GrL/ 2 X 1/ 2 7a0 Pr 1/ 4 (6.166) (6.167) where 1008(3 + 7 Pr)  a0 =   49 Pr 2   1/ 4 (6.168) Equations (6.166) and (6.167) are valid for the case that the effect of curvature is negligible. To consider the effect of curvature, Alamgir (1979) suggested replacing (Δ / X ) cot φ in eqs. (6.160) and (6.161) by the following: − Δ / X cot φ ≈ (Δ / X ) cot φ = 4a0 GrL 1/ 4 cot φ In this case, eqs. (6.160) and (6.161) become: (1 + 0.6a0ε ) UX d (U12 X Δ) = 84(1 + 0.5a0ε ) X Δ − 252 1 dX Δ d X Pr(1 + 0.5a0ε ) (U1 X Δ) = 84 dX Δ (6.169) (6.170) where ε = 2GrL−1/ 4 cot φ (6.171) is a curvature parameter. Equations (6.169) and (6.170) can be solved to obtain: X Δ = a1    GrL  1/ 4 (6.172) 548 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U1 = 336 1 GrL/2 X1/2 7a12 Pr(1 + 0.5a0ε ) 1/ 4 (6.173) where  3(1 + 0.6a0ε ) + 7(1 + 0.5a0ε ) Pr  a1 = 1008  (1 + 0.5a0ε )[7(1 + 0.5a0ε ) Pr]2   (6.174) The local wall heat flux is: q′′ = −k ∂T ∂y = −k y=0 Tw − T∞ ∂θ T − T∞ = 2k w L Δ ∂η η = 0 LΔ 1 At (6.175) The average heat flux is: q′′ =  A q′′dA (6.176) and where is the total surface of the cone, dA = 2π x sin φ dx . The average Nusselt number can be obtained as: Nu L = hL q′′L = k k (Tw − T∞ ) At = π L2 sin φ Substituting eqs. (6.175) and (6.176) into the above expressions, the final result becomes: Nu L = 1 4GrL / 4 a1  1 0 X 3/ 4 dX = 1/ 16GrL 4 7a1 (6.177) Horizontal Cylinder and Sphere When a horizontal cylinder or sphere with temperature Tw is immersed in a fluid with temperature T∞ ( Tw > T∞ ), a boundary layer develops along the curved surface. The boundary layer thickness is a function of the angle φ ( φ = 0 is at the bottom of the cylinder or sphere), as shown in Fig. 6.12. The similarity solution that worked for the case of the vertical plate does not work for natural convection over a horizontal cylinder or sphere. Merk and Prins (1953, 1954a, b) obtained an integral solution by assuming the momentum and thermal boundary layer thicknesses are identical. The local Nusselt numbers along the surfaces of the cylinder and sphere obtained by Merk and Prins (1954b) are shown in Fig. δt Tw D T∞ φ r Figure 6.12 Boundary layer for natural convection over horizontal cylinder or sphere Chapter 6 Natural Convection 549 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (a) horizontal cylinder (b) Sphere Figure 6.13 Local Nusselt number for natural convection over horizontal cylinder or sphere (Merk and Prins, 1954b) 6.13. It can be seen that the local Nusselt number is highest at the bottom where the boundary layer is thinnest. As the angle φ increases, the thickness of the boundary layer increases and the local Nusselt number decreases. Although the integral solution can yield results all the way to the top where φ = 180 and Nuφ = 0 , the result beyond φ = 165 is no longer applicable because boundary layer separation occurred and plume flow takes place. Based on the integral solution, Merk and Prins (1954b) recommended the following empirical correlation for natural convection over a horizontal cylinder and sphere: Nu D = C Ra D1/ 4 (6.178) where the characteristic length in the average Nusselt number and Rayleigh number is the diameter of the cylinder or sphere. The constant C is a function of Prandtl number and can be found from Table 6.1. Table 6.1 Constant in the empirical correlation for natural convection over horizontal cylinder and sphere Prandtl number Cylinder C Sphere 0.7 0.436 0.474 1.0 0.456 0.497 10.0 0.520 0.576 100.0 0.523 0.592 ∞ 0.523 0.595 550 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Practically, the empirical correlations based on experimental results are more useful. Churchill and Chu (1975) recommended the following correlation for horizontal cylinders:   0.387Ra1/ 6 D Nu D = 0.6 + 9/16 8/ 27  [1 + (0.559 / Pr) ]   2 (6.179) which has the same form as the correlation for vertical plate, eq. (6.122), except the characteristic length has been changed from vertical plate height to the diameter of the cylinder. Equation (6.179) covers all Prandtl number and Rayleigh number between 0.1 and 1012. For natural convection over a sphere, Churchill (1983) suggested that the following correlation can best fit the available experimental data: Nu D = 2 + 0.589Ra1/ 4 D [1 + (0.469 / Pr)9/16 ]4/9 (6.180) where the average Nusselt number, Nu D , and the Rayleigh number, RaD, are based on the diameter of the sphere. Equation (6.180) is valid for the case where Pr ≥ 0.7 and the Rayleigh number is less than 1011. 6.4.5 Free Boundary Flow The common feature of external convection discussed in the preceding subsections is that the flows are always near a heated (or cooled) solid wall. For some other applications, however, the thermally induced flow occurs without presence of a solid wall and is referred to as free boundary flow. When a point or line heat source is immersed into a bulk fluid, the fluid near the heat source is heated, becomes lighter and rises to form a plume [see Fig. 6.14 (a)]. The plume generated by a point heat source is axisymmetric, but a line heat source will create a two-dimensional plume. A thermal is a column of rising air near the (a) plume Figure 6.14 Examples of free boundary flows (b) thermal Chapter 6 Natural Convection 551 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press x, u T0 T∞ g y, v Heat source Figure 6.15 Physical model of natural convection from a line heat source ground due to the uneven solar heating of the ground surface. The lighter air near the ground rises and cools due to expansion [see Fig. 6.14(b)]. Let us consider now the two-dimensional free boundary flow induced by a line heat source (see Fig. 6.15). Since the problem is symmetric about x = 0, only half of the domain ( x > 0 ) needs to be studied. In the coordinate system shown in Fig. 6.15, the continuity, momentum and energy equations are the same as those for natural convection near a vertical flat plate, eqs. (6.17), (6.20) and (6.21), respectively. The boundary conditions at y → ∞ can still be described by eqs. (6.24) and (6.25). However, the boundary conditions at the centerline, x = 0 , should be changed to: ∂u ∂T = = 0, v = 0 at y = 0 ∂y ∂y (6.181) which indicates that the velocity component in the x-direction and the temperature at the centerline are at their respective maximum. Gebhart et al. (1970) assumed that the difference between the centerline temperature and the bulk fluid temperature is of the following power law form: T0 − T∞ = Nx n (6.182) where the index n will be determined using the overall energy balance. The continuity, momentum, and energy equations can be transformed to the following ordinary differential equations via the same similarity variables for natural convection over a vertical flat plate in Section 6.4.1 and considering eq. (6.182). F ′′′ + θ + (3 + n ) FF ′′ − (2n + 2)( F ′) 2 = 0 (6.183) θ ′′ + Pr[(n + 3) F θ ′ − 4nF ′θ ] = 0 (6.184) Equations (6.183) and (6.184) are applicable to any natural convection problem that satisfies eq. (6.182). To get the ordinary differential equations that are specific for natural convection induced by a line heat source, proper values of N and n, as well as appropriate boundary conditions, must be specified. The energy convected across any horizontal plane in the plume is 552 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press q′ = ρ c p  u(T − T∞ )dy −∞ ∞ (6.185) which is identical to the intensity of the line heat source and it can be expressed in terms of similarity variables 1/ 4 ∞  gβ N  (5 n + 3)/ 4 (6.186) q′ = 4νρ c p N  x 2 −∞ F ′θ dη  4ν  Since the line heat source is the only source of heating, the above q′ should be independent from x which can be true only if n is equal to −3 / 5 . Therefore, the ordinary differential equations for natural convection over a line heat source respectively become: F ′′′ + θ + 12 4 FF ′′ − ( F ′) 2 = 0 5 5 12 θ ′′ + Pr( F θ ′ + F ′θ ) = 0 5 (6.187) (6.188) which are subject to the following boundary conditions: F (η ) = F ′′(η ) = 0, and θ (η ) = 1 at η = 0 (6.189) F ′(η ) = θ (η ) = 0 at η → ∞ (6.190) Equations (6.187) – (6.190) can be solved numerically and the results are shown in Fig. 6.16. In contrast to Fig. 6.3 that shows the velocity at the wall is zero, the velocity at the centerline is at its maximum. With known intensity of the line heat source, the constant N in eq. (6.182) can be obtained from eq. (6.186):  q ′4 N =  64 gβρ 2 μ 2 c 4 I 4 p      1/5 (6.191) Figure 6.16 Velocity and temperature distributions of natural convection from a line heat source (a) velocity, (b) temperature (Gebhart et al., 1970; reproduced with permission from Elsevier) Chapter 6 Natural Convection 553 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where I =  F ′(η )θ (η )dη −∞ ∞ (6.192) is a function of Prandtl number. The values of I calculated at Pr = 0.7, 1.0, 6.7, and 10.0 by Gebhart et al. (1970) are 1.245, 1.053, 0.407 and 0.328, respectively. Example 6.4 Analyze natural convection over a point heat source to obtain the appropriate ordinary differential equations and the corresponding boundary conditions. Solution: The natural convection in this case will be axisymmetric and the governing equations must be given in a cylindrical coordinate system (Gebhart et al., 1988). Equations (6.136) – (6.138) are applicable for this case with the following boundary conditions ∂u ∂T = = 0, v = 0 at r = 0 ∂r ∂r u = v = 0, T = T∞ at r → ∞ (6.193) (6.194) Similar to the case of line heat source, the centerline temperature can also be expressed as a power function of x as indicated by eq. (6.182). Introducing Stokes stream function: 1 ∂ψ 1 ∂ψ , v=− u= (6.195) r ∂r r ∂x and defining the following similarity variables: η= T − T∞ r 1/ 4 Grx , ψ = ν xF (η ), θ = x T0 − T∞ (6.196) the following ordinary differential equations are obtained:  F ′ ′ 1 + n ( F ′) 2 =0 F ′′′ + ηθ + ( F − 1)   − η 2 η  (ηθ ′)′ + Pr( F θ ′ − nF ′θ ) = 0 (6.197) (6.198) ∞ The intensity of the point heat source is: q = ρ c p  u(T − T∞ )2π rdr = 2πνρ c p Nx n +1  F ′θ dη 0 0 ∞ (6.199) and it must be independent from x. Therefore, the value of n must be −1 and the ordinary differential equations for natural convection over a point heat source become:  F ′ ′ F ′′′ + ηθ + ( F − 1)   = 0 η  (ηθ ′)′ + Pr( F θ )′ = 0 (6.200) (6.201) The velocity components in the x- and r- directions are related to the dimensionless variables by the following equations:  F′ ν ν 1  F F′ u= Grx   , v = − Grx / 4  −  (6.202) x x η  η 2  554 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press which can be substituted into eqs. (6.193) and (6.194) to obtain the following boundary conditions for eqs. (6.200) and (6.201):  F ′ ′   = 0, F = 0, θ ′ = 0 and θ = 1 at η = 0 η  F′ = θ = 0 at η → ∞ (6.203) (6.204) η 6.5 Natural Convection in Enclosures Enclosures are referred to as finite spaces bounded by walls and filled with fluid (liquid or gas). Natural convection in enclosures, also known as internal convection, takes place in rooms and buildings, furnaces, cooling towers, fire safety, as well as electronic cooling systems. Internal natural convection is different from the cases of external convection, where a heated or cooled wall is in contact with the quiescent fluid and the boundary layer can be developed without any restriction. As will become evident later, internal convection usually cannot be treated using simple boundary layer theory because the entire fluid in the enclosure engages to the convection. Figure 6.17 illustrates several examples of both two- and three dimensional internal natural convections (Yang, 1987). Figures 6.17(a) and (b) depicts natural convection in a rectangular enclosure with left and right sides heated and cooled while the top and bottom walls are insulated. For a shallow enclosure ( H / L < 1 ), a boundary layer may exist near the heated and cooled wall and the buoyancy flow in the enclosure is of the recirculating type. For the case that H/L is greater than unity, natural convection in the enclosure can no longer be described using the boundary layer theory. Natural convection in an inclined rectangular enclosure can be considered as a generalized case for natural convection in a rectangular enclosure since Fig. 6.17(a) and (b) are special cases with γ = 90 . When γ = 180 , Fig. 6.17(c) becomes the case of natural convection in a rectangular enclosure heated from below. Natural convection in an enclosure can be suppressed by vertical partitions such as that shown in Fig. 6.17(d). Natural convection in an annular enclosure formed by two differentially heated concentric cylinders is illustrated in Fig. 6.17(e). The flow will be of the recirculating type and axisymmetric. Figure 6.17 (f) and (g) depict natural convection in a box enclosure and truncated annular enclosure, respectively. Both of them can be tilted relative to gravity and various heating and cooling boundary conditions can be imposed at different walls. It should be pointed out that the examples shown in Fig. 6.17 are illustrative, rather than inclusive. Natural convection can occur in different and more complex enclosures. In this section, scale analysis, analytical solutions, numerical solutions and empirical correlations pertaining to natural convection in enclosures will be presented. Chapter 6 Natural Convection 555 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press L L TH H TC TH H TC TH γ L H TC (a) shallow enclosure (b) tall enclosure (c) inclined enclosure Do Di L TH H TC (d) enclosure with vertical partitions a) ( (e) concentric annulus Do Di L TH γ TC D H (f) box enclosure (g) truncated annular enclosure Figure 6.17 Different configuration of natural convection in enclosures 6.5.1 Scale Analysis Two-dimensional natural convection in a rectangular enclosure with two differentially heated sides and insulated top and bottom surfaces (Fig. 6.18) will be considered. The fluid is assumed to be Newtonian and incompressible. The system is initially at a uniform temperature of zero and therefore, the fluid is motionless ( u = v = 0 ). At time zero the two sides are instantaneously heated and cooled to ΔT / 2 and −ΔT / 2 , respectively. The transient behavior of the system during the establishment of the natural convection (Patterson and Imberger, 1980; Bejan, 2004) is the subject of this analysis. It – is assumed that the fluid is single-component and that there is no internal heat generation in the fluid. Therefore, the governing equations for this internal convection problem can be obtained by simplifying eqs. (6.8), (6.13) and (6.14): 556 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press y, υ H g T=+ΔT/2 Thermal boundary layer Flow circulation direction T=-=–ΔT T ΔT/2 δT 0 0 L Figure 6.18 Two-dimensional natural convection in rectangular enclosure. x, u ∂u ∂v + =0 ∂x ∂y (6.205) (6.206) (6.207) (6.208)  ∂ 2u ∂2u  ∂u ∂u ∂u 1 ∂p +u +v =− +ν  2 + 2  ∂t ∂x ∂y ∂y  ρ ∂x  ∂x  ∂2v ∂2v  ∂v ∂v ∂v 1 ∂p +u +v =− + ν  2 + 2  − g[1 − β (T − T0 )] ∂t ∂x ∂y ∂y  ρ ∂y  ∂x  ∂ 2T ∂ 2T  ∂T ∂T ∂T +u +v =α 2 + 2  ∂t ∂x ∂y ∂y   ∂x Immediately after imposing the temperature difference ( t = 0+ ), the fluid is still motionless ( u = v = 0 ), hence the energy equation (6.208) reflects the balance between the thermal inertia and the conduction in the fluid. The scales of the two terms enclosed in the parentheses on the right-hand side of eq. (6.208) are ΔT / δ t2 and ΔT / H 2 , respectively. Since δ t  H , one can conclude that ∂ 2T / ∂y 2  ∂ 2T / ∂x 2 . The balance of scales for eq. (6.208) then becomes: ΔT ΔT α 2 t δt Thus, the scale of the thermal boundary layer thickness becomes: δ t ~ (α t )1/ 2 (6.209) Chapter 6 Natural Convection 557 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press To estimate the scale of the velocity, one can combine eqs. (6.206) and (6.207) by eliminating the pressure to obtain: ∂  ∂v ∂v ∂v  ∂  ∂u ∂u ∂u  +u +v   +u +v −  ∂x  ∂t ∂x ∂y  ∂y  ∂t ∂x ∂y   ∂  ∂ 2 v ∂ 2 v  ∂  ∂ 2u ∂ 2u  ∂T = ν   2 + 2  −  2 + 2   + gβ ∂y  ∂y  ∂x ∂y   ∂x  ∂x  ∂x (6.210) where the left-hand side represents the inertia terms, and the right-hand side represents the viscosity and buoyancy terms. The scales of these three effects are shown below Inertia v δt t Viscosity v δ δt To examine the relative strength of each effect, one can divide the above expression by the scale of viscosity effect to obtain 3 t ν Buoyancy gβ ΔT (6.211) Inertia 1 Pr Viscosity 1 Buoyancy gβΔT δ t2 νv where eq. (6.209) was used to simplify the inertia term. For the fluid with Pr > 1 , the momentum balance at t = 0+ requires a balance between the viscosity and buoyancy terms: gβ ΔT δ t2 1~ νv Substituting eq. (6.209) into the above expression and rearranging the resultant expression, the scale of vertical velocity at the initiation of the natural convection is obtained as follows: g β ΔT α t v~ (6.212) ν Let us turn our attention now to the energy equation (6.208). The first term on its left-hand side is the inertia term, and its scale is ΔT / t . The second and third terms on the left-hand side are the advection terms and both of them have scale vΔT / H . The right-hand side of eq. (6.208) is the conduction term, and its scale is αΔT / δ t2 . As time increases, the effect of the inertia term weakens, hence the effect of advection becomes stronger. This trend continues until a final time, tf, when the energy balance requires balance between the advection and conduction terms, i.e., v ΔT ΔT ΔT ~α 2 ~ H tf δt , f 1/ 2 Thus, the scale of tf becomes  νH  tf ~    gβΔT α  (6.213) 558 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The scale of thermal boundary layer thickness at time tf is: δ t , f ~ (α t f )1/ 2 ~ H Ra −1/ 4 (6.214) H At time tf, natural convection in the rectangular enclosure reaches steady-state and the thickness of the thermal boundary layer no longer increases with time. As the fluid in the thermal boundary layer near the hot wall is heated and rises due to buoyancy force, it drags the fluid outside of the thermal boundary layer upward forming a wall jet. Similar to the thermal boundary layer thickness, the wall jet thickness also increases with time until t = tf when the maximum wall jet thickness, δv,f, is reached (see Fig. 6.19). Outside the thermal boundary layer, the buoyancy force is absent and the thickness of the wall jet can be determined by balancing the inertia and viscosity terms in eq. (6.210): ~ν 3 δ vt δv which can be rearranged to obtain: δ v ~ (ν t )1/ 2 ~ Pr1/ 2 δ t (6.215) For t > tf , steady-state has been reached, and the wall jet thickness is related to the thermal boundary layer thickness by δ v, f ~ Pr1/ 2 δ t , f . v v Similarly, the condition to have distinct vertical wall jets or momentum boundary layers is δ v, f < L , or equivalently: H <Ra1/ 4 Pr −1/ 2 H L (6.216) T 0 v αt δt, f x 0 νt δ v, f x Figure 6.19 Two-layer structure near the heated wall. Chapter 6 Natural Convection 559 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press When the vertical wall jet encounters the horizontal wall, it will turn to the horizontal direction and become a horizontal jet. This horizontal jet will contribute to the convective heat transfer from the heated wall to the cooled wall: ′ qconv ~ ρ v f δ T , f c p ΔT Considering eqs. (6.212) and (6.214), the above scale of convective heat transfer becomes: ′ qconv ~ k ΔT Ra1/ 4 H When a warm jet is formed at the top and a cold jet is formed at the bottom, there will be a temperature gradient along the vertical direction. The heat conduction due to this temperature gradient is: ′ qcond ~ kL ΔT H The condition under which the horizontal wall jets can maintain their temperature identity is when the heat conduction along the vertical direction is negligible compared to the energy carried by the horizontal jets: kL ΔT < k ΔT Ra1/ 4 H H H > Ra −1/ 4 H L Table 6.2 Characteristics of natural convection in a rectangular enclosure heated from the left side Regimes Condition of occurrence Flow pattern Effect of flow on heat transfer Heat transfer mechanism Heat transfer I: Conduction II: Tall Systems III: Boundary layer IV: Shallow systems or equivalently (6.217) The characteristics of various heat transfer regimes are summarized in Table 6.2. Ra H < 1 Clockwise circulation Insignificant Conduction in horizontal direction H / L >Ra 1/ 4 H Ra −1/ 4 < H / L <Ra1/ 4 H H Boundary layer on all four walls. Core remains stagnant Significant Boundary layer convection H / L <Ra −1/ 4 H Two horizontal wall jets flow in opposite directions. Significant Conduction in vertical direction Distinct boundary layer on top and bottom walls Insignificant Conduction in horizontal direction q′  kH ΔT / L q′  kH ΔT / L q′  (k / δ T , f ) H ΔT q′  (k / δ T , f ) H ΔT 6.5.2 Rectangular Enclosures Heated from the Side Heat transfer in regimes I and II, shown in Table 6.2, is dominated by conduction and the fluid circulation plays an insignificant role. Therefore, heat transfer in these two regimes is practically equal to heat transfer rate estimated using Fourier’s law. 560 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (a) (b) (c) (d) Figure 6.20 Natural convection in a square enclosure (RaH =106, Pr=0.71). Steady-state distributions of dimensionless (a) temperature, (b) pressure, (c) horizontal velocity, and (d) vertical velocity (Wang et al., 2010) The flow pattern in regime III is of the boundary layer type and the core of the fluid is stagnant and stratified. The fluid flow and heat transfer in regime III can be obtained via boundary layer analysis. Figure 6.20 shows contours of dimensionless temperature, pressure, horizontal, and vertical velocities for natural convection in a square enclosure ( H / L = 1 ) heated from its left side and cooled from its right side while insulated from the top and bottom. The behavior of boundary layer type flow and heat transfer is apparent because both temperature and velocity gradient peak near the heated and cooled wall. The low gradients of temperature and velocity in the middle of the enclosure further indicate a stagnant and stratified core. For high Rayleigh number, the average Nusselt number for regime III, as obtained from boundary layer analysis, is (Bejan, 2004): Nu H = 0.364 L Ra1/ 4 H H (6.218) which indicates that natural convection heat transfer in a rectangular enclosure is dominated by both the Rayleigh number and aspect ratio. Chapter 6 Natural Convection 561 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press For laminar natural convection in a rectangular enclosure with aspect ratio H/L between 1 and 10, the following correlations were suggested by Berkovsky and Polevikov (1977):  Pr  Nu H = 0.22  Ra H   0.2 + Pr  which is valid for 2 < H / L < 10, Pr < 105 , 0.28 L  H Ra H < 1013 , and −0.13 0.09 (6.219) Catton (1978) found that the above correlations agreed well with the experimental data and the results from numerical solutions. For a rectangular enclosure with a larger aspect ratio, MacGregor and Emery (1969) recommended the following correlations: H Nu L = 0.42Ra Pr (6.221)  L which is valid for 10 < H / L < 40, 1< Pr < 2 × 104 , 104 <Ra L < 107 , and 1/ 4 L 0.012 −0.3  Pr  L (6.220) Nu H = 0.18  Ra H     0.2 + Pr  H which is valid for 1 < H / L < 2, 10-3 < Pr < 105 , 103 < [Pr/ (0.2 + Pr)]Ra H ( L / H )3 . 0.29 (6.222) which is valid for 1 < H / L < 40, 1< Pr < 20, 10 < Ra L < 10 . It should be pointed out that the characteristic length for the Rayleigh and Nusselt numbers used in eqs. (6.221) and (6.222) is the width of the enclosure, L. For the shallow enclosure represented by regime IV in Table 6.2, the horizontal wall jet flows from the left to right on the top wall and from right to the left on the bottom wall. Figure 6.21 shows streamlines and isotherms for natural convection in a shallow enclosure ( H / L = 0.1 ) heated from the right side Nu L = 0.046Ra1/3 L 6 9 Figure 6.21 Natural convection in a shallow enclosure: (a) RaH = 106, (b) RaH = 108 (H/L = 0.1, Pr = 1.0; Tichy and Gadgil, 1982). 562 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press and cooled from the left side (Tichy and Gadgil, 1982). It can be observed from these figures that there exist wall jets on both the top and bottom walls and that the high velocity in these wall jets drives the core flow. The entire enclosure can be divided into two end regions near the heated and cooled vertical walls and a core region between the two end regions. There are large temperature gradients along the horizontal direction in the end regions, but the core region is essentially stratified. Under fixed aspect ratio, increasing Rayleigh number makes the horizontal wall jets thinner. The core region extends almost the entire length of the shallow enclosure for high Rayleigh number. Bejan and Tien (1978) performed a scale analysis and obtained an analytical solution via analysis of fluid flow and heat transfer in the core region. As indicated by Fig. 6.22, their results compared favorably with experimental (Imberger, 1973) and numerical (Cormack et al., 1974) results. Bejan (2004) pointed out that Fig. 6.22 is not only useful for regime IV of Table 6.2 (shallow enclosure), but also for regime III (boundary layer regime). The success of Bejan and Tien’s analysis was attributed to the fact that the thermal resistances associated with the vertical end walls were taken into account. For a limited case, the asymptotic heat transfer results can be expressed as: 1 H  (6.223)  Ra H  362880  L  which is valid when ( H / L) 2 Ra H → 0 , and is shown as the dashed line in Fig. Nu H = 1 + 2 6.22. Eq. (6.223) Figure 6.22 Average Nusselt number for natural convection in a shallow enclosure (Bejan and Tien, 1978). Chapter 6 Natural Convection 563 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Heated from the Bottom For a rectangular enclosure filled with fluid and heated from the side, natural convection will be initiated as soon as the temperature difference between the two vertical walls is established. For a rectangular enclosure heated from below, natural convection may or may not occur depending on whether the temperature difference between the top and bottom walls exceeds a critical value. If the temperature difference is not sufficient to initiate natural convection, the mechanism of heat transfer is by conduction of the fluid (and radiation for the case where the fluid is a gas). The condition for the onset of natural convection can be expressed in terms of a critical Rayleigh number. For the case that the horizontal dimension is much larger than the height of the enclosure, the criterion for the onset of natural convection is: gβΔTH 3 RaH = > 1708 (6.224) να When the Rayleigh number just exceeds the above critical Rayleigh number, the flow pattern is two-dimensional counter rotating rolls – referred to as Bénard cells [see Fig. 6.23 (a)]. As the Rayleigh number further increases to one or two orders of magnitude higher than the above critical Rayleigh number, the twodimensional cells breakup to three dimensional cells whose top view is hexagonal [see Fig. 6.23(b)]. The function of the two-dimensional rolls and three-dimensional hexagonal cells is to promote heat transfer from the heated bottom wall to the cooled top wall. Globe and Dropkin (1959) suggested the following empirical correlation Nu H = 0.069Ra1/3 Pr 0.074 (6.225) H Figure 6.23 Rolls and hexagonal cells in natural convection in enclosure heated from below (Oosthuizen and Naylor, 1999). 564 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where all thermophysical properties are evaluated at (Th + Tc ) / 2 . Equation (6.225) is valid for 3 × 105 < Ra H < 7 × 109 . In addition, H/L must be sufficiently small so that the effect of the sidewalls can be negligible. Inclined Rectangular Enclosure When the rectangular enclosure heated from the side is tilted relative to the direction of gravity, additional unstable stratification and thermal instability, which are usually associated with natural convection in an enclosure heated from below, will affect the fluid flow and heat transfer. The tilt angle γ [see Fig. 6.17 (c)] becomes an additional parameter that plays a significant role. When γ is increased from 0° to 90°, the heat transfer mechanism changes from pure conduction at γ = 0 (heated from the top) to single-cell convection at γ = 90 (heated from the side), where both buoyancy effect and Nusselt number are maximized. As the inclination angle γ is further increased beyond 90°, the Nusselt number starts to decrease until a local minimum at γ = γ c is reached. The value of γ c depends on the aspect ratio of the enclosure (see Table 6.3). Beyond γ c , the stable two-dimensional single-cell convection becomes unstable and transits to three-dimensional rolls. The mechanism of Bénard convection sets in, which results in an increase of the Nusselt number. The variation of Nusselt number as a function of tilt angle is qualitatively shown in Fig. 6.24. Table 6.3 Critical inclination angle for different aspect ratio (Arnold et al., 1976) Aspect ratio, H/L 1 3 6 Critical tilt angle, γc 155° 127° 120° 12 113° >12 110° Nu H 0 90° γc 180° γ Figure 6.24 Effect of inclination angle on natural convection in enclosure Chapter 6 Natural Convection 565 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Isotherms Streamlines Figure 6.25 Natural convection in inclined squared enclosures (Zhong et al, 1983). Zhong et al. (1983) numerically solved the natural convection in an inclined square enclosure, and the isotherms and the streamlines for Ra = 105 are shown in Fig. 6.25. The transition from a conduction dominated regime at γ = 45 to single-cell convection at γ = 90 is apparent from the isotherms. At γ = 135 , which, according to Table 6.3, is less than the critical inclination angle, the isotherms start to exhibit some behaviors of thermally unstable conditions. They proposed the following correlation for natural convection of air in a square enclosure ( H / L = 1 ) in the region 0 < γ < 90 : Kγ = Nu H (γ ) − Nu H (0 ) 2 = γ sin γ Nu H (90 ) − Nu H (0 ) π (6.226) where Nu H (0 ) is for pure conduction. While eq. (6.226) is good for air in a squared enclosure, the following correlation can be applied to other situations: 566 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  L   1 +  H Nu H (90 ) − 1 sin γ L   Nu H (γ ) =   H L  Nu H (90 )(sin γ )1/ 4 H  0 < γ < 90 90 < γ < γ c  (6.227) which indicates that the Nusselt number for an inclined enclosure is related to the Nusselt number at γ = 90 and the actual inclination angle. Example 6.5 A rectangular cavity is formed by two parallel plates, each with a dimension of 0.5 m by 0.5 m, which are separated by a distance of 5 cm. The temperatures of the two plates are 37 °C and 17 °C, respectively. Find the heat transfer rate from hot plate to cold plate for the inclination angles of 0°, 45°, 90°, and 180°. Solution: The temperatures of the two plates are Th = 37  C and Tc = 17  C , respectively. The dimensions of the enclosure as shown in Fig. 6.26, are H = 0.05 m and L = 0.5 m, respectively. The depth of the enclosure is D = 0.5 m. For all heat transfer calculations, the thermophysical properties of the air will be evaluated at the mean temperature Tm = (Th + Tc ) / 2 = 27 o C = 300K . The required properties can be obtained from Table C.1: ν = 15.89 × 10−6 m 2 /s , k = 0.0263 W/m-K, α = 22.5 × 10−6 m 2 /s and Pr = 0.707. At the inclination angle 0°, we have heat transfer in a rectangular cavity heated from the above. Natural convection is absent in this case, and heat transfer is solely by conduction. The rate of heat transfer is q(0o ) = kDH Th − Tc 37 − 17 = 0.0263 × 0.5 × 0.5 × = 2.63W 0.05 L L H Tc Th Air γ Figure 6.26 Natural convection in inclined squared enclosure. Chapter 6 Natural Convection 567 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The Nusselt number for the case of pure conduction is Nu L (0o ) = 1 or Nu H (0o ) = ( H / L )Nu L (0o ) = 10 . Heat transfer for the case of γ = 45o can be calculated using eq. (6.227), which requires the Nusselt number for γ = 90o . Thus, the heat transfer rate for γ = 90o will be calculated first. The Rayleigh number is: gβ (Th − Tc ) H 3 9.807 × 1 / 300 × (37 − 17) × 0.53 Ra H = = = 2.286 × 108 να 15.89 × 10−6 × 22.5 × 10 −6 The Nusselt number for γ = 90o can be obtained from eq. (6.220):  Pr  Nu H (90o ) = 0.18  Ra H   0.2 + Pr  0.28 L  H −0.13  0.707  = 0.18 ×  × 2.286 × 108   0.2 + 0.707  0.28  0.05     0.5  −0.13 = 49.6 The heat transfer coefficient is: k Nu H (90o ) 0.0263 × 49.6 = = 2.61 W/m 2 -K H 0.5 Therefore, the heat transfer rate for γ = 90o is: h (90o ) = q(90o ) = h(90o ) DH (Th − Tc ) = 2.61× 0.5 × 0.5 × (37 − 17) = 6.53 W When the inclination angle is γ = 45o = π / 4 , eq. (6.227) yields: L L  Nu H (45o ) = 1 +  Nu H (90 ) − 1 sin γ H H   0.05  = 1+  × 49.6 − 1 sin 45o = 3.80 0.5   Thus, the Nusselt number is Nu H (45o ) = 38.0 and the corresponding heat transfer coefficient is: k Nu H (45o ) 0.0263 × 38.0 = = 2.00 W/m 2 -K H 0.5 The heat transfer rate for γ = 45o is therefore: h (45o ) = q(45o ) = h (45o ) DH (Th − Tc ) = 2 × 0.5 × 0.5 × (37 − 17) = 5.00W When the inclination angle is γ = 180o , the problem becomes natural convection in an enclosure heated from the below. The Rayleigh number is: gβ (Th − Tc ) L3 9.807 × 1 / 300 × (37 − 17) × 0.053 Ra L = = = 2.286 × 105 να 15.89 × 10−6 × 22.5 × 10−6 The Nusselt number in this case can be obtained from eq. (6.225): Nu L (180o ) = 0.069Ra1/3 Pr 0.074 H = 0.069 × (2.286 × 105 )1/3 × 0.707 0.074 = 4.11 The heat transfer coefficient is 568 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press k Nu L (180o ) 0.0263 × 4.11 = = 2.16 W/m 2 -K L 0.05 The heat transfer rate for γ = 180o is therefore, h (180o ) = q(180o ) = h(180o ) DH (Th − Tc ) = 2.16 × 0.5 × 0.5 × (37 − 17) = 10.80 W 6.5.3 Annular Space between Concentric Cylinders and Spheres Natural convection in spaces between long horizontal concentric cylinders [see Fig. 6.17 (e)] or between spheres is very complicated. The only practical approach to analyze the problem is via numerical solution. The physical model of natural convection in annular space between concentric cylinders is shown in Fig. 6.27. Since the problem is axisymmetric, one only needs to study the right half of the domain. In the coordinate system shown in Fig. 6.27, the governing equations are 1 ∂u ∂v v + + =0 r ∂θ ∂r r  1 ∂ 2 u 1 ∂u ∂ 2 u u 2 ∂v  ∂u uv u ∂u 1 ∂p +v + =− +ν  2 + + −+  2 ∂r r r ∂θ r ∂r ∂r 2 r 2 r 2 ∂θ  ρ r ∂θ  r ∂θ − gβ (T − Tref ) sin θ  1 ∂ 2 v 1 ∂v ∂ 2 v v 2 ∂u  ∂v u 2 u ∂v 1 ∂p +v − =− +ν  2 + + −−  2 ∂r r r ∂θ r ∂r ∂r 2 r 2 r 2 ∂θ  ρ ∂r  r ∂θ +gβ (T − Tref ) sin θ (6.228) (6.229) (6.230) δ ro θ, u r, v ri g Ti To Figure 6.27 Natural convection in horizontal annular space between concentric cylinders Chapter 6 Natural Convection 569 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  1 ∂ 2T ∂ 2T 1 ∂T  ∂T u ∂T +v =α 2 + 2+  2 ∂r ∂r r ∂θ r ∂r   r ∂θ (6.231) where Tref is a reference temperature and: (6.232) peff = p + ρ g sin θ − ρ g cos θ is the effective pressure. Equations (6.228) – (6.231) are subject to the following boundary conditions: u = 0, (6.234) u = v = 0, T = To , at r = ro (6.235) Date (1986) solved this problem numerically using a modified SIMPLE algorithm for δ / Di = 0.8 and 0.15 . Figure 6.28 shows the contours of the stream functions and isotherms at Ra = 4.7 × 104 . The dimensionless stream function, dimensionless temperature and Rayleigh number are defined as the following T − To g β ΔT δ 3 ψ Ψ= , Θ= , Ra = (6.236) α να Ti − To where the stream function, ψ , is defined as ∂ψ 1 ∂ψ u=− , v= (6.237) ∂r r ∂θ It can be seen that isotherms concentrate near the lower portion of the surface of the inner cylinder and the upper portion of the surface of the outer cylinder, which are indications of the development of thermal boundaries near the heated and cooled surfaces. While the contours of the streamlines have kidney-like shapes with the center of the flow rotation moves upward due to effect of natural convection. The streamlines and isotherms shown in Fig. 6.28 agreed very well ∂v ∂T = = 0, at θ = 0 or π ∂θ ∂θ u = v = 0, T = Ti , at r = ri (6.233) Ψ Θ Figure 6.28 Natural convection in a horizontal annulus (Pr =0.7, δ/Di=0.8, Ra =4.7×104; Date, 1986) 570 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press with results obtained by Kuehn and Goldstein (1976) using vortex-stream function method. The rate of heat transfer per unit length of the annulus can be calculated by the following correlation (Raithby and Hollands, 1975; Bejan, 2004): 2.425k (Ti − To )  Pr Ra Di  q′ ≅   [1 + ( Di / Do )3/5 ]5/ 4  0.861 + Pr  1/ 4 (6.238) where the Rayleigh number is defined as gβ (Ti − To ) Di3 Ra D = (6.239) να Equation (6.238) is valid for 0.7 < Pr < 6000 and Ra < 107 . The thermophysical properties of the fluid should be evaluated at the mean temperature (Ti + To ) / 2 . The scales of the thermal boundary layer on the inner surface of the outer cylinder and on the outer surface of the inner cylinder are: − − δ o ~ Do RaD1/ 4 , δ i ~ Di RaD1/ 4 i o i It is obvious that δ o > δ i since Do > Di. Equation (6.238) will be valid only if the boundary layer thickness is less than the gap between the two cylinders, i.e. only if δo < Do – Di. Under lower Rayleigh numbers, on the other hand, we have Do Ra −1/ 4 > Do − Di (6.240) D o and the heat transfer mechanism between two cylinders will approach pure conduction. In this case, the heat transfer rate obtained from eq. (6.238) may be less than that obtained from the pure conduction model, which does not make sense. Instead of using eq. (6.240) to check the validity of eq. (6.238), an alternative method is to calculate the heat transfer rate via eq. (6.238) and pure conduction model, and the larger of the two heat transfer rates should be used. For natural convection in the annulus between two concentric spheres, the trends for the evolution of the flow pattern and isotherms are similar to the concentric cylinder except the circulation between concentric spheres has the shape of a doughnut. The empirical correlation for the heat transfer rate is (Raithby and Hollands, 1975; Bejan, 2004): 2.325kDi (Ti − To )  Pr Ra Di  q≅   [1 + ( Di / Do )7 /5 ]5/ 4  0.861 + Pr  1/ 4 (6.241) where the definition of Rayleigh number is the same as for eq. (6.239). Equation (6.241) is valid for 0.7 < Pr < 4000 and Ra < 104. The thermophysical properties of the fluid should be evaluated at the mean temperature, (Ti + To ) / 2 . Similar to the annulus between two cylinders, one should also make use of both eq. (6.241) and the pure conduction model and retain the larger of the two values obtained. Chapter 6 Natural Convection 571 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 6.6 Natural Convection in Melting and Solidification Interest in utilizing clean energy is growing due to environmental considerations. Solar energy is one of the most promising clean energies, but due to its periodic nature, a heat storage device is needed. A thermal energy storage system using a Phase Change Material (PCM) is an attractive method because a large amount of latent heat can be stored or released during the melting or solidification process of the PCM. During the thermal energy storage process, the hot working fluid from the solar energy collector flows into the thermal energy storage device to melt the PCM and store thermal energy in the form of latent heat. During the nighttime, the PCM solidifies and the latent heat is released. Natural convection plays a significant role in both melting and solidification processes due to the temperature difference in the liquid phase. While melting is promoted by the presence of natural convection, solidification is slowed down by its presence. In this section, natural convection during solidification around a horizontal tube (Zhang et al., 1997) and during melting in an enclosure heated from the side (Zhang, 2002) will be discussed. 6.6.1 Solidification around Horizontal Cylinder The theoretical model to be employed in this study is shown in Fig. 6.29. At the very beginning of the process (t = 0), the tube, which has a radius of R0, is surrounded by a liquid phase change material with uniform temperature Tf > Tm. The temperature of the working fluid inside the tube is Ti, and the convective heat transfer coefficient between the working fluid and the internal tube wall is hi. Both hi and Ti are kept constant throughout the process. The thickness of the tube is assumed to be very thin, so the thermal resistance of the tube wall can be neglected. The phase change material can be treated as if it were directly in contact with the coolant inside the tube. The liquid adjacent to the cooled surface R0 R θ hi, Ti δ x g y Figure 6.29 Solidification around a horizontal tube. 572 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press will be solidified, and the temperature difference between the solid-liquid interface and the otherwise quiescent liquid will drive natural convection in the liquid region. The liquid region is assumed to be sufficiently large so that the convection can be treated as natural convection between the surface and an extensive fluid medium. Conduction is the only heat transfer mode in the solidified region. The analysis is based on several additional assumptions: 1. The liquid is Newtonian and Boussinesq as well as incompressible. The Prandtl number of the liquid phase change material is greater than unity. 2. The solid is homogeneous and isotropic. 3. The liquid motion induced by volumetric variation during solidification is neglected, i.e., the density of the solid is equal to the density of the liquid. In addition, the properties of the phase change material are constant in the liquid and solid regions. 4. The solid-liquid interface is assumed to be a smooth cylinder concentric with the cooled tube. 5. The effects of natural convection are restricted to within the boundary layer and the bulk liquid has a uniform temperature, Tf. Compared to the thickness of the solidified layer, the thickness of the natural convective boundary layer on the solid-liquid interface is very thin, except at the onset of freezing (the boundary layer thickness in Fig 6.29 is significantly exaggerated so as to be clearly visible). However, since the onset of freezing is a very short period compared to the whole solidification process, it is reasonable to neglect the curvature effect in the equations of the boundary layer. The boundary layer equations of the problem can be written as follows: ∂u ∂v + =0 ∂x ∂y (6.242) (6.243) (6.244) ν ∂2u + gβ (T − Tm ) sin θ = 0 ∂y 2 ∂T ∂ (uT ) ∂ (vT ) ∂ 2T + + =αf ∂t ∂x ∂y ∂y 2 Since the Prandtl number of the liquid is much larger than unity, the inertia terms in the momentum equation have been neglected. These equations were solved by an integral method. Assuming a polynomial profile, the temperature and velocity profiles inside the boundary layer of thickness δ are expressed as  y T = T f − (T f − Tm ) 1 −   δ u =U y y 1− δ δ   2 2 (6.245) (6.246) where U is a characteristic velocity, δ is the boundary layer thickness, and both of them are functions of time t and angle θ. Integrating both sides of eq. (6.243) Chapter 6 Natural Convection 573 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press with respect to y within the interval (0, δ ), the expression for U in eq. (6.246) should be gβ sin θ U= (T f − Tm )δ 2 (6.247) 3ν Integrating eq. (6.244) in the same manner, one obtains 1 ∂δ dR 1 gβ sin θ ∂ 3δ 2α f − + − (T f − Tm ) + =0 (6.248) 3 ∂t dt 90 ν R ∂θ 3 δ Heat transfer in the solidified layer is dominated by conduction. The governing equation of the solid layer and the corresponding boundary conditions are as follows: 1 ∂  ∂T r r ∂r  ∂r  1 ∂T =  α s ∂t R0 < r < R t > 0 (6.249) (6.250) (6.251) (6.252) ks ∂T = hi (T − Ti ) r = R0 t > 0 ∂r Ts (r, t ) = Tm r = R t > 0 − kf ∂T ∂y = ρ hs y =0 ks ∂Ts ∂r r=R dR dt r=R t >0 Equations (6.248) – (6.252) can be nondimensionalized and solved using an integral approximation method. When Bi → ∞ , it corresponds to boundary conditions of the first kind, i.e. the tube wall temperature is Ti and is kept steady throughout the process. Wang et al. (1991) experimentally investigated the solidification process around a horizontal cooled tube. A comparison of the predicted solidification rate, V / V0 = ( R / R0 ) 2 and the experimental results is shown in Fig. 6.30. When Ra = 0, i.e., no superheat exists in the liquid region or the solidification process is dominated by conduction, the predicted value is 18% Figure 6.30 Comparison of predicted solidification rate with experiments (Zhang et al., 1997). 574 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press lower than the experimental data. During the conduction-dominated freezing process, the front of the freezing layer is a dendritic layer. Therefore, it is believed that the solid-liquid interface is extended by the dendritic layer. As the Rayleigh number increases, natural convection occurs in the liquid region, and the solid-liquid interface becomes smooth because of the natural motion of the liquid. When Ra = 1.8 × 105 , the predicted value is only 8% lower than the experimental data, so the agreement is satisfactory. 6.6.2 Melting in a Rectangular Enclosure Heated from the Side Natural convection controlled melting in a rectangular cavity under boundary condition of the first and second kinds has been investigated extensively (Ho and Viskanta, 1984; Gadgil and Gobin, 1984; Bejan, 1989; Cao and Faghri, 1990; Zhang and Bejan, 1989; Zhang and Chen, 1994; Wang et al., 2010). In a solar energy utilization system, the working fluid absorbs heat in the solar energy collector and transfers this heat to the latent heat thermal energy storage system via convection. Therefore, melting in the latent heat thermal energy storage system is under boundary condition of neither the first kind nor the second kind. It is instead under boundary condition of the third kind. In order to thoroughly understand the mechanism of the latent heat thermal energy storage unit, it is necessary to study melting in an enclosure under the boundary condition of the third kind. The physical model under consideration is shown in Fig. 6.31. The enclosure with height of H and width of L is filled with a solid PCM at its melting point, Tm. The top, bottom and the left side of the enclosure are adiabatic. At time t > 0, the right side of the enclosure is exposed to the working fluid at temperature of Te. The convective heat transfer coefficient between the working fluid and the right wall is constantly equal to he. It is assumed the thickness of the right-side wall is very small and therefore, the conduction thermal resistance of the wall is negligible. The melting process can be divided into two stages: conduction and convection stages. During the conduction stage, conduction is the dominant mode of heat transfer. The problem can be simplified as a 1-D melting problem under boundary condition of the third kind and the analytical solution is available (see Chapter 3). The melting process enters the convection stage when the thickness of the melting layer increases and natural convection becomes the dominant mode of heat transfer. Although the process is still unsteady, the transient terms in the governing equations describing the melting process have negligible effect and the process can be considered to be quasi-steady (Bejan, 1989). When the melting is in the early stage of the convection regime, the solidliquid interface can be assumed to be sufficiently straight and vertical. In addition, the following assumptions are made: Chapter 6 Natural Convection 575 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press y y λ T = Tm v H u Te he δ x L 0 xr Figure 6.31 Melting in an enclosure heated from the side (Zhang, 2002). 1. The properties of liquid PCM are constants except for density, which is a function of temperature (Boussinesq assumption). The liquid is Newtonian and incompressible. 2. The Prandtl number of the liquid phase change material is much greater than 1. This assumption is valid for most PCMs (Zhang and Bejan, 1989). 3. The liquid motion induced by volumetric variation during melting is neglected, i.e. the density of the solid is equal to that of the liquid. The liquid phase can be divided into three regions: (1) a cold boundary layer near the solid-liquid interface, (2) a warm boundary layer near the heated vertical wall, and (3) the core region between these two boundary layers (Bejan, 1989). The analysis of the cold boundary layer is performed in a frame (x-y system) that is attached to the solid-liquid interface. Since the Prandtl number of the liquid is much larger than unity, the inertia terms in the momentum equation can be neglected. The advection induced by the solid- liquid interface motion in the boundary layer equations can also be neglected since the solid-liquid interface velocity is typically several orders of magnitude smaller than the natural convective velocity in the boundary layer (Gadgil and Gobin 1984). The momentum and energy equations for the cold boundary layer are: ∂ 3 v gβ ∂T + =0 (6.253) ∂x 3 ν ∂x 576 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u ∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂x (6.254) which can be solved by an integral method. By following the procedure similar to Bejan (1989) and Zhang and Bejan (1989), one obtains an ordinary differential equation for cold boundary layer thickness, δ: 2k (Tc − Tm ) 2 gβ (Tc − Tm ) d δ 3 (Tc − Tm )  +  ρ hsδ 36ν dy  T − Tm gβ d δ 3 (Tc − Tm ) 2  = −2α c − (6.255)   60ν dy δ where Tc is the core temperature. Since the cold boundary layer starts from the top of the enclosure, the boundary condition of eq. (6.255) is δ (H ) = 0 (6.256) The solid-liquid interface velocity, u0, can be obtained by energy balance at the interface. Considering the assumed temperature profile in the cold boundary layer used in the integral solution, the solid-liquid interface velocity is u0 = 2k (Tc − Tm ) ρ hsδ (6.257) The warm boundary layer is attached to the heated wall, i.e. on to a stationary frame (xr-y). The analysis of the warm boundary layer is very similar to that of the cold boundary layer. The momentum and the energy equations are the same as eqs. (6.253) and (6.254) except (x, v, T) are replaced with (xr, vr, Tr). After performing the integral solution, the equation of the warm boundary layer thickness becomes: gβ (Tw − Tc )λ 3 dTc T − Tc gβ d  λ 3 (Tw − Tc ) 2  + = −2α w (6.258)   90ν dy 36ν dy λ which is subject to the following boundary condition: λ (0) = 0 (6.259) The wall temperature, Tw, can be obtained by considering energy balance at the right side wall of the enclosure: (6.260) = he (Te − Tw ) λ The temperature of the core region, Tc, can be obtained by considering the fact that the horizontal entrainment velocities of the two boundary layers represent the same core flow (Bejan, 1989; Zhang and Bejan, 1989). δ 3 (Tc − Tm ) = λ 3 (Tw − Tc ) (6.261) The core temperature is the lowest at the bottom and highest at the top of the enclosure: Tc (0) = Tm (6.262) Tc ( H ) = Tw ( H ) (6.263) Defining the following dimensionless variables: 2k (Tw − Tc ) Chapter 6 Natural Convection 577 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Ste = c p (Te − Tm ) hs , Ra * = gβ (Te − Tm ) H 3 να ,Y= 1 1 y δ λ , Δ = Ra * 4 , Λ = Ra * 4 H H H T − Tm hH u0 H ρ hs , Ra = θ w Ra * , Bi = e , U 0 = θ= Te − Tm k k (Te − Tm ) Ra*1/ 4 (6.264) the following dimensionless equations are obtained: θ2 θ d 2θ 1d 2Ste c + c (Δ 3θ c ) − (Δ 3θ 2 c ) = − c Δ 36 dY 60 dY Δ 3 Λ (θ w − θ c ) dθ c 2(θ w − θ c ) 1d  Λ 3 (θ w − θ c ) 2  + =−  90 dY  36 Λ dY 3 3 Δ θ c = Λ (θ w − θ c ) 2(θ w − θ c ) = BiRa*−1/ 4 (1 − θ w ) Λ (6.265) (6.266) (6.267) (6.268) and their corresponding boundary conditions are: Y =0 Λ=0 θc = 0 (6.269) Y =1 Δ = 0 θc = θ w (6.270) The dimensionless solid-liquid interface velocity is: U 0 = 2θ c / Δ (6.271) The set of ordinary differential equations (6.265) and (6.266), with boundary conditions specified by eq. (6.269) and (6.270) is a boundary value problem. The boundary value problem was solved by using a shooting method in the interval of Y = 0 and Y = 1. The objective of the shooting method was to satisfy boundary condition Δ(1) = 0 by properly adjusting Δ(0). Equations (6.265) and (6.266) were discretized using an implicit scheme. The discretized equations and eqs. (6.267) and (6.268) are solved using an iteration method. The converged solution for the entire problem was obtained when shooting error for Δ(1) was less than 10-9. 1.0 1.0 0.8 0.8 0.6 0.6 Y 0.4 Y 0.4 Bi/Ra 1/4=0.1 * * * 0.2 Bi/Ra 1/4 * * * =0.1 Bi/Ra 1/4=1.0 0.2 Bi/Ra 1/4=1.0 0.0 0.0 0.5 1.0 Bi/Ra 1/4=10 4 Bi/Ra 1/4=104 0.0 2.0 1.5 Bejan(1989) 0.0 0.5 1.0 1.5 2.0 θw U0 (a) (b) Figure 6.32 Dimensionless wall temperature and interfacial velocity (Zhang, 2002). 578 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Figure 6.32(a) shows the distribution of the dimensionless wall temperature at different external convective heat transfer conditions. The Stefan number is usually very small for the solid-liquid phase change system, and it has a negligible effect on the results (Bejan, 1989). Thus, the results in Fig. 6.32 were obtained for the case of Ste = 0. The dimensionless wall temperature is a constant when Bi/Ra*1/4 = 104, which agrees with the case of the boundary condition of the first kind (Tw = Te). The dimensionless wall temperature increases with increasing Y for lower Bi/Ra*1/4. The wall temperature increases rapidly near Y= 0 and Y = 1, and it is a linear function of Y in the middle portion along the height direction. When Bi/Ra*1/4 = 0.1, the result is similar to that of constant heat flux heating in Zhang and Bejan (1989). The variation of the dimensionless solidliquid interface velocity along the Y direction for different external heat transfer conditions is shown in Fig. 6.32(b). For the case of Bi/Ra*1/4 = 104, the present results agreed very well with the results of Bejan (1989), who investigated melting in an enclosure under boundary conditions of the first kind. As can be seen from Fig. 6.32(b), the solid-liquid interface velocity decreases with decreasing Bi/Ra*1/4, since the heat absorbed by the PCM decreases with decreasing Bi/Ra*1/4. The heat transfer characteristics and the rate of melting are shown in Fig. 6.33. The heat transfer characteristics is evaluated by using the following parameter: 2 (θ w − θ c ) Nu (6.272) = 1/ 4 Λθ w Ra* where Nu = hL/k is Nusselt number. It should be pointed out that the temperature difference in Ra* is (Te−Tm), which is really not the driving force of natural convection. The driving force for natural convection in the liquid phase should be (Tw−Tm). Therefore, the Rayleigh number based on (Tw−Tm), Ra, is defined and its 0.5 0.4 0.3 0.2 0.1 R Nu/Ra*1/4 Nu/Ra1/4 0.0 0.1 1 10 100 1000 10000 Bi/Ra*1/4 Figure 6.33 Heat transfer and Fig. 5 Heat transfer and melting rate melting rate (Zhang, 2002). Chapter 6 Natural Convection 579 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press relationship with Ra* is shown in eq. (6.264). Substituting eq. (6.264) into eq. (6.272), the following parameter is obtained; it is a parameter that accurately reflects the characteristics of the natural convection: 2 (θ w − θ c ) Nu (6.273) = 1/ 4 1/ 4 Ra Λθ w θ w It can be seen from Fig. 6.33 that the difference between the parameters defined by eqs. (6.272) and (6.273) vanishes and they become a constant equal to 0.374 when Bi/Ra*1/4>102, which is the same as the case of constant temperature heating in Bejan (1989). Therefore, the melting problem can be treated as boundary condition of the first kind when Bi/Ra*1/4 > 102. Nu/Ra1/4 increases with decreasing Bi/Ra*1/4 and its value approaches the results of a constant heat flux heating when Bi/Ra*1/4 = 0.1. The melting rate can be obtained by R= 1 V / V0 L =  U 0 (Y )dY 0 SteFo H (6.274) where V/V0 is the ratio of the liquid volume over the total volume of the enclosure. As can be seen from Fig. 6.33, the rate of melting is the same as that of boundary condition of the first kind when Bi/Ra*1/4 > 102. The rate of melting decreases with decreasing Bi/Ra*1/4 since the wall temperature decreases with decreasing Bi/Ra*1/4. 6.7 Instability Analysis of Natural Convection Natural convection in an enclosure has been a classic problem in fluid mechanics and heat transfer due to its fundamental and practical significance. It has attracted extensive experimental, analytical and numerical studies. The majority of the earlier studies have focused on the steady-state solutions, particularly for natural convection in a rectangular cavity with differentially heated sidewalls. In this case, the basic flow is a simple, unicellular circulation. When the horizontal temperature gradient increases, however, the base flow undergoes instability, and eventually becomes unsteady and turbulent. Computation of flow instabilities is one of the most fascinating and challenging aspects of the present computational fluid dynamics. The aim of instability analysis is to predict the values of the parameters at which instabilities occur, the nature of the instability (steady or unsteady, supercritical or subcritical), the type of solution that results from the instability, the modification of the base flow solution that results from the nonlinear interactions of the unstable modes and the subsequent route to chaos (Gadoin, et al., 2001). The instability analysis of natural convection in a vertical fluid layer (Bergholz, 1978), for example, shows that if the value of the Prandtl number is in the low to moderate range (Pr < 12.7), there is a transition from stationary to travelling-wave instability if the vertical temperature stratification exceeds a 580 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press certain magnitude. However, if the Prandtl number is large, the transition with increasing stratification, is from travelling-wave to stationary instability. These predictions of the instability analysis are in good agreement with experimental observations for both stationary and travelling-wave instabilities (e.g., Elder, 1965). Xin and Le Quere (2001) developed a general methodology for investigating the stability of a two-dimensional base solution with respect to both two- and three-dimensional perturbations by combining the Arnoldi’s method, the preconditioned Newton’s iteration and the preconditioned continuation method. They applied this methodology to investigate in detail the stability of natural convection in a differentially heated square cavity with conducting horizontal walls: 2D base solutions with respect to both 2D and 3D perturbations for a large range of Prandtl number. Their results show that for 2D perturbations, the critical Rayleigh number corresponding to the onset of unsteadiness increases with the Prandtl number. Table 6.4 lists the critical Rayleigh numbers and angular frequencies under 2D perturbations for Pr = 7 and 0.015, respectively. The corresponding base solutions and unstable modes are shown in Figs. 6.34 and 6.35, respectively. Figure 6.34 Base solution at Ra=5.5×106 (top) and eigenfunctions (temperature) of the first three unstable 2D modes (bottom) for Pr=7 (Xin and Le Quere, 2001). Figure 6.35 Base solution at Ra=4.5×104 (left and middle) and eigenfunction (temperature) of the unstable 2D mode (right) for Pr=0.015 (Xin and Le Quere, 2001). Chapter 6 Natural Convection 581 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Figure 6.36 Neutral curves (left) and angular frequencies (right) of 3D modes for Pr=0.0015 (Xin and Le Quere, 2001). Table 6.4 Critical Rayleigh numbers and angular frequencies under 2D perturbations Mode 1 Mode 2 Mode 3 Pr Rac ωc Rac ωc Rac ωc 7 5,573,177±1 2.466713 5,694,340±1 2.681573 6,305,615±1 2.883617 0.015 40,524±1 0.712348 – – – – The numerical results for 3D perturbations show that (1) for Pr=1, the base solutions are more unstable to 3D perturbations, and the most unstable 3D mode is connected to the third most unstable 2D mode; (2) for larger Prandtl numbers (3, 7 and 20), the most dangerous perturbations are two-dimensional; (3) for smaller Prandtl numbers (0.71, 0.1, and 0.015) the base solutions are also more unstable to 3D perturbations, but the most unstable 3D mode corresponds to a stationary mode, instead of connecting to the most unstable 2D modes. Moreover, the critical Rayleigh numbers for Pr=0.015 and 0.1 are approximately 40 times smaller than those found for 2D perturbations. For Pr=0.015, the neutral curves and angular frequencies of the 3D modes are shown in Fig. 6.36. As far as a differentially heated square cavity with conducting horizontal walls is concerned, the study by Xin and Le Quere (2001) suggests that for large Prandtl number two-dimensional flows, two-dimensional studies make sense at least up to the onset of unsteadiness, while for other Prandtl numbers, it is necessary to perform three-dimensional investigations. The need to do so is even stronger for small Prandtl numbers because the two-dimensional base flows are unstable to three-dimensional perturbations at very low Rayleigh numbers. 582 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press References Alamgir, M., 1979, “Overall Heat Transfer from Vertical Cones in Laminar Free Convection: An Approximate Method,” ASME Journal of Heat Transfer, Vol. 101, pp. 174-176. Arnold, J.N., Catton, I., and Edwards, D.K., 1976, “Experimental Investigation of Natural Convection in Inclined Rectangular Regions of Differing Aspect Ratios,” ASME J. Heat Transfer, Vol. 98, pp. 67-71. Bejan, A., 1989, “Analysis of Melting by Natural Convection in an Enclosure,” Int. J. Heat and Fluid Flow, Vol. 10, pp. 245-252. Bejan, A., 2004, Convection Heat Transfer, 3rd ed., John Wiley & Sons, New York. Bejan, A, and Lage, J.L., 1990, “Prandtl Number Effect on the Transition in Natural Convection along a Vertical Surface,” ASME Journal of Heat Transfer, Vol. 112, pp. 787 – 790. Bejan, A., and Tien, C.L., 1978, “Laminar Natural Convection Heat Transfer in a Horizontal Cavity with Different End Temperatures,” ASME J. Heat Transfer, Vol. 100, pp. 641-647. Bergholz, R.F., 1978, “Instability of Steady Natural Convection in a Vertical Fluid Layer,” J. Fluid Mech., Vol. 84, pp. 743-768. Berkovsky, B.M. and Polevikov, V.K., 1977, “Numerical Study of Problems on High-Intensive Free Convection,” Heat Transfer and Turbulent Buoyant Convection, Eds. Spalding, D.B., and Afgan, N., Vol. II, pp. 443-455, Hemisphere Publishing Corporation, Washington, DC. Cao, Y., and Faghri, A., 1990, “A Numerical Analysis of Phase Change Problem Including Natural Convection,” ASME Journal of Heat Transfer, Vol. 112, pp. 812-815. Caton, I., 1978, “Natural Convection in Enclosures,” Proceedings of the 6th International Heat Transfer Conference, Vol. 6, pp. 13-43, Toronto, Canada. Churchill, S.W., 1983, “Free Convection around Immersed Bodies,” Heat Exchanger Design Handbook, Schlunder, E.U., ed., Hemisphere Publishing, New York, NY. Churchill, S. W., and Chu, H.H.S., 1975, “Correlating Equations for Laminar and Turbulent Free Convection from a Vertical Plate,” Int. J. Heat Mass Transfer, Vol. 18, pp. 1323-1329. Cormack, D. E., Leal, L. G., and Seinfeld, J. H., 1974, “Natural Convection in a Shallow Cavity with Differentially Heated End Walls: Part 2. Numerical Solutions,” J. Fluid Mech., Vol. 65, pp. 231-246. Chapter 6 Natural Convection 583 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Date, A.W., 1986, “Numerical Prediction of Natural Convection Heat Transfer in Horizontal Annulus,” International Journal of Heat and Mass Transfer, Vol. 29, pp. 1457-1464. Eckert, E.R.G., and Jackson, T.W., 1951, “Analysis of Turbulent Free Convection Boundary Layer on a Flat Plate,” NACA Report 1015. Ede, A.J., 1964, “Advances in Free Convection,” Advances in Heat Transfer, Vol. 4, pp. 1-64, Academic Press, New York. Elder, J.W., 1965, “Laminar Free Convection in a Vertical Slot,” J. Fluid Mech., Vol. 23, pp. 77-98. Gadgil, A., and Gobin, D., 1984, “Analysis of Two Dimensional Melting in Rectangular Enclosures in Presence of Convection,” AMSE J. Heat Transfer, Vol. 106, pp. 20-26. Gadoin, E., Le Quere, P., and Daube, O. 2001, “A General Methodology for Investigating Flow Instabilities in Complex Geometries: Application to Natural Convection in Enclosures,” Int. J. Numer. Meth. Fluids, Vol. 37, pp. 175-208. Gebhart, B., Jaluria, Y., Mahajan, R.L., and Sammakia, B., 1988, BuoyancyInduced Flows and Transport, Hemisphere Publishing Co., Washington, DC. , Gebhart, B., Pera, L., and Schorr, A.W., 1970, “Steady Laminar Natural Convection Plumes above a Horizontal Line Heat Source,” International Journal of Heat and Mass Transfer, Vol. 13, pp. 161-171 Globe, S., and Dropkin, D., 1959, “Natural Convection Heat Transfer in Liquids Confined by Two Horizontal Plates Heated from Below,” ASME J. Heat Transfer, Vol. 81, pp. 24-28. Ho, C.J. and Viskanta, R., 1984, “Heat Transfer During Melting from an Isothermal Vertical Wall,” ASME J. Heat Transfer, Vol. 106, pp. 12-19. Imberger, I., 1973, “Natural Convection in a Shallow Cavity with Differentially Heated End Walls, Part 3. Experimental Results,” J. Fluid Mech., Vol. 65, pp. 247-260. Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S. 2007, Fundamentals of Heat and Mass Transfer, 6th ed., John Wiley & Sons, New York. Jaluria, Y., 2003, “Natural convection,” in Heat Transfer Handbook, Eds. A. Bejan and A.D. Kraus, Wiley, NY, pp. 525-571. Kays, W.M., Crawford, M.E., and Weigand, B., 2005, Convection Heat Transfer, 4th ed., McGraw-Hill, New York, NY. Kuehn, T.H., and Goldstein, R.J., 1976, “An Experimental and Theoretical Study of Natural Convection in the Annulus between Horizontal Concentric Cylinders,” J. Fluid Mech., Vol. 74, pp. 695-719. 584 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Le Fevre, E.J., 1956, “Laminar Free Convection from a Vertical Plane Surface,” Proceedings of the 9th International Congress of Applied Mechanics, Brussels, Vol. 4, pp. 168-174. Le Fevre, E.J., and Ede, A.J., 1956, “Laminar Free Convection from the Outer Surface of a Vertical Cylinder,” Proceedings of the 9th International Congress of Applied Mechanics, Brussels, Vol. 4, pp. 175-183. MacGregor, R.K., and Emery, A.P., 1969, “Free Convection through Vertical Plane Layers - Moderate and High Prandtl Number Fluids,” ASME J. Heat Transfer, Vol. 91, pp. 391-403 McAdams, W.H., 1954, Heat Transmission, 3rd Ed. McGraw-Hill, New York, NY. Merk, H.J., and Prins, J.A. 1953, Thermal Convection Laminar Boundary Layer I, Applied Scientific Research, Vol. A4, pp. 11-24. Merk, H.J., and Prins, J.A. 1954a, Thermal Convection Laminar Boundary Layer II, Applied Scientific Research, Vol. A4, pp. 195-206. Merk, H.J., and Prins, J.A. 1954b, Thermal Convection Laminar Boundary Layer III, Applied Scientific Research, Vol. A4, pp. 207-221. Oosthuizen, P.H., and Naylor, D., 1999, Introduction to Convective Heat Transfer Analysis, WCB/McGraw-Hill, New York. Ostrach, S., 1953, An Analysis of Laminar Free Convection Flow and Heat Transfer about a Flat Plate Parallel to the Direction of the Generating Body Force, NASA Report, 1111, pp. 63–79. Patterson, J., and Imberger, J., 1980, “Unsteady Natural Convection in a Rectangular Cavity,” J. Fluid Mech., Vol. 100, pp. 65-86. Pera, L., and Gebhart, B., 1973, “Natural Convection Boundary Layer Flow over Horizontal and Slightly Inclined Surfaces,” Int. J. Heat Mass Transfer, Vol. 16, pp. 1131-1146. Quon, C., 1983, “Effects of Grid Distribution on the Computation of High Rayleigh Number Convection in a Differentially Heated Cavity,” Numerical properties and Methodologies in Heat Transfer, Eds., Shih, T.M., pp. 261-281, Hemisphere, Washington, DC. Raithby, G.D., and Hollands, K.G.T., 1975, “A General Method of Obtaining Approximate Solutions to Laminar and Turbulent Free Convection Problems,” Irvine, T.F., and Hartnett, J.P., Eds., Advances in Heat Transfer, Vol. 11, pp. 265-315, Academic Press, New York, NY. Rotem, Z., and Claassen, L., 1969, “Natural Convection above Unconfined Horizontal Surface,” J. Fluid Mech., Vol. 39, pp. 173-192. Chapter 6 Natural Convection 585 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Tichy, J., and Gadgil, A., 1982, “High Rayleigh Number Laminar Convection in Low Aspect Ratio Enclosures with Adiabatic Horizontal Walls and Differentially Heated Vertical Walls,” ASME J. Heat Transfer, Vol. 104, pp. 103-110. Vliet, G.C., 1969, “Natural Convection Local Heat Transfer on Constant-HeatFlux Inclined Surface,” ASME J. Heat Transfer, Vol. 91, pp. 511-516. Wang, Q.J., Wang, C., and Chen, Z.Q., 1991, “Heat Transfer during Superheated Solidification around Bare Tube and Finned Tubes,” Experimental Heat Transfer, Fluid Mechanics and Thermodynamics 1991, Elsevier Science Publishing Co., Inc, pp. 1171-1176. Wang, S., Faghri, A., and Bergman, T.L., 2010, “A Comprehensive Numerical Model for Melting with Natural Convection,” Int. J. Heat Mass Transfer, Vol. 53, pp. 1986-2000. Warner, C.Y., and Arpaci, V.S., 1968, “An Experimental Investigation of Turbulent Natural Convection in Air at Low Pressure along a Vertical Heated Flat Plate,” Int. J. Heat Mass Transfer, Vol. 11, pp. 397-406. Xin, S., and Le Quere, P., 2001, “Linear Stability Analyses of Natural Convection Flows in a Differentially Heated Square Cavity with Conducting Horizontal Walls,” Physics of Fluids, Vol. 13, pp. 2529-2542. Yang, K.T., 1987, “Natural Convection in Enclosures,” Handbook of Single Phase Convective Heat Transfer, Eds. Kakac, S., Shah, R., and Aung, W., Wiley & Sons, New York, NY. Zhang, Y., 2002, “Analysis of Melting in an Enclosure under Boundary Condition of the Third Kind,” Proceedings of the 12th International Heat Transfer Conference, Vol. 3, pp. 315-320, Grenoble, France, August 18-23, 2002. Zhang, Y., and Chen, Z.Q., 1994, “Effect of Wall Conduction on Melting in an Enclosure Heated by Constant Rate,” Int. J. Heat Mass Transfer, Vol. 37, pp. 340-343. Zhang, Y., Chen, Z.Q., and Faghri, A., 1997, “Heat Transfer During Solidification around a Horizontal Tube with Internal Convective Cooling,” ASME Journal of Solar Energy Engineering, Vol. 119, pp. 44-47. Zhang, Z., and Bejan, A., 1989, “Melting in an Enclosure Heated at Constant Rate,” Int. J. Heat Mass Transfer, Vol. 32, pp. 1063-1076. Zhong, Z.Y., Lloyd, J.R., Yang, K.T., 1983, “Variable-Property Natural Convection in Tilted Square Cavities,” Numerical Method in Thermal Problems, Vol. III, Eds. Lewis, R.W., Johnson J.A, and Smith W.R., pp. 968-979, Pineridge Press, Swansea, UK. 586 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Problems 6.1. For an ideal gas, show that the volumetric coefficient of thermal expansion defined in eq. (6.10) is 1/T. 6.2. The scale analysis presented in Section 6.3 is for the case where the wall temperature is constant. In many applications, the wall heat flux is constant ( q′′ = const ). For the case where Pr > 1, what are the scales of the thickness of the thermal boundary layer and local Nusselt number? 6.3. Estimate the scales of the thickness of the thermal boundary layer and the local Nusselt number for natural convection of low Prandtl number (Pr < 1) over a vertical plate subject to constant heat flux. 6.4. Simplify equation (6.83) for the cases of Pr → 0 and Pr → ∞ . Express the local Nusselt number in terms of the dimensionless variables defined in Example 6.1. 6.5. Write a computer program to solve the simplified ordinary differential equation from the previous problem and eq. (6.84). Compare your results with eq. (6.85). 6.6. The average Nusselt number for laminar natural convection and the local Nusselt number at the trailing edge are related by eq. (6.86). Prove this relationship. 6.7. The empirical correlation for laminar natural convection that covers a wide range of Prandtl numbers was given by eq. (6.82). Simplify this correlation for the cases that Pr → 0 and Pr → ∞ . Compare your results with eq. (6.85). 6.8. Assume the velocity profile in the velocity boundary layer is a third degree polynomial function: u = a0 + a1 y + a2 y 2 + a3 y3 and use appropriate boundary condition to determine the unspecified constants. Compare your result with eq. (6.89). 6.9. The temperature profile in the thermal boundary layer can be assumed to be a second degree polynomial function. Show that the final temperature profile is eq. (6.90) after all constants are determined. 6.10. Assume U = C1 x1/ 2 and δ = C 2 x1/ 4 in Example 6.2, and determine C1, C2 and q from eqs. (6.101), (6.102) and (6.103). Show that the local Nusselt number can be expressed as eq. (6.104). 6.11. Determine the local Nusselt number for natural convection over a vertical plate subject to constant heat flux via an integral solution. The velocity and temperature profiles given by eqs. (6.89) and (6.90) can be used. Compare your result with eq. (6.106). Chapter 6 Natural Convection 587 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 6.12. The empirical correlation for laminar natural convection over a vertical plate under constant heat flux is given by eq. (6.106). Simplify this correlation for the cases that Pr → 0 and Pr → ∞ . 6.13. Substitute eq. (6.118) into eqs. (6.116) and (6.117) and show that the values of m and n must be 0.5 and 0.7, respectively. 6.14. Determine U and δ for turbulent natural convection using eqs. (6.116) and (6.117). Show that the local Nusselt number can be calculated by eq. (6.119). 6.15. Show that the average Nusselt number and the local Nusselt number at the trailing edge for turbulent natural covection are related by eq. (6.120). 6.16. The average Nusselt number for natural convection of hot surface with constant heat facing down can be obtained from eq. (6.134), provided that the Rayleigh number is defined using the average surface temperature. Since the average surface temperature is unknown a priori, it is desirable to obtain the correlation in term of flux Rayleigh number Ra *L = gβ q′′L4 / (αν k ) . Rewrite eq. (6.134) in term of Ra*L. 6.17. Derive eqs. (6.141) and (6.142) by integrating eqs. (6.136) – (6.138) with respect to r in the interval of (R, R + δ). 6.18. For integral solution of natural convection over a thin vertical cylinder, derive eqs. (6.146) – (6.149) by substituting eq. (6.145) into eqs. (6.143) and (6.144). 6.19. Prove eqs. (6.150) and (6.151). 6.20. For natural convection over a vertical cone with negligible effect of curvature, show that the dimensionless characteristic velocity and the boundary layer thickness can be expressed by eqs. (6.166) and (6.167). 6.21. When the effect of curvature on natural convection over a vertical cone is not negligible, the dimensionless characteristic velocity, U1, and boundary layer thickness, Δ, satisfy eqs. (6.169) and (6.170). Assume that U1 and Δ are power functions of X and prove eqs. (6.172) and (6.173). 6.22. Show that if the centerline temperature for natural convection over a line heat source satisfies eq. (6.182), the dimensionless stream function and dimensionless temperature satisfy eqs. (6.187) and (6.188). 6.23. Repeat the previous problem for natural convection over a point heat source and show that the dimensionless stream function and dimensionless temperature satisfy eqs. (6.197) and (6.198). 6.24. Show that the scales of the inertia, viscosity and buoyancy terms in eq. (6.210) satisfy eq. (6.211). 588 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 6.25. Perform a scale analysis for natural convection in a shallow enclosure ( H  L ) and estimate the scale of the horizontal and vertical velocity components, u and v. 6.26. A vertical double-paned window has a dimension of 1 m by 1 m and air is filled in the gap between the two layers of glasses that are separated by a distance of 25 mm. The air temperature inside the room is 20 °C and the temperature in the outside air is 34 °C. Assume that the radiation between the two layers of glasses, as well as natural convection between each pane and its surroundings can be neglected. What is the heat transfer rate through this window from the outside to inside of the room? 6.27. An annular space filled with air is formed by two concentric cylinders with length of 0.2 m and respective radii of 50 cm and 60 cm. The temperatures of the inner and outer cylinders are kept 80 °C and 24 °C, respectively. Find the rate of natural convection heat transfer from the inner cylinder to the outer cylinder. 6.28. Repeat the above problem for the annular space formed by two concentric spheres. 6.29. Integrate eqs. (6.243) and (6.244) in the interval of (0, δ) and use eqs. (6.245) and (6.246) to show eqs. (6.247) and (6.248). 6.30. Integrate eqs. (6.253) and (6.254) within the cold boundary layer and assume second-degree polynomial velocity and temperature profiles in the cold boundary layer to show eq. (6.255). Chapter 6 Natural Convection 589 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press