Back to Advanced Heat and Mass Transfer Home
Table of Contents
4 EXTERNAL CONVECTIVE HEAT AND MASS TRANSFER 4.1 Introduction Convective heat and mass transfer is the process of heat and mass transport through fluid flow. Understanding the fundamentals of fluid mechanisms is an important part of predicting and modeling convective transfer problems. Convective heat and mass transport can be generally classified as internal or external. In internal cases, the flow is constrained on all sides by solid boundaries except where there are inlets or outlets, such as flow in ducts and channels. In external cases, the fluid has at least one side that is not bounded by a solid surface. Flow over a flat plate is a classic example of external flow. In external flow, the temperature of the fluid far away from the surface as well as a condition on the surface such as temperature, heat flux, or a combination are usually known. In internal flow cases such information at the inlet and on the wall bounded by the fluid are normally known. Convective transport phenomena can also be classified as forced or natural. Forced convection transport is characterized by a convective force that originates in an external source, such as a fan or pump. Natural convection transport is characterized by fluid motion that is caused primarily by body forces resulting from properties such as density and mass concentration variation. In some cases, both forced and natural convection may exist together, however, in most cases one or the other is the dominant mechanism in the transport process. Note that transport phenomena in porous media is considered forced convection, even though it is created naturally by capillary forces. So far, we have briefly discussed convective heat transfer primarily as a boundary condition, in terms of the heat transfer coefficient, which is a function of flow, fluid properties, and geometry of flow. Convective heat and mass transfer describe momentum, energy, and/or mass transfer between a surface and a fluid in contact with that surface. Convective heat and mass transfer is caused by bulk fluid motion (advection) and the random motion of fluid molecules (conduction and diffusion). The goal of this chapter is to understand both the physical phenomena and mechanics, as well as methods to calculate the heat and mass transfer coefficients. Chapter 4 External Convective Heat and Mass Transfer 339 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u ω1,∞ T∞ U∞ T ω1 y, v vw x, u Tw ω1,w Figure 4.1 Convective transport of mass, momentum, and heat between a surface and a fluid. As an engineer designing a thermal fluid system consisting of fluid and solid components, the most basic and fundamental questions are: What is the force exerted by the fluid onto the component/system? What is the heat and mass transferred between the object and the fluid? For simplicity, consider that the component in question is part of a heat exchanger surrounded by a stream of fluid as shown in Figure 4.1. There is the possibility of both heat and/or mass transfer to and from the object and the fluid stream. Consider the convective mass, momentum, and heat transfer between a surface and a fluid as shown in Fig. 4.1, where a stream of fluid with uniform velocity U∞, temperature T∞, and mass fraction ω1,∞ flows over the surface. For simplicity, assume the mass transfer is binary and there is a fixed proportion of species 1 and 2 so that the mass fraction in the free stream ω1,∞ is constant, similar to T∞ and U∞. We further assume that there is a vertical velocity at the surface due to either mass diffusion of species 1 to the main stream, and/or forced convection (injection or suction) independent of mass diffusion. The net force exerted by the fluid on the plate is obtained by integrating the shear stress, τw, at the surface over the area of the surface F =  τ w dA (4.1) A where τ w = −μ ∂u ∂y (4.2) y=0 The mass flux of species 1 to the main stream due to convection and  diffusion is m1'' . ∂ω  = ρ1vw + ρ hm (ω1, w − ω1,∞ ) m1′′ = ρ1vw − ρ D12 1 (4.3) ∂y y =0 It is customary for the mass transfer coefficient, hm, to be defined based on diffusive flux; therefore it is important that the effect of induced convection (ρ1vw) be added separately when calculating the total mass flux. The total mass transfer rate between the surface and the fluid is obtained by integrating the mass flux over the surface area: 340 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press   m1 =  m1′′dA A (4.4) The heat transfer rate, q , between the fluid and flat plate is obtained by integrating the heat flux over the surface area ′′ q =  qw dA (4.5) A where ′′ qw = −k ∂T ∂y = h ( Tw − T∞ ) y =0 if radiative heat flux and axial heat conduction in the wall are ignored. Furthermore, the convective heat transfer due to injection or suction is neglected. The goal of this chapter is to present fundamental models, analytical and numerical solutions of both laminar and turbulent external forced convections. Sections 4.2 – 4.4 introduce the concept, approximations, and governing equations of the boundary layer theory, followed by discussions on similarity solution and integral solutions of laminar boundary layer problems in Sections 4.5 – 4.7. The numerical solution of forced convection problem based on the full Navier-Stokes equation using finite volume method is discussed in detail in Section 4.8.; it is then followed, in Section 4.9, by application of computational method, and analogies and differences in various transport phenomena processes. This chapter is closed by a very detailed discussion on the turbulence for external forced convection in Section 4.11. Internal forced convection and natural convection (for both internal and external) are presented in Chapters 5 and 6, respectively. 4.2 Concepts of the Boundary Layer Theory Unless the geometry is very simple, it is difficult to solve for the complete viscous fluid flow around a body. A full domain numerical solution is time consuming and impractical, because one needs to solve the full Navier-Stokes equations in the full domain, which are nonlinear, elliptic, and complex. In 1904, Prandtl discovered that for most practical applications, the influence of viscosity is observed in a very thin domain close to the object, as shown in Fig. 4.2. Outside this region one can assume the flow is inviscid (µ = 0) (Prandtl, 1904; Schlichting and Gersteu, 2000; White, 2005). The thin region where the effect of viscosity is dominant is called the momentum or viscous boundary layer. The solution of boundary layer analysis can be simplified due to the fact that its thickness is much smaller than the characteristic dimension of the object. The fluid adjacent to the surface of the body has zero relative velocity; ufluid – usurface = 0. This is also called the no slip boundary condition. One of the assumptions of boundary layer approximations is that the fluid velocity at the wall is at rest relative to the surface. This is true except when the fluid pressure is very low, and therefore, the Kundsen number, Kn= λ/L, of the Chapter 4 External Convective Heat and Mass Transfer 341 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U∞ U∞ Potential Flow Domain δ y, v x, u Viscous Boundary Layer Domain Figure 4.2 Viscous or momentum boundary layer. fluid molecules is much larger than 1. For external flow, the flow next to an object can be divided into two parts. The larger part is related to a free stream of fluid, in which the effect of viscosity is negligible (potential flow theory). The smaller region is a thin layer next to the surface of the body, in which the effects of molecular transport (such as viscosity, thermal conductivity and mass diffusivity) are very important. Potential flow theory neglects the effect of viscosity, and therefore, significantly simplifies the Navier-Stokes equations, which provides the solution for the velocity distribution. A disadvantage of the potential theory is that since second order terms are neglected, the effect of viscosity, no slip, and impermeability boundary conditions at the surface cannot be accounted for. In general, potential flow theory predicts the free stream field accurately, despite its simplicity. Boundary layer thickness is defined as the distance within the fluid in which most of the velocity change occurs. This thickness is usually defined as the thickness in which the velocity reaches 99% of the free stream velocity, u = 0.99U∞. Figure 4.3 illustrates how the momentum boundary layer thickness changes along the plate. Flow is laminar at relatively small values of x, where it is shown Turbulent region Laminar region U∞ y x Transition region Laminar sublayer δ Figure 4.3 Laminar and turbulent boundary layer flow over a flat plate. 342 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press as the laminar boundary layer region. As x increases, the fluid motion begins to fluctuate. This is called the transition region. The boundary layer may be either laminar or turbulent in this region. Finally, above a certain value of x, the flow will always be turbulent. There is a very thin region next to the wall in the turbulent region where the flow is still laminar, called the laminar sublayer. For flow over a flat plate, experimental data indicates the flow is laminar Rex ≤ 2×105 the flow is in transition 2×105 <Rex <3×106 the flow is turbulent Rex ≥ 3×106 where ρU ∞ x Re x = μ 4.3 Boundary Layer Approximation If one assumes that the boundary layer thickness, δ, is very small compared to the characteristic dimension of the object, one can make the assumption that δ is significantly less than L (δ  L), where L is the characteristic dimension of the object. Using scale analysis discussed in Chapter 1, for flow over a flat plate with constant free stream velocity, one can show, in a steady two-dimensional laminar boundary layer representation, that the following conditions are met within the boundary layer region, assuming there are no body forces: ∂ 2u ∂ 2u 2 ∂y 2 ∂x uv ∂p dp ≈ ∂x dx ∂p ≈0 ∂y (4.6) (4.7) (4.8) (4.9) A similar concept exists when there is heat and/or mass transfer between a fluid and the surface of an object. Again, the region in which the effect of temperature or concentration is dominant is, in the most practical case, a region very close to the surface, as shown in Fig. 4.4. In general, for the case of flow over a flat plate, one can expect three different boundary layer regions with thicknesses δ, δT, and δC, corresponding to momentum, thermal, and concentration boundary layers, respectively. δ, δT, and δC are not necessarily the same thickness and their values, as will be shown later, depend on the properties of the fluid such as kinematic viscosity (ν ), thermal diffusivity (α = k/ρcp), specific heat, c, and mass diffusion coefficient, D. Chapter 4 External Convective Heat and Mass Transfer 343 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U∞ T∞ ω∞ U∞ y δ x u=0 (a) velocity profile U∞ T∞ ω∞ T∞-Tw y δT x T = Tw (b) temperature profile ω∞-ωw U∞ T∞ ω∞ y δC x ω=ωw (c) mass fraction profile Figure 4.4 Boundary layer concept over a flat plate. Using boundary layer approximation, scale analysis, and order of magnitude, one can show that similar approximations exist for thermal and mass concentration boundary layer analysis. ∂T ∂T  , ∂y ∂x ∂ 2T ∂ 2T 2 ∂y 2 ∂x (4.10) (4.11) ∂ω ∂ω  , ∂y ∂x ∂ 2ω ∂ 2ω 2 ∂y 2 ∂x where ω is the mass fraction. 4.4 Governing Equations for Boundary Layer Approximation As noted before, one can obtain all pertinent information related to momentum, heat, and mass transfer between a surface and a fluid by focusing on the thin region (boundary layer) next to the surface, and solving the governing equations 344 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press including the boundary conditions. This provides significant simplification irrespective of whether an analytical or numerical approach to solve the physical problem is used. For most practical applications, the effect of molecular transport due to mass, momentum, energy, and species is dominant in this thin region. The purpose of this section is to develop the transport phenomena equations for boundary layer approximation, starting from the original differential conservation equations developed in Chapter 2. Consider flow over a flat plate as shown in Figure 4.4 with constant free stream velocity, temperature, and mass fraction of U∞, T∞, and ω∞, respectively. The surface wall is kept at constant temperature and concentration. Starting from the differential conservation equations for mass, momentum, energy, and species (as presented in Chapter 2), and assuming two-dimensional, steady, laminar flow, and constant properties in a Cartesian coordinate system, the following results are obtained: Continuity equation ∂u ∂v + =0 ∂x ∂y (4.12) Momentum equation in the x-direction u  ∂ 2u ∂ 2u  ∂u ∂u 1 ∂p +v =− +ν  2 + 2  ρ ∂x ∂x ∂y ∂y   ∂x  ∂2v ∂2v  ∂v ∂v 1 ∂p +v =− +ν  2 + 2  ρ ∂y ∂x ∂y ∂y   ∂x (4.13) Momentum equation in the y-direction u (4.14) 2 Energy equation u  ∂ 2T ∂ 2T ∂T ∂T +v =α 2 + 2 ∂x ∂y ∂y  ∂x u +v  ν  ∂u  +    cp  ∂y  (4.15) Species equation for a binary system  ∂ 2ω1 ∂ω1 ∂ω1 ∂x ∂ 2ω1  = D12  2 +  ∂y ∂y 2   ∂x (4.16) Assuming there is a simultaneous transfer of momentum, heat and mass, the following boundary conditions are one possibility for conventional applications: u ( x , 0 ) = 0 No slip boundary condition (4.17) vw = 0 Impermeable wall v( x , 0) =  vw > 0 Injection and vw < 0 suction u ( x , ∞ ) = U ∞ , U ∞ is the free stream velocity Tw = uniform temperature at the wall  heat transfer=  ∂T ′′ = qw = constant −k ∂y y=0  T ( x , ∞ ) = T∞ , T∞ is the free stream temperature (4.18) (4.19) (4.20) (4.21) Chapter 4 External Convective Heat and Mass Transfer 345 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ω1 ( x , 0 ) = ω1, w , ω1, w is the constant mass fraction at the wall ω1 ( x , ∞ ) = ω1,∞ , ω1,∞ is the free stream mass fraction (4.22) (4.23) At the wall, the temperature or heat flux (or a combination of the temperature and heat flux variation) which may change along the wall, should be known. Equation (4.20) presents the two limiting cases of constant wall temperature or constant heat flux at the wall. The normal velocity at the wall is zero for the case of no mass transfer from the wall; however, there are three classes of problems in which vw ≠ 0 at the wall. 1. Mass transfer due to phase change, such as condensation or evaporation. In such cases, temperature, mass fraction, and vw are coupled at the wall through mass and energy balance at the surface of the wall. 2. Injection or blowing of the same fluid through a porous wall (for example to protect the surface from a very high temperature main stream). In these cases, vw is positive, known, and independent of temperature. 3. Suction of the same fluid through a porous wall (for example to prevent boundary layer separation because of an adverse pressure gradient). In these cases, vw is also known and independent of temperature, but it is negative. In general, not all mechanisms of transport phenomena (mass, momentum, and energy) occur simultaneously, even though it is possible in many practical applications. If the transport phenomena is due to momentum and heat and not mass diffusion, then eqs. (4.12) to (4.15) are uncoupled for the case of constant properties. If there is mass transfer by diffusion from the wall however, then vw = f (ω1 ) ≠ 0 and eqs. (4.12) through (4.16) are coupled, even if the properties are constant. The occurrence of mass transfer at the wall may cause vw to be nonzero and can be calculated from Fick’s Law for binary mass diffusion. The coupling effect between mass and momentum in binary mass diffusion can be easily observed by calculating mass flux at the wall as shown below:  m1′′ = ρ1, w vw − ρ D12 ∂ω1 ∂y  ∂ω = ρ  ω1, w vw − D12 1  ∂y y=0     y =0  (4.24) where ω1,w and vw are mass fraction and the y-component of the velocity at the wall, respectively. There are five dependent variables (u, v, T, p, ω1) and two independent variables (x and y). There are five equations (4.12)-(4.16) and five unknowns with appropriate boundary conditions; therefore, it is a well posed problem. Using scale analysis, presented in Chapter 1, order of magnitude analysis, and considering boundary layer assumptions, as noted in Section 4.3, eqs. (4.12) through (4.16) are reduced to the following forms: Continuity equation: ∂u ∂v + =0 ∂x ∂y (4.25) 346 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Moment equation in the x-direction: u ∂u ∂u ∂2u +v =ν 2 ∂x ∂y ∂y (4.26) Energy equation: u ∂T ∂T ∂ 2 T ν ∂ 2T +v =α 2 = ∂x ∂y Pr ∂y 2 ∂y +v ∂ω1 ∂ 2ω1 ν ∂ 2ω1 = D12 = Sc ∂y 2 ∂y ∂y 2 (4.27) Species equation: ∂ω1 u ∂x (4.28) Using boundary layer approximation for flow over a flat plate with constant free stream velocity, the momentum equation in the y-direction and the pressure gradient in the x-direction were eliminated. Eqs. (4.25) through (4.28) are now transformed into a parabolic form rather than the original elliptic form (Chapter 2). There are now four dependent variables (u, v, T, ω1) and four equations. This is much easier to solve analytically or numerically because the axial diffusion terms (∂ 2 u / ∂x 2 , ∂ 2T / ∂x 2 and ∂ 2ω1 / ∂x 2 ) disappear. Example 4.1: Using boundary layer approximation and order of magnitude analysis, reduce the continuity and momentum equations (4.12), (4.13) and (4.14) to eqs. (4.25) and (4.26). Equations (4.25) and (4.26) correspond to two-dimensional, steady, laminar flow with constant free stream velocity and constant properties for flow over a flat plate. Solution: The characteristic dimension is assumed to be the length of the plate and δ is the order of magnitude of the distance in which u changes from zero at the wall to approximately U∞ at the free stream. The continuity and momentum equations are made dimensionless by using the following dimensionless variables: x∗ = x y u v , y∗ = , u∗ = , v∗ = , L L U∞ U∞ δ∗ = δ L << 1, P ∗ = LU ∞ P , Re = ρU ∞ 2 ν (4.29) The dimensionless form of the continuity and momentum equations in the x and y directions are: Continuity: ∂u∗ ∂v∗ + =0 ∂x ∗ ∂y∗ (4.30) Momentum (x-direction): u∗ ∂u∗ ∂u∗ ∂P ∗ 1  ∂ 2 u∗ ∂ 2 u∗  + v∗ ∗ = − ∗ + +   ∂x ∗ ∂y ∂x Re  ∂x ∗2 ∂y∗2  (4.31) Chapter 4 External Convective Heat and Mass Transfer 347 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Momentum (y-direction): u∗ ∂v∗ ∂v∗ ∂P ∗ 1  ∂ 2 v∗ ∂ 2 v∗  + v∗ ∗ = − ∗ + +   ∂x ∗ ∂y ∂y Re  ∂x ∗2 ∂y∗2  (4.32) To determine the significance of each term, the above dimensionless equations are subjected to an order of magnitude analysis. The symbol “O” is used to indicate the “order of” various terms in the above equations. The order of magnitudes of various dimensionless variables are: u∗  O (1), x ∗  O(1) and y∗  O (δ ∗ ) where δ ∗  1 (4.33) Continuity equation (4.30) in the form of scaling is shown below: O (1) O (1) + O (δ ∗ ) v∗ ~0 (4.34) From above * v∗ ~ 0 (δ ∗ ) << 1 (4.35) Clearly v is very small but nonzero, otherwise the boundary layer concept does not hold. The momentum equation in the x-direction, eq. (4.31), in the form of an order of magnitude analysis is shown below: O (1) O (1) O (1) + O (δ ∗ ) O (δ ∗ ) O (1) ~  O (1) 1  O (1) +  O (12 ) O δ ∗ 2 Re   ()  ∗  − ∂p  ∂x ∗   (4.36) It is easy to see that the first term in the brackets on the right hand side of the above equation is much smaller than the second term in brackets, and therefore it can be neglected, i.e., ∂2u ∂2u << 2 ∂x 2 ∂y (4.37) It is also evident that Re ~ O δ∗ () 2 1 (4.38) Similarly, the y-direction momentum eq. (4.32) is presented in the form of the order of magnitude analysis. O (1) O (δ ∗ ) O (1) + O (δ ∗ ) O (δ ∗ ) O (δ ∗ ) ~ O δ∗ ( )  O (1)  2 O δ ∗ ()   + O (δ ∗ )  ∂p∗ − 2 ∗ O δ ∗  ∂y   () (4.39) It can be easily concluded that ∂p∗ ∂p  O (δ * )  =0 ∗ ∂y ∂y (4.40) For the case of flow over a flat plate, U∞ is constant, and upon applying the x-direction momentum equation in the free stream region (µ = 0), we can conclude 348 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press dp =0 dx (4.41) As discussed before, in general, the free stream velocity outside the boundary layer, U∞, is obtained from the inviscid momentum equation outside the boundary (potential flow theory) ρ ∂U ∞  ∂U ∞ +U∞ ∂x  ∂t ∂p   = − ∂x  (4.42) where the free stream velocity is either known or determined from potential flow theory. One can assume ∂p / ∂x is known in the boundary layer equation. For flow over a flat plate, U∞ = constant. The pressure gradient, dp / dx , is zero for laminar flow over a flat plate. With constant free stream velocity, however, in most flow situations, a pressure gradient does exist. Pressure gradient plays a major role in momentum, heat, and mass transfer and it is important to provide some basic understanding of the role of dp / dx in fluid flow. The momentum boundary layer equation at the wall, y = 0, where u = v = 0, becomes μ ∂2u ∂y 2 = y=0 dp dx (4.43) which relates the velocity gradient at the surface to the pressure gradient. Figure 4.5 shows the variation of u, ∂u / ∂x , and ∂ 2 u / ∂y 2 across the boundary layer for dp / dx = 0 , dp / dx < 0 and dp / dx > 0 . When dp / dx = 0 , ∂ 2 u / ∂y 2 = 0 and therefore the velocity profile is linear near the wall. ∂u / ∂y is maximum at the wall and decreases to zero near the free surface. This decrease in the velocity gradient, y dp >0 ∂x y y dp <0 ∂x dp <0 ∂x dp =0 ∂x dp <0 ∂x u dp >0 ∂x dp =0 ∂x dp =0 ∂x ∂u ∂y dp >0 ∂x ∂ 2u ∂y 2 Figure 4.5 The variation in velocity and its first and second order derivatives across the momentum boundary layer for various pressure gradients. Chapter 4 External Convective Heat and Mass Transfer 349 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U∞ δ Separation point dp / dx > 0 . Separation region Figure 4.6 Velocity profile for separated flow when ∂u / ∂y , means ∂ 2 u / ∂y 2 < 0 . The second derivative of velocity is zero at the wall, negative within the boundary layer, and approaches zero at the outer edge of the negative side. The profiles of velocity variation for dp / dx < 0 are very similar to the case of zero pressure gradient, except that they produce a negative ∂ 2 u / ∂y 2 at the wall. The value of ∂ 2 u / ∂y 2 at the wall is positive when dp / dx > 0 ; however, it is zero at some point within the boundary layer. Within the boundary layer, a zero second derivative velocity is associated with an inflection point, which is related to flow separation as shown in Figure 4.6. Separation occurs when the velocity gradient of the fluid adjacent to the wall is either zero or negative. Positive pressure gradient is necessary for flow separation, but is not a sufficient condition. A negative pressure gradient is usually called a favorable pressure gradient compared to when the pressure gradient is positive and is called an adverse pressure gradient. 4.5 Laminar Boundary Layer Solutions for Momentum, Heat, and Mass Transfer There are three basic approaches to solve boundary layer equations for momentum, heat, and mass transfer: 1. Similarity solutions 2. Integral methods 3. Full numerical solution First we will consider the similarity approach, since it was the original classic method developed to solve boundary layer problems analytically. In this approach, using the fact that in some circumstances velocity is geometrically similar along the flow direction, the conservation partial differential equations are converted to ordinary differential equations. Similarity methods historically 350 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press provided significant insight and information about the physical boundary layer phenomenon when computer and numerical methodologies for the solution of partial differential equations were non-existent. However, there are limitations to the use of similarity methods in terms of applications. In general, they are only applicable to two-dimensional laminar flow for conventional geometry and boundary conditions with constant properties. The integral methods will be presented second. Integral methods are approximate solutions and provide closed form solutions by assuming a profile for velocity, temperature, and concentrations. Finally, there is the full numerical solution. The full umerical solution, along with the availability of proven numerical schemes, practical simulation codes, and high speed computers, are more efficient than other methods for complex problems. 4.6 Similarity Solutions The form of the velocity, temperature, and concentration profiles for flow over a flat plate are presented in Figure 4.4, and are based on qualitative measurements. These presentations indicate the possibility that these profiles are geometrically similar to each other along the flow direction for dependent variables such as velocity, temperature, and concentration, if the coordinates are stretched properly. For example, the velocity profiles are geometrically similar along the flow in the x-direction, differing only by a stretching factor (similarity factor) if the coordinates are properly sketched. As noted before, not all boundary layer flow configurations have similar geometric profiles (similarity solution), but some do (Burmeister, 1993), especially for more conventional geometries and boundary conditions. The steady, two-dimensional, laminar, momentum boundary layer equation is a second order, nonlinear, partial differential equation. The two-dimensional, steady, laminar flow, energy, and species boundary layer equations are linear, second order, partial differential equations for most constant properties, and/or for decoupled heat and mass transfer problems. If a similarity solution exists for a given situation, a mathematical transformation of coordinate systems can be performed to reflect this fact. A similarity technique converts conservation partial differential equations into ordinary differential equations and therefore, makes the solution much simpler. However, the analytical solutions may still require numerical integrations. The conservation equations are uncoupled when each equation and its boundary condition can be solved independently of one another, except for continuity and momentum, which must be solved simultaneously. Coupled transport phenomena can occur in some applications because of coupled conservation governing equations and/or boundary conditions for mass diffusion, momentum, or heat transfer. Chapter 4 External Convective Heat and Mass Transfer 351 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 4.6.1 Uncoupled Mass, Momentum, and Heat Transfer Problems Similarity Solution for Flow over a Flat Plate We first apply the similarity solution to the classic problem of forced convective flow over a flat plate with constant free stream velocity, where mass, momentum, and heat transfer equations are uncoupled. Next, similarity solutions of eqs. (4.25) to (4.27) with appropriate boundary conditions corresponding to Fig. 4.4 (a) and (b) are presented below. The viscous dissipation term is neglected for simplicity even though it is possible to obtain a similarity solution with viscous dissipation for this simple configuration. Assuming geometrically similar velocity profiles are possible along the flow over a flat plate using the following non-dimensional variables (Kays et al., 2005, Bejan, 2004, Burmeister ,1993): u = f (η ) (4.44) η = y ⋅ g( x ) (4.45) where η is a non-dimensional independent variable (function of both x and y) that fulfills the similarity requirement. Mass and momentum boundary layer equations for flow over a flat plate [eqs. (4.25) and (4.26)] can be converted to a new coordinate system η , f, and g by using the notation f ′ = dg df d2 f , g′ = and f ′′ = 2 . dη dx dη ∂u ∂f ∂f ∂η = = = f′ g ∂y ∂y ∂η ∂y (4.46) (4.47) (4.48) (4.49) (4.50) ∂u ∂f ∂f ∂η yg = = = f′ ′ ∂x ∂x ∂η ∂x ∂2u ∂2 f ∂  ∂f  ∂  ∂f   ∂η  = = g2  =    = f ′′ ∂y 2 ∂y 2 ∂y  ∂y  ∂η  ∂y   ∂y  Substituting eqs. (4.46)-(4.48) into (4.25) and (4.26) yields: ν f ′′ 2 = ff ′ ′ + vf ′ g yg g ∂v f ′ ′+ yg =0 ∂y Combining eqs. (4.49) and (4.50) to eliminate v and separating variables gives: 1 d  f ′′  1 g′  = f dη  f ′  ν g3 (4.51) Considering that both η and x are independent variables, the left side of eq. (4.51) is a function of η , and the right side is a function of x, then each side must be a constant. 1 g′ ν g3 = −c1 (4.52) 352 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Integrating (4.52) gives: − 1 = −c1ν x + c2 2g2 (4.53) At the leading edge of the flat plate ( x = 0 ), the boundary condition is u = U ∞ which can be satisfied only at η → ∞ . In other words, η must be infinite when x = 0 . It follows from eq. (4.45) that g(0) must be infinite and consequently, c2 in eq. (4.53) must be zero. Therefore η= g= y 2c1ν x 1 2c1ν x (4.54) (4.55) Since c1 is constant, the axial velocity is of the following form: y u= f   x (4.56) From eq. (4.51) 1 d  f ′′    = −c1 f dη  f ′  (4.57) Integrating eq. (4.57) gives: f ′′ = −c1  fdη + c3 f′ (4.58) The following boundary conditions are used to evaluate constants c1 and c3: at y = 0, u = v = 0 or at η = 0, f = 0 (4.59) Using eqs. (4.58) and (4.59): f ′′ = 0 at η = 0 (4.60) then, using eqs. (4.58) and (4.60): c3 = 0 therefore, η f ′′ = −c1  fdη 0 f′ (4.61) Let the function ς ′ = dς / dη represent the non-dimensional axial velocity, u / U ∞ : ς′ = u f = U∞ U∞ (4.62) Then, using eq. (4.62) f ′′ ς ′′′ = f ′ ς ′′ (4.63) Substituting (4.63) in (4.61) gives: η ς ′′′ dς = −c1  U ∞ dη = −c1U ∞ς 0 ς ′′ dη (4.64) Chapter 4 External Convective Heat and Mass Transfer 353 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1.0 0.8 0.6 u =f U∞ 0.4 v Rex1/2 U∞ 0.2 0 1 2 3 η = Rex y x 4 1 2 5 6 7 Figure 4.7 Dimensionless velocity distribution in a laminar boundary layer. or ς ′′′ + c1U ∞ςς ′′ = 0 (4.65) From (4.65) we conclude that c1U∞ is a constant, non-dimensional variable. Let c1U ∞ = 1 / 2 , eq. (4.65) becomes ς ′′′ + ςς ′′ = 0 1 2 (4.66) u U∞ where η= y ν x / U∞ and ς ′(η ) = (4.67) Equation (4.66) is known as the Blasius equation and is a nonlinear, third order, ordinary differential equation, which requires three boundary conditions. ς ′(0) = 0, ς ′(∞) = 1, ς (0) = 0 (4.68) Equation (4.66) can be numerically integrated using boundary conditions given by eq. (4.68). Blasius clearly integrated eq. (4.66) by hand without using a computer (Blasius, 1908). The numerical results using the shooting method are presented in Fig. 4.7, which shows that both u and v as functions of the independent variable η . A number of conclusions can be made from the numerical results presented in Fig. 4.7. 1. The edge of the momentum boundary layer is approximately at η =5. 2. The vertical velocity component v is non-zero at the edge of the boundary layer at η =5, even though the free stream flow is parallel to the plate. 3. The momentum boundary layer thickness is approximately 354 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ x = 5 Re x U∞ x (4.69) where ν Using the definition of shear stress at the wall, τw, and the skin friction coefficient, c f , we obtain: τw = μ ∂u ∂y = cf y=0 2 ρU ∞ Re x = 2 (4.70) Rearranging the above equation in terms of ς gives: ς ′′ y = 0 cf ν = ς ′′ η = 0 = 2 U∞ x cf 2 Re x = 0.332 Re x (4.71) (4.72) For the purposes of design and comparison with experimental measurements, it is more practical to obtain the average value of the skin friction, cf, which is usually determined by 2 2 x ρU ∞ ρU ∞ dx = c f xW total drag =  τ w dA = W  c f (4.73) 0 A 2 2 where W is the width of the plate. From eq. (4.72), cf 2 = 1 x cf dx = 0.664 Re −1/ 2 x x 0 2 (4.74) It should be emphasized that the results presented above are for laminar flow parallel to a flat plate with constant free stream velocity. The approach to obtain the velocity distribution presented above can be extended to obtain the temperature distribution for laminar boundary layer flow over a flat plate with constant properties, constant free stream velocity, and constant temperature. The constant property assumption is important to decouple the momentum and energy equations in this configuration in order to obtain velocity first and temperature second. The non-dimensional temperature, θ , is defined as: θ= T − Tw T∞ − Tw (4.75) The non-dimensional form of the energy equation (4.27) (neglecting viscous dissipation) and boundary conditions and using the variables η and ς noted in eqs. (4.67) and (4.62) are as shown below. θ ′′ + Pr ςθ ′ = 0 2 θ (0) = 0 (4.76) (4.77) (4.78) θ (∞ ) = 1 Chapter 4 External Convective Heat and Mass Transfer 355 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press with the notation that θ ′ = dθ d 2θ and θ ′′ = 2 . dη dη Equation (4.76) can be directly integrated: dθ ′ Pr + ςθ ′ = 0 dη 2 dθ ′ Pr + ς dη = 0 θ′ 2 θ = c1  exp  −  ς dη   dη + c2 0  20   Applying boundary conditions (4.77) and (4.78) gives c2 = 0 and c1 = η (4.79) (4.80) (4.81)   Pr η   ∞ 0 1   Pr η  exp  − 2 0 ς dη   dη    (4.82) Therefore, we obtain the following equation for the dimensionless temperature, θ:   Pr η  exp  − 2 0 ς dη   dη   θ (η ) =  ∞ Pr η   0 exp  − 2 0 ς dη  dη     η 0 (4.83) The energy equation (4.76) is a second order, linear, ordinary differential equation. The dimensionless temperature, θ , is a function of η and Pr. Equation (4.83) can be integrated numerically for a given Pr. The numerical results of θ as a function of η for different Pr values are presented in Fig. 4.8. 1 0.75 0.5 Pr = 0.5 1 3 10 50 1000 = 0.25 0 0 1 2 = 3 Re 1 ⁄2 ∞ − − 4 5 6 Figure 4.8 Dimensionless temperature distribution in a laminar boundary layer over a flat plate. 356 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The following conclusions can be made upon reviewing the results of the momentum and thermal boundary layer equations, as presented in Figs. 4.7 and 4.8. 1. Temperature profile is identical to velocity for Pr = 1. 2. Temperature profile has strong dependence on the Prandtl number. 3. For Pr<1, the thermal boundary layer thickness is greater than the momentum boundary layer thickness, δ T > δ . 4. For Pr>1 the thermal boundary layer thickness is less than the momentum boundary layer thickness, δ T < δ . The local heat transfer coefficient can be easily obtained from the numerical results of eq. (4.83) for a given Prandtl number by using Fourier’s law of heat conduction. ∂T dθ ∂η −k −k ∂y y = 0 dη ∂y w ′′ qw (4.84) = = hx = Tw − T∞ Tw − T∞ Tw − T∞ The local Nusselt number is presented in terms of dimensionless variables as follows: x θ ′ η =0 hx = Re1/ 2 θ ′ ( 0 ) Nux = x = (4.85) x k ν x / U∞ or by using eq. (4.83): Nux = Re1/ 2 x ∞  Pr η  0 exp  − 2 0 ς dη   dη    (4.86) The numerical results can be approximated for the range of Prandtl numbers from 0.5 to 15 by the following equation: Nux = 0.332 Pr1/3 Re1/ 2 (4.87) x The ratio of thermal boundary layer thickness, δT, to momentum boundary layer thickness, δ, for flow over a flat plate can also be approximated over this range of Prandtl numbers according to the relationship below: δT = Pr −1/3 , 0.5 ≤ Pr < 10 (4.88) δ For both design purposes and comparison with experimental measurement it is more convenient to calculate the average or mean heat transfer coefficient, h . q =  h ( Tw − T∞ ) dA = W  h ( Tw − T∞ )dx = h W ( Tw − T∞ ) x A 0 x (4.89) From above, the mean heat transfer coefficient and mean Nusselt number become: h= 1x k hx dx = 0.664 Re x1/ 2 Pr1/3 x 0 x hx = 2 Nux = 0.664 Re x1/ 2 Pr1/3 k (4.90) (4.91) Nu = Chapter 4 External Convective Heat and Mass Transfer 357 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press This shows that the average heat transfer coefficient is twice the local heat transfer coefficient, just as the mean skin friction coefficient is twice the local skin friction coefficient for flow over a flat plate. It should be noted that these conclusions concerning the relation between the local and average transfer coefficients can not be extended to other configurations. Example 4.2: Obtain the local Nusselt number for flow over a flat plate with constant free stream velocity and temperature for the limiting case of very low Prandtl numbers. Solution: For very low Prandtl numbers, such as with liquid metals, the thermal boundary layer will develop much faster than the momentum boundary layer, as shown in Figure 4.9. δT  δ Also, as shown in Figure 4.9, one can safely assume the velocity inside the thermal boundary is constant and equal to U∞ for very low Prandtl numbers. δT U∞ y Velocity profile T∞ x δ Figure 4.9 Laminar momentum and thermal boundary layer for a fluid with a very low Prandtl number, Pr<<1, over a flat plate. Equation (4.76) can be simplified using ς ′ = u / U ∞ = 1 to obtain the following equation: d (θ ′′ / θ ′ ) Pr + =0 (4.92) dη 2 The above equation needs to be integrated three times with respect to η and therefore requires three boundary conditions. The additional boundary condition is θ ′′(0) = 0 , which can be obtained from the energy equation. The numerical result of eq. (4.92) yield: Nux = 0.565 Pr1/ 2 Re1/ 2 for Pr  1 (4.93) x Example 4.3: Obtain the local Nusselt number for flow over a flat plate when the Prandtl number is much greater than 1. 358 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press U(x) y x U∞ β Figure 4.10 Wedge flow configuration. Solution: When the Prandtl number is very high, the momentum boundary layer is much thicker than the thermal boundary layer; therefore, the thermal boundary layer is within the momentum boundary layer. We can assume a linear velocity near the wall within the thermal boundary layer. From the numerical results presented in Figure 4.7, one can conclude that ς ′′ ( 0 ) = 0.332 . Similar results are obtained through the basis of linear velocity. ς ′′ = 0.332 (4.94) Integrating the above equation twice with respect to η and applying ς ′(0) = ς (0) = 0 gives: ς = 0.166η 2 (4.95) Substituting eq. (4.95) into eq. (4.83) and repeating the numerical integration yields the Nusselt number from eq. (4.85) Nux = 0.339 Pr1/3 Re1/ 2 for Pr  1 (4.96) x Similarity Solution for Flow over a Wedge Similarity solutions also exist for some other conventional geometries with simple boundary conditions. Some of these cases are discussed below. For flow over a wedge with the oncoming velocity of U∞, as shown in Fig. 4.10, the free stream velocity varies according to potential (inviscid) flow theory: U = cx m (4.97) where U is the free stream velocity at the outer wedge surface and m is related to the wedge angle β by: β /π x dU m= = (4.98) 2 − ( β / π ) U dx  When β (measured in radians) is positive, the free stream velocity increases along the wedge surface. Chapter 4 External Convective Heat and Mass Transfer 359 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press    When β is negative, the free stream velocity decreases along the wedge surface. β = 0 corresponds to flow parallel to a flat plate. β = π corresponds to flow perpendicular to the walls. Equation (4.97) can be also applied near the leading edge of a blunt object. For example, the free stream velocity on the surface of a cylinder and a sphere near a stagnation point can be predicted by potential flow theory according to the equations below (Schlichting and Gersteu, 2000): (cylinder) (sphere) x x U = 2U ∞ sin   ≈ 2U ∞   R R 3 x 3 x  U = U ∞ sin   ≈ U ∞ R 2  R 2 (4.99) (4.100) where R is the radius of the cylinder or sphere, U∞ is the oncoming velocity and U is the free stream velocity just outside the boundary layer. The pressure gradient term is related to free stream velocity by potential flow theory. This can be done using the inviscid momentum equation ( μ = 0 ) 1 dp dU U 2m = −U = −cx m cmx m −1 = − dx x ρ∞ dx (4.101) The two dimensional, steady, constant property boundary layer equations for mass, momentum, and energy of flow over a wedge neglecting viscous dissipation are: Continuity equation: ∂u ∂v + =0 ∂x ∂y (4.102) Momentum (x-direction): u ∂u ∂u ∂ 2 u 1 dP +v =ν 2 − ρ dx ∂x ∂y ∂y (4.103) Energy equation: u ∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y (4.104) The above equations can be easily obtained by beginning with the general equations from Chapter 2 and making appropriate boundary layer approximations using scaling and order of magnitude analysis. Using eq. (4.101) for the pressure and comparable dimensionless similarity parameters to flow over a flat plate, we get the following ordinary differential equation for momentum and energy for constant wall temperature with no blowing or suction at the wall. ς ′′′ + 1 ( m + 1) ςς ′′ + m (1 − ξ ′2 ) = 0 2 (4.105) 360 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press θ ′′ + Pr ( m + 1) 2 ςθ ′ = 0 (4.106) The boundary conditions for velocity and temperature are the same as in the case of flow over a flat plate. Assuming that heating starts at the leading edge, the ordinary differential equations (4.105) and (4.106) for dimensionless velocity and temperature can be solved using the given boundary condition. The friction coefficient and Nusselt number can be calculated from the numerical results by the following equations: 2ς ′′(0) c f ,x = (4.107) 1/ 2 Re x Nux = A Re1/ 2 x (4.108) where the constant A depends on m and Pr. Solutions of eqs. (4.105) and (4.106) for m > 0 are unique. For m < 0 two groups of solutions exist for a limited range of m values corresponding to negative β . These two groups correspond to a decelerating mainstream to the point of incipient separation or laminar boundary layer flows after separation. Falkner and Skan (1931) originally developed the similarity solution for flow over a wedge with U∞ = cxm and an impermeable wall. The numerical results for the friction coefficient for selected values of m and β are presented in Table 4.1. Table 4.1 The local friction coefficient for laminar boundary layer flow over a wedge, with U = cxm and an impermeable wall. β m ∞ 1.0 0.333 0.111 0 -0.0654 -0.0826 -0.0904 ς ′′(0) = c f Re x1/ 2 1 2 2π π π/2 π /5 0 -0.14 -0.18 -0.199 ∞ 1.233 Stagnation 0.757 0.512 0.332 Flat plate 0.164 0.082 0 Separation Table 4.2 Local values of NuxRex-1/2 for laminar boundary layer flow over a wedge with constant wall and free stream temperature and an impermeable wall with U∞ = cxm β m -0.512 0 π/5 π/2 π 8π/5 -0.0753 0 0.111 0.333 1.0 4 0.7 0.242 0.292 0.331 0.384 0.496 0.813 0.8 0.253 0.307 0.348 0.403 0.523 0.858 Pr 1 0.272 0.332 0.378 0.440 0.570 0.938 5 0.45 0.585 0.669 0.792 1.043 1.736 10 0.570 0.730 0.851 1.013 1.344 2.236 Chapter 4 External Convective Heat and Mass Transfer 361 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press m = β = 0 corresponds to the case of flow parallel to a flat plate with constant free stream values (Blasius Solution). Eckert (1942) was the first to obtain the numerical results for eqs. (4.105) and (4.106) for constant free stream and wall temperature as presented in Table 4.2. In similarity solutions for both flow over a flat plate and over a wedge, it is assumed that there is no blowing or suction at the wall. It can be shown (see Problem 4.7) that similarity solutions with blowing or suction exist only if the vertical velocity at the wall changes along the flow according to the following equation: vw ∝ x ( m −1)/ 2 (4.109) where vw is the blowing or suction velocity. Clearly this requirement restricts the use of the similarity solution to special cases that satisfy eq. (4.109). 4.6.2 Coupled Mass, Momentum, and Heat Transfer Problems Coupled transport phenomena occur frequently in problems associated with evaporation, sublimation, absorption, or combustion. Sublimation over a flat plate can find its application in an analogy between heat and mass transfer and will be used here as an example to show the methodology used in solving these coupled problems. Figure 4.11 Phase diagram for solid-liquid and solid-vapor phase change. 362 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press When the pressure and temperature of ice are above the triple point and the ice is then heated, melting occurs. However, when the ice is exposed to moist air with a partial pressure of water below its triple point pressure, heating of the ice will result in a phase change from ice directly to vapor without going through the liquid phase. Spacecrafts and space suits can reject heat by sublimating ice into the vacuum of space. Another application for sublimation of ice is the preparation of specimens using freeze-drying for a scanning electronic microscope (SEM) or a transmission electronic microscope (TEM). This type of phase change is referred to as sublimation. The opposite process is deposition, which describes the process of vapor changing directly to solid without going through condensation. The phase-change processes related to solids can be illustrated by a phase diagram in Fig. 4.11. Sublimation and deposition will be the subjects of this sub-section. When a subcooled solid is exposed to its superheated vapor, as shown in Fig. 4.12(a), the vapor phase temperature is above the temperature of the solid-vapor interface and the temperature of the solid is below the interfacial temperature. The boundary condition at the solid-vapor interface is ∂T dδ ks s − hδ (T∞ − Tδ ) = ρ s hsv (4.110) ∂x dt where hδ is the convective heat transfer coefficient at the solid-vapor interface, hsv is the latent heat of sublimation, and δ is the thickness of the sublimated or deposited material. The interfacial velocity dδ/dt in eq. (4.110) can be either positive or negative, depending on the direction of the overall heat flux at the interface. While a negative interfacial velocity signifies sublimation, a positive interfacial velocity signifies deposition. When the vapor phase is superheated, as shown in Fig. 4.12(a), the solid-vapor interface is usually smooth and stable. (a) Subcooled solid exposed to superheated vapor (b) Superheated solid exposed to supercooled vapor Figure 4.12 Temperature distrbution in sublimation and deposition. Chapter 4 External Convective Heat and Mass Transfer 363 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press In another possible scenario, as shown in Fig. 4.12 (b), the solid temperature is above the interfacial temperature and the vapor phase is supercooled. The interfacial energy balance for this case can still be described by eq. (4.110). Depending on the degrees of superheating in the solid phase and supercooling in the vapor phase (the relative magnitude of the first and second terms in eq. (4.110)) both sublimation and deposition are possible. During sublimation, a smooth and stable interface can be obtained. During deposition, on the other hand, the interface is dendritic and unstable, because supercooled vapor is not stable. The solid formed by deposition of supercooled vapor has a porous structure. During sublimation or deposition, the latent heat of sublimation can be supplied from or absorbed by either the solid phase or the vapor phase, depending on the temperature distributions in both phases. Naphthalene sublimation is also a case whereby a heat transfer coefficient can be obtained through the measurement of a mass transfer coefficient and the analogy between heat and mass transfer (Eckert and Goldstein, 1976). The significant advantages of this method include its high accuracy and the simplicity of the experimental apparatus. In addition, the local heat transfer coefficient can be obtained by measuring the local sublimed depth of the specimen. Figure 4.13 shows the physical model of a sublimation problem, where a flat plate is coated with a layer of sublimable material and is subject to constant heat flux underneath (Zhang et al., 1996). A gas with the ambient temperature and mass fraction of sublimable material flows over the flat plate at a velocity of U∞. The heat flux applied from the bottom of the flat plate will be divided into two parts: one part is used to supply the latent heat of sublimation, and the other is transferred to the gas through convection. The sublimated vapor is injected into the boundary layer and is removed by the gas flow. The following assumptions are made in order to solve the problem: 1. The flat plate is very thin, so the thermal resistance of the flat plate can be neglected. 2. The gas is incompressible, with no internal heat source in the gas. 3. The sublimation problem is two-dimensional steady state. v U ∞ , T∞ , ω∞ y x u Boundary layer ′′ qw = constant Figure 4.13 Sublimation on a flat plate with constant heat flux. 364 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The governing equations for mass, momentum, energy, and species of the problem are ∂u ∂v + =0 ∂x ∂y (4.111) (4.112) (4.113) (4.114) u ∂u ∂u ∂2u +v =ν 2 ∂x ∂y ∂y u u ∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y ∂ω ∂ω ∂ 2ω +v =D 2 ∂x ∂y ∂y A no slip condition at the surface of the flat plat requires that u = 0 at y = 0 . For a binary mixture that contains the vapor, sublimable substance, and gas, the molar flux of the sublimable substance at the surface of the flat plate is [see eq. (1.109)] ρ D ∂ω  , y=0 (4.115) m ′′ = − 1 − ω ∂y Since the mass fraction of the sublimable substance in the mixture is very low, i.e., ω  1, the mass flux at the wall can be simplified to ∂ω  m ′′ = − ρ D , y=0 (4.116) ∂y  Sublimation at the surface causes a normal blowing velocity, vw = m ′′ / ρ . The normal velocity at the surface of the flat plate is therefore ∂ω v = vw = − D , y=0 (4.117) ∂y The energy balance at the surface of the flat plate is ∂T ∂ω ′′ −k − ρ hsv D = qw , y = 0 ∂y ∂y (4.118) Another reasonable, practical, representable boundary condition at the surface of the flat plate emerges by setting the mass fraction at the wall as the saturation mass fraction at the wall temperature. The mass fraction and the temperature at the surface of the flat plate have the following relationship (Kurosaki, 1974): ω = aT + b , y = 0 (4.119) where a and b are constants that depend on the sublimable material and its temperature. As y → ∞, the boundary conditions are u → U ∞ , T → T∞ , ω → ω∞ (4.120) Introducing the stream function ψ , ∂ψ ∂ψ u= v=− (4.121) ∂y ∂x Chapter 4 External Convective Heat and Mass Transfer 365 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the continuity equation (4.111) is automatically satisfied, and the momentum equation in terms of the stream function becomes ∂ψ ∂ 2ψ ∂ψ ∂ 2ψ ∂ 3ψ (4.122) − =ν 3 2 ∂y ∂x ∂y ∂x ∂y ∂y Similarity solutions for eq. (4.122) do not exist unless the injection velocity vw is proportional to x-1/2, and the incoming mass fraction of the sublimable substance, ω∞, is equal to the saturation mass fraction corresponding to the incoming temperature, T∞ (Kurosaki, 1974; Zhang et al., 1996). The governing equations cannot be reduced to ordinary differential equations. The local nonsimilarity solution proposed by Zhang et al. (1996) will be presented here. Defining the following similarity variables: U∞ x ψ ξ = , η=y , f= L 2ν Lξ 2ν U ∞ Lξ (4.123) ρ hsv D(ω − ω∞ ) k (T − T∞ ) , ϕ= θ= ′′ ′′ qw 2ν Lξ / U ∞ qw 2ν Lξ / U ∞ eqs. (4.122) and (4.113) – (4.114) become f ′′′ + ff ′′ = 2ξ ( f ′ ′ − f ′′ ) F F (4.124) θ ′′ + Pr( f θ ′ − f ′θ ) = 2 Pr ξ ( f ′Θ − θ ′F ) (4.125) (4.126) where prime ' represents partial derivative with respect to η, and all upper case variables represent partial derivatives of primary similarity variables with respect to ξ. ∂f ∂θ ∂ϕ F= , Θ= , Φ= (4.127) ∂ξ ∂ξ ∂ξ It can be seen from eqs. (4.124) – (4.126) that the similarity solution exists only if F = Θ = Φ = 0. In order to use eqs. (4.124) – (4.126) to obtain a solution for the sublimation problem, supplemental equations for F, Θ, and Φ must be obtained. Taking partial derivatives of eqs. (4.124) – (4.126) with respect to ξ and neglecting the higher order term, one obtains F ′′′ + Ff ′′ + F ′′f = 2 ( f ′ ′ − f ′′ ) F F (4.128) Θ′′ + Pr( F θ ′ + f Θ′ − F ′θ − f ′Θ) = 2 Pr ( f ′Θ − θ ′F ) ϕ ′′ + Sc( f ϕ ′ − f ′ϕ ) = 2Sc ξ ( f ′Φ − ϕ ′F ) (4.129) (4.130) The boundary conditions of eqs. (4.124) – (4.126) and eqs. (4.128) – (4.130) are f ′(ξ , 0) = 0, η = 0 2 f (ξ , 0) = − B ξ 1/ 2ϕ ′(ξ , 0) − ξ 3/ 2 Φ ′(ξ , 0)  , η = 0  3 f ′(ξ , ∞) = 1, η = ∞ F ′(ξ , 0) = 0, η = 0 Φ ′′ + Sc( F ϕ ′ + f Φ ′ − F ′ϕ − f ′Φ ) = 2Sc ( f ′Φ − ϕ ′F ) (4.131) (4.132) (4.133) (4.134) 366 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1 1  F (ξ , 0) = − B  ξ −1/ 2ϕ ′(ξ , 0) − ξ 1/ 2 Φ ′(ξ , 0)  , η = 0 3 2  F ′(ξ , ∞) = 0, η = ∞ θ ′(ξ , 0) + ϕ ′(ξ , 0) = −1, η = 0 (4.135) (4.136) (4.137) (4.138) (4.139) (4.140) (4.141) (4.142) (4.143) (4.144) (4.145) θ (ξ , ∞) = 0, η = ∞ Θ′(ξ , 0) + Φ ′(ξ , 0) = 0, η = 0 Θ(ξ , ∞) = 0, η = ∞ ϕ (ξ , 0) = ahsv 1 θ (ξ , 0) + ϕ sξ −1/ 2 , η = 0 c p Le ϕ (ξ , ∞) = 0, η = ∞ ah 1 ϕs Φ (ξ , 0) = sv Θ(ξ , 0) − 3 / 2 , η = 0 c p Le 2ξ Φ (ξ , ∞) = 0, η = ∞ where B= ′′ qw ρ hsvν 2ν L U∞ reflects the effect of injection velocity at the surface due to sublimation, and ρ h D (ωsat ,∞ − ω∞ ) ϕs = sv (4.146) ′′ qw 2ν L / U ∞ represents the effect of the mass fraction of the sublimable substance in the incoming flow. ωsat ,∞ is saturation mass fraction corresponding to the incoming temperature: ωsat ,∞ = aT∞ + b (4.147) The set of ordinary differential equations (4.124) – (4.126) and (4.128) – (4.130) with boundary conditions specified by eqs. (4.131) – (4.144) are boundary value problems that can be solved using a shooting method (Zhang et al., 1996). Figure 4.14 shows typical dimensionless temperature and mass fraction profiles obtained by numerical solution. It can be seen that the dimensionless temperature and mass fraction at different ξ are also different, which is further evidence that a similarity solution does not exist. Once the converged solution is obtained, the local Nusselt number based on the total heat flux at the bottom of the flat plate is Nux = ′′ hw x [qw / (Tw − T∞ )]x = = k k Re1/ 2 x 2θ (ξ , 0) (4.148) and the Nusselt number based on convective heat transfer is hx x  ∂T  θ ′(ξ , 0) Nu* = x = − Re1/ 2   =− x x k Tw − T∞  ∂y  y = 0 2θ (ξ , 0) (4.149) Chapter 4 External Convective Heat and Mass Transfer 367 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press θ ahsv = 0.1 cp ahsv = 1.0 cp 0.1 0.5 ahsv = 1.0 cp ahsv = 0.1 cp Figure 4.14 Temperature and mass fraction distributions (Zhang et al. 1996). ahsv = 1.0 cp ahsv = 1.0 cp Figure 4.15 Nusselt number based on convection and Sherwood number (Zhang et al. 1996). The Sherwood number is Shx = hm x ∂ω x =− D ωw − ω∞ ∂y =− y =0 ϕ ′(ξ , 0) Re1/ 2 x 2ϕθ (ξ , 0) (4.150) Figure 4.15 shows the effect of blowing velocity on the Nusselt number based on convective heat transfer and the Sherwood number for ϕsat ,∞ = 0 , i.e., the mass fraction of sublimable substance is equal to the saturation mass fraction corresponding to the incoming temperature. It can be seen that the effect of blowing velocity on mass transfer is stronger than that on heat transfer. 368 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 4.7 Integral Methods As noted before, the similarity solution provides an exact analytical solution for laminar boundary layer conservation equations; however, there are limitations in terms of geometry and boundary conditions, as well as laminar flow restrictions. The integral method gives approximately closed form solutions, which have much less limitation in terms of their geometry and boundary conditions. It can also be applied to both laminar and turbulent flow situations. The integral method easily provides accurate answers (not exact) for complex problems. Using the integral method, one usually integrates the conservative differential boundary layer equation over the boundary layer thickness by assuming a profile for velocity, temperature, and concentration, as needed. The better the approximate shape for the profile, such as velocity and temperature, the better the prediction of drag force and heat transfer (friction coefficient or heat transfer coefficient). The integral methodology has been applied to a variety of configurations to solve transport phenomena problems (Schlichting and Gersteu, 2000). To illustrate the integral methodology, it will be applied to flow and heat transfer over a wedge with non-uniform temperature and blowing at the wall. Consider two dimensional laminar steady flow with constant properties over a wedge, as shown in Figure 4.16, with an unheated starting length, x0. δT δ U∞ T∞ y x x0 β Tw Figure 4.16 Momentum and heat transfer over a wedge with an unheated starting length. Chapter 4 External Convective Heat and Mass Transfer 369 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The governing boundary layer equations for mass, momentum and energy for constant property, steady state, and laminar flow as well as boundary conditions for convective heat transfer over a wedge are presented below: Continuity equation ∂u ∂v + =0 ∂x ∂y (4.151) Momentum equation u ∂u ∂u ∂2u dU +v =ν 2 +U ∂x ∂y ∂y dx u ∂T ∂T ∂ 2T +v =α 2 ∂x ∂y ∂y (4.152) Energy equation (4.153) Boundary conditions u ( 0, x ) = 0 T ( ∞, x ) = T∞ T ( 0, x ) = T∞ T ( 0, x ) = Tw 1 ∂p u (δ , x ) = U ( x ) for x < x0 for x > x0 dU dx v ( 0, x ) = vw (4.154) It should be noted that U is known from potential flow theory: − ρ ∂x =U If U is constant, then dU =0 dx (4.155) As for the case of flow over a flat plate, If x0 = 0, U = constant and vw = 0, the problem will be similar to the case presented using the similarity solution before. In integral methods, it is customary to assume a profile for u and obtain v from the continuity eq. (4.151). Let us integrate eq. (4.151) with respect to y from y = 0 to y = δ. The velocity and temperature field outside δ is uniform.   δ 0 δ 0 δ ∂v ∂u dy +  dy = 0 0 ∂y ∂x (4.156) The second term can be easily integrated. ∂v dy = v y =δ − v y = 0 = vδ − vw ∂y (4.157) Combining (4.156) and (4.157) we get:  δ 0 ∂u dy = vw − vδ ∂x (4.158) Applying Leibnitz’s formula to the left hand side of (4.158) yields: 370 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂δ dδ udy − u ( x , δ ) = vw − vδ ∂x 0 dx (4.159) which can be rearranged to obtain: dδ ∂ vδ = vw + U dx − ∂x (  udy ) δ 0 (4.160) The momentum equation (4.152) can be rearranged to the following form:  ∂u ∂v  ∂u 2 ∂vu dU ∂ 2u + − u +  =U +ν 2 dx ∂x ∂y ∂y  ∂x ∂y  (4.161) where term in parenthesis on the left-hand side is zero because of the continuity equation. The integration of eq. (4.161) from y = 0 to y = δ gives:   δ 0 δ 0 2 δ ∂vu δ δ∂ u ∂u 2 dU dy +  dy =  U dy + ν  dy 0 ∂y 0 0 ∂y 2 ∂x dx (4.162) Upon further integration and simplification, the above equation reduces to δ ∂u 2 dU ∂u ∂u dy + vu δ − vu 0 =  U dy + ν −ν 0 ∂x ∂y δ ∂y 0 dx (4.163) Using eq. (4.160) for vδ, the no slip boundary condition at the wall u(0,x) = 0, and assuming no velocity gradient at the outer edge of the boundary layer at y = δ we get: 2 δ ∂u δ τw ∂δ dU 2 dδ (4.164) 0 ∂x dy + Uvw + U dx − U ∂x 0 udy = − ρ + 0 U dx dy where ( ) τw = μ ∂u ∂y 0 (4.165) is the shear stress at the wall. Applying Leibnitz’s rule and rearranging will provide the final form. δ δ τ w + ρUvw ∂δ dU dU (4.166) 0 u ( u − U ) dy  + 0 udy dx − 0 U dx dy = − ρ    ∂x  The only dependent unknown variable in the above equation is u since v is eliminated using continuity. U, τw, and vw should be known quantities. Equation (4.166) can be further rearranged. ∂  2 δ u u    δ u   dU τ w + ρUvw (4.167) U 0 U 1 − U  dy  +  0 1 − U  dy  U dx = ∂x  ρ     It is customary to assume a third order polynomial equation for the velocity profile in order to obtain a reasonable result, u = c1 + c2 y + c3 y 2 + c4 y3 (4.168) where c1, c2, c3, and c4 are constants and can be obtained from boundary conditions for velocity and shear stress at the wall and outer edge. Once the constants are obtained, they are substituted into the momentum integral equation (4.167) and are used to solve for the momentum boundary layer thickness, δ. Similarly, the energy equation (4.153) can be rearranged into the following form to make the integration process easier: ( ) Chapter 4 External Convective Heat and Mass Transfer 371 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  ∂u ∂u  ∂uT ∂vT ∂ 2T + −T  +  =α 2 ∂x ∂y ∂y  ∂x ∂y  (4.169) Let’s integrate the above equation from y = 0 to y = δT, knowing that the term in parenthesis on the left hand side is zero because of the continuity equation.  δT 0 2 δ T ∂vT δT ∂ T ∂uT dy +  dy = α  dy 0 0 ∂y 2 ∂x ∂y (4.170) Using integration and the continuity equation to obtain vδ, eq. (4.160), and assuming there is no temperature gradient at the outer edge of the thermal boundary layer, we get the following equation: δ∂ dδ T ∂δ −k ∂T (4.171) 0 ∂x ( uT ) dy − T∞ ∂x 0 udy + T∞ vw − vwTw + T∞ uδ dx = ρ c p ∂y y =0 T T T Using Leibnitz’s rule and rearranging gives: ∂  δT u ( T − T∞ ) dy  +    ∂x  0 ( δT 0 udy ) dT = ρqc′′ dx ∞ w p + vw ( Tw − T∞ ) (4.172) Once again, the integral form of the energy equation is in terms of the unknown temperature, assuming the velocity profile is known. Similar to the momentum integral equation, a temperature profile should be assumed and substituted into the integral energy equation (4.172) to obtain δT. To illustrate the procedure, we use the above approximation to solve the classic problem of flow and heat transfer over a flat plate when U = U∞ = constant, with no blowing or suction at the wall, and with constant wall and flow stream temperature. The momentum and energy integral equations (4.167) and (4.172) reduce to the following forms using the above assumptions: ∂  2 δ u u   τw (4.173) U 0 U 1 − U  dy  = ρ ∂x    ′′ qw ∂  δT 0 u (T − T∞ ) dy  = ρ c p    ∂x  (4.174) Let’s assume a polynomial velocity profile is a third degree polynomial function with the following boundary conditions. u ( 0) = 0 (4.175) u (δ ) = U y =δ (4.176) (4.177) (4.178) ∂u ∂y =0 ∂ 2u ∂y 2 =0 y =δ It is assumed that shear stress at the boundary layer edge is zero, which is a good approximation for this configuration. Equation (4.178) is obtained by using eq. (4.177) and applying the x-direction momentum equation at the boundary layer edge. Upon applying eqs. (4.175) – (4.178) in eq. (4.168), one obtains four 372 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press equations containing four unknowns (c1, c2, c3, and c4). After solving these four eauations, the final velocity profile can be obtained as: u 3 y 1 y =  −   , U 2δ  2δ  3 y≤δ (4.179) Shear stress at the wall, τw, is calculated using eq. (4.179) ∂u 3μU = τw = μ ∂y y = 0 2δ (4.180) Substituting eqs. (4.179) and (4.180) into eq. (4.173), and performing the integration we get d  39U 2δ  3ν U (4.181)  = dx  280  2δ Integrating the above equation and assuming δ = 0 at x = 0, we get 1/ 2  280ν x  δ = (4.182)   13U  = 4.64 Re1/ 2 x or δ x (4.183) ∂u ∂y The friction coefficient is found as before cf = τω U2 ρ∞ 2 cf 2 = μ = y =0 U2 ρ∞ 2 (4.184) or using eqs. (4.180) and (4.182) 0.323 Re1/ 2 x (4.185) The predictions of the momentum boundary layer thickness and friction coefficient, cf, by the integral method are 7% and 3% lower than the exact solution obtained using the similarity method, respectively. Using the same general third order polynomial equation for a temperature profile with the following boundary conditions: T ( 0 ) = Tw = constant (4.186) T (δ T ) = T∞ = constant (4.187) (4.188) (4.189) ∂T ∂y =0 y =δ T ∂ 2T ∂y 2 =0 y =0 the temperature profile can be obtained as: T − T∞ 3 y 1 y  = 1− +   , y ≤ δT 2 δT 2  δT  Tw − T∞ 3 (4.190) Chapter 4 External Convective Heat and Mass Transfer 373 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press and the heat flux at the wall is ′′ qw = −k ∂T ∂y = y =0 3k (Tw − T∞ ) 2 δT (4.191) Substitution of eqs. (4.189), (4.190) and (4.191) into eq. (4.174), and approximate integration for δ T / δ < 1 yields: 1 u  T −T ∂ δT ∂ ∞ u ( T − T∞ ) dy = U ∞ ( Tw − T∞ ) δ T 0  U ∞  Tw − T∞ ∂x 0 ∂x    δ2 ∂  3 δT 2  1 − T 2  U ( Tw − T∞ )  =   ∂x  20 δ  14δ     ′′ qw 3α = = (Tw − T∞ ) ρ c p 2 δT   y   d     δT   (4.192) Upon further simplification we get 2 2 δ T   10ν 1 d δT   1 −  = dx  δ  14δ 2   Pr U δ T   The only unknown in the above equation is δT since δ is known. Assuming ζ = δ T / δ , eq. (4.193) becomes d  2  ζ 2   10 ν 1 ζ δ 1 −  = dx   14   Pr U ∞ ζδ (4.193) (4.194) where term ζ 2 /14 can be neglected since it is much less than 1. The solution of eq. (4.194) for ζ = 0 at x = x0 yields ζ= Pr −1/3   x0  1 − 1.026   x   −k ∂T ∂y 3/ 4 1/3     (4.195) The local heat transfer coefficient can now be calculated since δT is known. hx = ′′ qw 3k y =0 = = 2 δT Tw − T∞ Tw − T∞ hx = 3k 2 ζδ T (4.196) or (4.197) Using ζ from eq. (4.195) and δ from eq. (4.69) to calculate the local Nusselt number, Nu x = hx x 0.332 Pr1/3 Re1/ 2 x = 1/3 k   x0 3/ 4  1 −     x    (4.198) The above equation without unheated starting length (x0 = 0) reduces to the exact solution obtained by the similarity solution. Nu x = 0.332 Pr1/3 Re1/ 2 (4.199) x 374 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 4.8 Computational Methodologies for Forced Convection The convection problems that can be solved analytically – some of which were discussed in the preceding sections – are limited to the cases for which the boundary layer theory is valid. Even for cases where the boundary layer theory can be used to simplify the governing equations, analytical solutions can be obtained only if the similarity solution exists or if integral approximate solutions can be obtained. The analytical solution can be very complicated for cases with variable wall temperature or heat flux. When the boundary layer theory cannot be applied, or a simple analytical solution based on boundary layer theory cannot be obtained, a numerical solution becomes the desirable approach. Convection problems are complicated by the presence of advection terms in the governing equations. For incompressible flow, the problem is further complicated by the fact that there is not a governing equation for pressure. This section will cover (1) the numerical solution of the convection-diffusion equation with a known flow field, and (2) the algorithm to determine the flow field. The governing equation for a convection problem in the Cartesian coordinate system can be expressed into the following generalized form: ∂ ( ρφ ) ∂ ( ρ uφ ) ∂ ( ρ vφ ) ∂ ( ρ wφ ) ∂t + ∂x + ∂y + ∂z ∂  ∂φ  ∂  ∂φ  ∂  ∂φ  = Γ  + Γ  + Γ  + S ∂x  ∂x  ∂y  ∂y  ∂z  ∂z  (4.200) where φ is a general variable that can represent the directional components of velocity (u, v, or w), temperature, or mass concentration. The general diffusivity Γ can be viscosity, thermal conductivity, or mass diffusivity. S is the volumetric source term. Compared with the governing equation for heat conduction, the convection terms seem to be the only new terms introduced into eq. (4.200). As will become evident later, the contribution of the convection term on the overall heat transfer depends on the relative scale of convection over diffusion. Therefore, these two effects are always handled as one unit in the derivation of the discretization scheme (Patankar, 1980; Patankar, 1991). While analytical solutions of convection problems presented in the preceding sections based on boundary layer theory are few, the numerical solutions of convection based on boundary layer theory have been abundant (Tao, 2001; Minkowycz et al., 2006). With significant advancement of computational capability, the numerical solution of convection problems can now be performed by solving the full Navier-Stokes equation. The algorithms that will be discussed in this section are therefore based on the solution of the full Navier-Stokes equation and the energy equation, not on the boundary layer equations. Chapter 4 External Convective Heat and Mass Transfer 375 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 4.8.1 One-Dimensional Steady-State Convection and Diffusion Exact Solution The objective of this subsection is to introduce various discretization schemes of the convection-diffusion terms through discussion of the one-dimensional steady state convection and diffusion problem. For a one-dimensional steady-state convection and diffusion problem, the governing equation is d ( ρ uφ ) d  dφ  (4.201) = Γ  dx dx  dx  where the velocity, u, is assumed to be known. Both density, ρ, and diffusivity, Г, are assumed to be constants. The continuity equation for this one-dimensional problem is d ( ρ u) =0 (4.202) dx Equation (4.201) is subject to the following boundary conditions: φ = φ0 , x = 0 (4.203) φ = φL , x = L (4.204) By introducing the following dimensionless variables φ − φ0 x Φ= , X= (4.205) − φ0 L φL the one-dimensional steady-state convection and diffusion problem can be nondimensionalized as Pe dΦ d 2 Φ = dX dX 2 Φ = 0, X = 0 Φ = 1, X = 1 (4.206) (4.207) (4.208) (4.209) where Pe = ρ uL Γ is the Peclet number that reflects the relative level of convection and diffusion. Pe becomes zero for the case of pure diffusion and becomes infinite for the case of pure advection. The exact solution of eqs. (4.206) - (4.208) can be obtained as φ − φ0 exp(PeX ) − 1 (4.210) Φ= = exp(Pe) − 1 φL − φ0 which will be used as a criterion to check the accuracy of various discretization schemes. Central Difference Scheme Integrating eq. (4.201) over the control volume P (shaded area in Fig. 4.17), one obtains 376 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (δx)w (δx)w- (δx)w+ w (δx)e (δx)e- (δx)e+ e Control Volume W P (Δx)P Figure 4.17 Control volume for one-dimensional problem. E  dφ  ( ρ uφ )e − ( ρ uφ ) w =  Γ   dx e  dφ  −Γ   dx  w (4.211) The right-hand side of eq. (4.211) can be obtained by assuming the distribution of φ between any two neighboring grid points is piecewise linear (see Chapter 3), i.e., φ E − φP  dφ  Γ  = Γe (δ x )e  dx e φP − φW  dφ  Γ  = Γw (δ x ) w  dx w where Γe and Γw are the diffusivities at the faces of the control volume. To ensure that the flux of φ across the faces of the control volume is continuous, the harmonic mean diffusivity at the faces should be used. To evaluate the left hand side of eq. (4.211), it is necessary to know the values of φ at the faces of the control volume. If the piecewise linear profile of φ is chosen, it follows that φ +φ φ +φ φe = E P , φw = P W 2 2 Therefore, eq. (4.211) becomes φ +φ φ −φ φ + φP φ − φP ( ρ u) e E − ( ρ u) w P W = Γ e E − Γw P W 2 2 (δ x )e (δ x ) w Defining the mass flux and diffusive conductance F = ρ u, D = Γ δx (4.212) (4.213) (4.214) (4.215) (4.216) (4.217) eq. (4.212) can be rearranged as aP φP = aE φE + aW φW where 1 aE = De − Fe 2 1 aW = Dw + Fw 2 aP = aW + aE + ( Fe − Fw ) Chapter 4 External Convective Heat and Mass Transfer 377 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press This scheme is termed the central difference scheme because the values of φ at the faces of the control volume are taken as the averaged value between two grid points. The continuity equation (4.202) requires that Fe = Fw and therefore, eq. (4.217) reduces to aP = aW + aE To evaluate the performance of the central difference scheme, let us consider the case of a uniform grid, i.e., (δ x )e = (δ x ) w = δ x , for which case eq. (4.212) can be rearranged as φP = 1 − Δ  φE + 1 + Δ  φW  2  2 2   1  Pe   Pe   (4.218) where Pe Δ = ρ uδ x Γ = F D (4.219) is the Peclet number using grid size as the characteristic length, which is referred to as the grid Peclet number. The grid Pe is a ratio of the strength of convection over diffusion. To ensure stability of the discretization scheme, the value of φP should always fall between φE and φW , which requires that the coefficients, φE and φW , are positive, i.e., Pe Δ ≤ 2 (4.220) This is the criterion for stability of the central difference scheme. It can be demonstrated that the central difference becomes unstable if eq. (4.220) is violated (see Problem 4.22). The fact that the central difference scheme is stable under small grid Peclet number indicates that the central difference scheme is accurate only if the convection is not very significant. Upwind Scheme The central difference scheme assumes that the effects of the values of φ at two neighboring grid points on the value of φ at the face of the control volume are equal. This assumption is valid only if the effect of diffusion is dominant. If, on the other hand, the convection is dominant, one can expect that the effect of the grid point upwind is more significant than that of the point downwind. If we can assume that the value of φ at the face of the control volume is dominated by the value of φ at the grid point at the upwind side and that the effect of the value of φ at the downwind side can be neglected, the two terms on the left hand side of eq. (4.211) can be expressed as  F φ , Fe > 0 ( ρ uφ )e =  e P  FeφE , Fe < 0  F φ , Fw > 0 ( ρ uφ ) w =  w W  FwφP , Fw < 0 The above two equations can be expressed in the following compact form: 378 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ( ρ uφ )e = φP  Fe , 0 − φE  − Fe , 0 ( ρ uφ ) w = φW  Fw , 0 − φP  − Fw , 0 where the operator  A, B  denotes the greater of A and B (Patankar, 1980). Substituting the above expression into the left hand side of eq. (4.211) and using central difference for the right hand side of eq. (4.211), the discretized equation becomes aP φP = aE φE + aW φW (4.221) where aE = De +  − Fe , 0 (4.222) (4.223) aP = aW + aE + ( Fe − Fw ) (4.224) The above scheme is referred to as the upwind scheme because the value of φ at the grid point on the upwind side was used as the value of φ at the face of the control volume to discretize the convection term. The upwind scheme ensures that the coefficients in eq. (4.221) are always positive so that a physically unrealistic solution can be avoided. Hybrid Scheme aW = Dw +  Fw , 0 The upwind scheme uses the value of φ from the grid point at the upwind side as the value of φ at the face of the control volume regardless of the grid Peclet number. While this treatment can yield accurate results for cases with high Peclet number, the result will not be accurate for cases where the grid Peclet number is near zero; for which cases the central difference scheme can produce better results. Spalding (1972) proposed a hybrid scheme that uses the central difference scheme when Pe Δ ≤ 2 and the upwind scheme when PeΔ > 2 . To observe the difference between the central difference and upwind schemes, the coefficient for the east neighboring grid point, eqs. (4.215) and (4.222), can be rewritten as 1 aE / De = 1 − Pe Δe , Central difference scheme 2 aE / De = 1 +  − Pe Δe , 0, Upwind scheme (4.225) (4.226) The hybrid scheme can then be expressed as  − Pe Δe 1  aE / De = 1 − Pe Δe 2 0  Pe Δe < −2 − 2 ≤ Pe Δe ≤ 2 Pe Δe > 2 which can be rewritten in the following compact form   1 aE / De =  − Pe Δe ,1 − Pe Δe , 0   2   (4.227) Chapter 4 External Convective Heat and Mass Transfer 379 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The coefficient for the west neighbor grid point can be obtained using a similar approach.   1 aW / Dw =  Pe Δw ,1 + Pe Δw , 0   2   (4.228) The above hybrid scheme combines the advantages of the central difference and upwind schemes to yield better results for cases where Pe Δ → ∞ or Pe Δ  0 . However, there is still room for improvement of the solution when Pe Δ is near 2 (see Problem 4.23). Exponential and Power Law Schemes Since the exact solution of eq. (4.201) exists, one can reasonably expect that an accurate scheme can be derived if the result of the exact solution, eq. (4.210), is utilized. Equation (4.201) can be rewritten as d dφ  (4.229)  ρ uφ − Γ =0 dx  dx  Defining the total flux of φ due to convection and diffusion dφ J = ρ uφ − Γ dx (4.230) eq. (4.229) becomes dJ =0 dx (4.231) Integrating eq. (4.231) over the control volume P (shaded area in Fig. 4.17), yields Je = Jw (4.232) Instead of assuming piecewise linear distribution of φ as with central difference scheme or assuming φ at the face of the control volume is equal to the value of φ at the grid point on the upwind side in the upwind scheme, the distribution of φ between grid points can be taken as that obtained from the exact solution, eq. (4.210). Applying eq. (4.210) between grid points E and P, we have φ ( x ) − φP exp[Pe Δe ( x − x P ) / (δ x )e ] − 1 (4.233) = φ E − φP exp(Pe Δe ) − 1 Substituting eq. (4.233) into eq. (4.230) and evaluating the result at x = xe , the total flux of φ at the face of control volume becomes  φP − φ E  J e = Fe φP +  exp(Pe Δe ) − 1   (4.234) Similarly, the total flux at the west face of the control volume is  φW − φP  J w = Fw φW +   exp(Pe Δw ) − 1  (4.235) 380 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 5 Hybrid scheme, eq. (4.227) 4 aE / De 3 2 Upwind scheme, eq. (4.222) Central difference scheme, eq. (4.215), -2 < Pe < 2 1 Exact solution, eq. (4.240) 0 -5 -4 -3 -2 -1 0 Pe 1 2 3 4 5 Figure 4.18 Comparison among different schemes. Substituting eqs. (4.234) and (4.235) into eq. (4.232) and rearranging the resulting equation yields aP φP = aE φE + aW φW (4.236) where aE = aW = Fe exp(Pe Δe ) − 1 Fw exp(Pe Δw ) exp(Pe Δw ) − 1 (4.237) (4.238) (4.239) Equations (4.237) and (4.238) can be rewritten in a format similar to that of eqs. (4.225) – (4.228), i.e., aP = aW + aE + ( Fe − Fw ) aE / De = aW / Dw = Pe Δe exp(Pe Δe ) − 1 (4.240) (4.241) Pe Δw exp(Pe Δw ) exp(PeΔw ) − 1 The comparison of aE / De for different schemes is shown in Fig. 4.18. It can be seen that the hybrid scheme can be viewed as an envelope of the exponential scheme. The hybrid scheme is a good approximation if the absolute value of the grid Peclet number is either very large or near zero. While the exponential scheme is accurate, the computational time is much longer than for the central difference, upwind or hybrid schemes. Patankar (1981) proposed a power law scheme that has almost the same accuracy as the exponential scheme but a substantially shorter computational time. The coefficient of the neighbor grid point on the east side can be obtained by Chapter 4 External Convective Heat and Mass Transfer 381 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press PeΔe < −10  − Pe Δe (1 + 0.1Pe )5 − Pe − 10 ≤ Pe Δe < 0  Δe Δe aE / De =  (1 − 0.1Pe Δe )5 0 ≤ PeΔe ≤ 10   0 Pe Δe > 10  which can be rewritten in the following compact form  aE / De =  0, (1 − 0.1 Pe Δe   )  + 0, −Pe   5 Δe (4.242) A Generalized Expression The above discretization schemes can be expressed in a single generalized form. The total flux J at the interface between two grid points that were defined in eq. (4.230) can be used to define: J dφ J* = = Pe Δφ − (4.243) d(x / δ x) Γ/δx which relates to the values of φ at grid points i and i+1 (see Fig. 4.19). The first term on the right side of eq. (4.243) will be related to some weighted average of φi and φi +1 , and the second term will be related to the difference between φi and φi +1 . Thus, one can express J * the total flux as (Patankar, 1980) J * = B (Pe Δ )φi − A(Pe Δ )φi +1 (4.244) where A and B are dimensionless coefficients that are functions of the grid Peclet number. If the field of φ is uniform, we will have dφ / dx = 0 and eq. (4.243) becomes J * = Pe Δφi = Pe Δφi +1 (4.245) Comparing eqs. (4.244) and (4.245) yields (4.246) For the grid system shown in Fig. 4.19, if we reconsider the problem in a reversed coordinate system x ' ( x ' = − x ), the grid Peclet number will become −Pe Δ and J * becomes J * = B (−Pe Δ )φi +1 − A(−Pe Δ )φi (4.247) B (Pe Δ ) − A(Pe Δ ) = Pe Δ δx i J i+1 x Figure 4.19 Total flux at the interface between grid points i and i+1. 382 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The symmetric properties of A and B can be obtained by comparing eqs. (4.244) and (4.247), i.e., A(−Pe Δ ) = B (Pe Δ ) (4.248) B (− Pe Δ ) = A(Pe Δ ) (4.249) * For the exponential schemes discussed above, one can obtain J from eq. (4.234) or (4.235), i.e.,  φi − φi +1  J * = Pe Δ φi +   exp(Pe Δ ) − 1  = exp(Pe Δ )Pe Δ Pe Δ φi − φi +1 exp(Pe Δ ) − 1 exp(Pe Δ ) − 1 exp(Pe Δ )Pe Δ Pe Δ , A= exp(Pe Δ ) − 1 exp(Pe Δ ) − 1 Comparing the above expression with eq. (4.244), one obtains B= It can be verified that the above A and B satisfy eqs. (4.246), and (4.248) – (4.249). The implication of the above properties of A and B is that if the function A(PeΔ) for the case that PeΔ > 0 is known, the expressions of A and B for all PeΔ can be obtained. For example, if PeΔ < 0 , eq. (4.246) can be used to obtain A(Pe Δ ) = B (Pe Δ ) − Pe Δ Substituting eq. (4.248) into the above equation yields A(Pe Δ ) = A(− Pe Δ ) − Pe Δ Considering −Pe Δ = PeΔ for the case that Pe Δ < 0 , the above expression can be rewritten as A(Pe Δ ) = A( Pe Δ ) + Pe Δ for Pe Δ < 0 Since A(PeΔ ) = A( Pe Δ ) for Pe Δ > 0 the following expression for A under any grid Peclet number can be expressed as (4.250) A(Pe Δ ) = A( Pe Δ ) +  −Pe Δ , 0 Similarly, the expression of B for any grid Peclet number can be expressed as (see Problem 4.24). B (Pe Δ ) = A( Pe Δ ) +  Pe Δ , 0 (4.251) Therefore, different discretization schemes for the convection-diffusion terms can be characterized by different A( Pe Δ ) . To derive the generalized formula for different discretization schemes, let us begin from eq. (4.232), i.e., * * J e De = J w D w (4.252) The total fluxes at the faces of the control volumes can be obtained from eq. (4.244), i.e., * J e = B (Pe Δe )φP − A(Pe Δe )φE (4.253) * J w = B (Pe Δw )φW − A(Pe Δw )φP (4.254) Chapter 4 External Convective Heat and Mass Transfer 383 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Substituting the above expressions into eq. (4.252) and rearranging the resulting equation yields [ De B(PeΔe ) + Dw A(PeΔw )]φP = De A(PeΔe )φE + Dw B(PeΔw )φW which can be rearranged as aP φP = aE φE + aW φW (4.255) where aE = De A(Pe Δe ) = De { A( Pe Δe ) +  − Pe Δe , 0} (4.256) aW = Dw B (PeΔw ) = Dw { A( Pe Δw ) +  PeΔw , 0} aP = aE + aW + ( Fe − Fw ) (4.257) (4.258) In arriving at eqs. (4.256) and (4.257), A and B were obtained from eqs. (4.250) and (4.251). At this point, it is apparent that different discretization schemes can be characterized by different expressions for A(|PeΔ|). By comparing eqs. (4.256) and (4.257) with different expressions of aE and aW for different schemes, the corresponding A(|PeΔ|) for different schemes can be summarized in Table 4.3 and plotted in Fig. 4.20. It should be noted that the difference between 1 0.8 Upwind 0.6 A(|PeΔ|) Power law 0.4 0.2 Exponential Hybrid 0 0 -0.2 1 2 3 4 |PeΔ| 5 6 7 Central difference Figure 4.20 Comparison of A(|PeΔ|) for different schemes. Table 4.3 Summary of A(|PeΔ|) for different schemes Scheme Central difference A(|PeΔ|) 1 − 0.5 Pe Δ 1  0,1 − 0.5 Pe  Δ  Upwind Hybrid Exponential Power law Pe Δ / [exp( Pe Δ ) − 1]  0,(1 − 0.1 Pe )5    Δ   384 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the power law and exponential scheme is exaggerated for clear presentation. The generalized formula represented by eqs. (4.255) – (4.258) will be very helpful to develop a generalized computer code for all schemes. A special module or subroutine can be written for different schemes. 4.8.2 Multidimensional Convection and Diffusion Problems Two-dimensional problem The heat transfer problems discussed in the preceding subsection are steady-state convection-diffusion problems with the general variable φ varying in one dimension only. We now turn our attention to the unsteady state two-dimensional convection-diffusion problem which includes a source term S. The problem is described by ∂ ( ρφ ) ∂ ( ρ uφ ) ∂ ( ρ vφ ) ∂  ∂φ  ∂  ∂φ  + + = (4.259) Γ  + Γ  + S ∂t ∂x ∂y ∂x  ∂x  ∂y  ∂y  + ∂J y ∂y =S which can be rewritten as ∂ ( ρφ ) ∂J x ∂t + ∂x (4.260) where ∂φ ∂x ∂φ J y = ρ vφ − Γ ∂y J x = ρ uφ − Γ N (δx)w Jn n Δy W w Jw P e Je E (δx)e (4.261) (4.262) (δy)n Control Volume Js Δx S s (δy)s Figure 4.21 Control volume for two-dimensional problem. Chapter 4 External Convective Heat and Mass Transfer 385 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Integrating eq. (4.260) with respect to t in the interval of (t, t+Δt) and over the control volume P (see Fig. 4.21), we have  s e n w t + Δt t t + Δt n e ∂J t + Δt w n ∂J ∂ ( ρφ ) dtdxdy + t s w x dxdydt + t e s y dydxdt ∂t ∂x ∂y e = t +Δt t  s n w ( SC + S P φ )dxdydt where the source term is treated as a linear function of φ . Assuming the total fluxes are uniform on all faces of the control volume and employing fully-implicit scheme, the above equation becomes 00 e w n s ( ρ P φP − ρ P φP )Δx Δy + ( J x − J x )ΔyΔt + ( J y − J y )Δx Δt = ( SC + S P φP )Δx ΔyΔt where the superscript 0 represents the values at the previous time step. e w n Introducing the integrated total fluxes J e = J x Δy , J w = J x Δy , J n = J y Δx and s J s = J y Δx , and dividing the above equation by Δt yields 00 ( ρ P φP − ρ P φP ) Δx Δy + ( J e − J w ) + ( J n − J s ) = ( SC + S P φP )Δx Δy Δt (4.263) Defining the following integrated total flux * Je = Je J J J * * * , Jw = w , Jn = n , Js = s De Dw Dn Ds where De = Γe Γw Γn Γs Δy, Dw = Δy, Dn = Δx , Ds = Δx (δ x )e (δ x ) w (δ y)n (δ y) s the integrated fluxes at the east and west faces of the control volume can be evaluated using eqs. (4.253) and (4.254), i.e., * J e = De J e = De [ B (Pe Δe )φP − A(Pe Δe )φE ] * J w = Dw J w = Dw [ B (Pe Δw )φW − A(Pe Δw )φP ] Substituting eq. (4.246) into the two equations above yields * J e = De J e = De [ A(Pe Δe )φP + Pe ΔeφP − A(Pe Δe )φE ] * J w = Dw J w = Dw [ B (Pe Δw )φW − B(Pe Δw )φP + PeΔwφP ] Similarly, the integrated total flux at the north and south faces of the control volume can be expressed as * J n = Dn J n = Dn [ A(Pe Δn )φP + Pe ΔnφP − A(Pe Δn )φN ] * J s = Ds J s = Ds [ B (Pe Δs )φS − B (Pe Δs )φP + Pe ΔsφP ] Substituting the above four integrated total fluxes into eq. (4.263), we have {ρ P Δx Δy / Δt − S P Δx Δy + De A(PeΔe ) + Dw B(PeΔw ) + Dn A(Pe Δn ) + Ds B(Pe Δs ) + ( Fe − Fw ) + ( Fn − Fs )} φP = De A(Pe Δe )φE + Dw B (Pe Δw )φW + Dn A(Pe Δn )φ N 0 0 + Ds B (Pe Δs )φS + SC Δx Δy + ( ρ P Δx Δy / Δt )φP (4.264) where Fe = De Pe Δe = ( ρ u)e Δy, Fw = Dw Pe Δw = ( ρ u) w Δy 386 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Fn = Dn Pe Δn = ( ρ v)n Δx , Fs = Ds Pe Δs = ( ρ v) s Δx are the mass flow rates at the four faces of the control volume. Equation (4.264) can be rearranged to obtain the final form of the following discretized equation: aP φP = aE φE + aW φW + aN φ N + aS φS + b (4.265) where aE = De A(Pe Δe ) = De { A( Pe Δe ) +  − Pe Δe , 0} (4.266) aW = Dw B (PeΔw ) = Dw { A( Pe Δw ) +  PeΔw , 0} aN = Dn A(Pe Δn ) = Dn { A( Pe Δn ) +  −Pe Δn , 0} (4.267) (4.268) (4.269) (4.270) aS = Ds B (Pe Δs ) = Ds { A( Pe Δs ) +  Pe Δs , 0} 00 b = a P φ P + S C Δ x Δy aP = aE + aW + aN + aS + ( ρ P / ρ )a 0 P 0 aP = 0 ρ P Δx Δy 0 P − S P Δx Δy + ( Fe − Fw ) + ( Fn − Fs ) Δt (4.271) (4.272) If the continuity equation is satisfied, eq. (4.271) can be simplified as (see Problem 4.4) 0 aP = aE + aW + aN + aS + aP − S P Δx Δy (4.273) Similar to the case of one-dimensional convection-diffusion, different discretization schemes for the discretized equations (4.265) – (4.272) can be obtained by using different expressions for A(|PeΔ|) from Table 4.3. Three-dimensional problem The discretized equation for a transient three-dimensional convection-diffusion problem can be obtained by integrating the conservation equation with respect to t in the interval of (t, t+Δt) and over the three-dimensional control volume P (formed by considering two additional neighbors at top, T, and bottom, B). The final form of the governing equation is (Patankar, 1980) aP φP = aE φE + aW φW + aN φ N + aS φS + aT φT + aBφB + b (4.274) where aE = De { A( Pe Δe ) +  − Pe Δe , 0} (4.275) aN = Dn { A( Pe Δn ) +  −Pe Δn , 0} aT = Dt { A( PeΔt ) +  −PeΔt , 0} aB = Db { A( Pe Δb ) +  Pe Δb , 0} b = a φ + SC Δx ΔyΔz 00 PP 0 P aW = Dw { A( Pe Δw ) +  Pe Δw , 0} aS = Ds { A( Pe Δs ) +  Pe Δs , 0} (4.276) (4.277) (4.278) (4.279) (4.280) (4.281) (4.282) aP = aE + aW + aN + aS + aT + aB + a − S P Δx ΔyΔz Chapter 4 External Convective Heat and Mass Transfer 387 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 0 aP = ρΔx ΔyΔz Δt (4.283) The expressions for conductance at the faces of the control volume are De = Γ e Δy Δ z Γ Δy Δz Γ Δx Δz , Dw = w , Dn = n (δ x )e (δ x ) w (δ y)n Γ Δx Δ z Γ Δx Δy Γ Δx Δy Ds = s , Dt = t , Db = b (δ y) s (δ z )t (δ z )b (4.284) and the flow rates are: Fe = De Pe Δe = ( ρ u)e ΔyΔz , Fw = Dw Pe Δw = ( ρ u) w ΔyΔz (4.285) Fn = Dn Pe Δn = ( ρ v)n Δx Δz , Fs = Ds Pe Δs = ( ρ v) s Δx Δz Ft = Dt Pe Δt = ( ρ w)t Δx Δy, Fb = Db Pe Δb = ( ρ w)b Δx Δy The Different discretization schemes for the above three-dimensional problem can be obtained by using different expressions for A(|PeΔ|) from Table 4.3. In addition to the six first order discretization schemes described above, some researchers have used higher order schemes such as second order upwind (Leonard et al., 1981) and QUICK (Quadratic Upwind Interpolation of Convective Kinetics; Leonard, 1979) schemes to overcome the false diffusion problem, which is referred to as error caused by using the discretization scheme with accuracy less than the second order (Patankar, 1980). The error due to false diffusion could potentially be severe for (1) transient problems, (2) multidimensional steady-state problems, or (3) problems with non-constant source terms (Tao, 2001). While the accuracies of these higher order schemes are better than the first order schemes, their computational time is much greater than that of the first order schemes. 4.8.3 Numerical Solution of Flow Field Special Difficulties The discussions in the previous subsection are for the solution of a general variable φ with a known flow field. We now turn our attention to the solution of the flow field. For a two-dimensional incompressible flow problem without a body force, the continuity and momentum equations in the Cartesian coordinate system are ∂u ∂v + =0 ∂x ∂y (4.286) (4.287) (4.288) ρu ρu ∂u ∂u ∂  ∂u  ∂  ∂u  ∂p + ρv = − μ + μ ∂x ∂y ∂x  ∂x  ∂y  ∂y  ∂x ∂v ∂v ∂  ∂v  ∂  ∂v  ∂p + ρv = − μ + μ ∂x ∂y ∂x  ∂x  ∂y  ∂y  ∂y 388 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press It is observed that the momentum equations in each direction could be expressed in the same format as eq. (4.200) by replacing the general variable φ by the components of velocity in that particular (x-, or y-) direction. The source term for the momentum equations in the x- and y- directions can be expressed as Su = −∂p / ∂x and S v = −∂p / ∂y, respectively. It seems that we can directly apply the discretization schemes discussed in the preceding subsection in order to obtain the solution of the flow field. The only additional work necessary is to include the pressure gradient in each momentum equation. With the aid of Fig. 4.17, the pressure gradient in the x-direction can be expressed as pe − pw  ∂p   = ∂x  P Δx  (4.289) If central difference is employed, the pressures at the faces of the control volume become pe = p + pW p E + pP , pw = P 2 2 Substituting the above expressions into eq. (4.289), the pressure gradient in the x-direction becomes pE − pW  ∂p   = ∂x  P 2 Δx  (4.290) Similarly, the pressure gradient in the y-direction can be expressed as  ∂p  pN − pS  = ∂y  P 2Δy  (4.291) It can be seen that the pressure gradient at point P is related to the pressures of the neighbor grid points and is not related to its own pressure. Thus, the 50 10 50 10 50 10 50 100 30 100 30 100 30 100 50 10 50 10 50 10 50 100 30 100 30 100 30 100 50 10 50 10 50 10 50 100 30 100 30 100 30 100 50 10 50 10 50 10 50 Figure 4.22 Checkerboard pressure field. Chapter 4 External Convective Heat and Mass Transfer 389 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press pressure gradient is obtained by using the pressures of two alternative points, so that the accuracy of the pressure gradient for grid size Δx (or Δy) is equivalent to that obtained for coarse grid size 2Δx (or 2Δy). Another consequence of eqs. (4.290) and (4.291) is that if we have a pressure distribution as shown in Fig. 4.22 – referred to as a checkerboard pressure field – the discretization scheme represented by eqs. (4.290) and (4.291) will obtain ∂p / ∂x = 0 and ∂p / ∂y = 0 throughout the computational domain, i.e., eqs. (4.290) and (4.291) cannot recognize the difference between a checkerboard pressure field and a uniform pressure field. In other words, if a checkerboard pressure field is supposed to be the real pressure field for any reason, the discretization scheme based on eqs. (4.290) and (4.291) can neither detect nor eliminate such an unrealistic effect. While it is straightforward to use the momentum equations (4.287) and (4.288) to solve for the velocity components in the x- and y-directions, the pressure only appears as a source term in the momentum equations and there is not a separate equation for pressure itself. Since the continuity equation, (4.286), is not used, it would be logical to use the continuity equation to get the pressure field. Although a correct pressure field will yield a velocity field that satisfies the continuity equation, an algorithm must be developed to obtain the correct pressure field. The checkerboard pressure field problem and the lack of an equation for pressure are associated with the coupling between the pressure and velocity. An unrealistic (e.g., checkerboard) pressure field can be obtained if such a coupling relationship is not reflected in the algorithm (e.g., using central difference scheme to discretize the pressure gradient in the momentum equations). Staggered grid To overcome the checkerboard pressure field and develop an effective algorithm for the pressure field, one can use a staggered grid system (Patankar and Spalding, 1972; Patankar, 1980) that stores the pressure and all other variables on the main grid but calculates the velocity at the face of the control volume (see Fig. 4.23). The control volume for the velocity component in the x-direction, u, N uw W P vs vn ue E W uw P vs N N vn ue E W uw P S vs n vn ue E S S (a) (b) (c) Figure 4.23 Staggered grid: (a) control volume for all other variables, (b) control volume for u, and (c) control volume for v. 390 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press is staggered from the control volume for all other variables to the right direction by half a grid. Similarly, the control volume for v is staggered up by a half grid. The discretized equations for u and v will be obtained by integrating the momentum equation in their control volumes shown in Fig. 4.23 (b) and (c), respectively. The discretization equations in the staggered grid system for all variables except velocity are the same as those presented in the preceding subsection. Since velocity is defined at the face of the control volume, the grid Peclet numbers can be directly calculated from the velocity component, i.e., no assumption on the velocity profile between the grid points is needed. For the momentum equation in the x-direction, we will need to integrate eq. (4.287) for the control volumes of u shown in Fig. 4.23 (b). The result can be expressed as ae ue =  anb unb + b + Ae ( pP − pE ) (4.292) where the first term on the right hand side represents the summation of all the neighbor points of e. The effect of pressure has been separated from the source term, and the pressures at P and E were used to calculate the velocity ue. For a two-dimensional problem, the area on which the pressure difference acts on is Ae = Δy . Similarly, the discretization equation for v can be obtained by integrating eq. (4.288) for the control volume of vn shown in Fig. 4.23(c), i.e., (4.293) an vn =  anb vnb + b + An ( pP − pN ) where the area on which pressure difference acts on is An = Δx for a twodimensional problem. For a three-dimensional problem, eqs. (4.292) and (4.293) are still valid except the neighbor points from the top and bottom should be included in the first term on the right-hand side. The areas on which the pressure difference act on should be modified to Ae = ΔyΔz and An = Δx Δz . In addition to eqs. (4.292) and (4.293), another equation for the velocity component in the z-direction, wt, is also needed, i.e., at wt =  anb wnb + b + At ( pP − pT ) (4.294) which is obtained by integrating the momentum equation in the z-direction for the control volume of wt that are staggered to the positive z-direction by a half grid. The area on which the pressure difference, pP − pT , acts on is At = Δx Δy . Pressure Correction Equation The ability (of one) to solve eqs. (4.292) – (4.294) to obtain the velocity field relies on the accuracy of the pressure field, which is unknown a priori. If we assume that the pressure field is p* – which can be obtained from the previous time step, iteration, or guessing – the corresponding velocity field can be obtained by solving eqs. (4.292) – (4.294). * * ae ue =  anb unb + b + Ae ( p* − p* ) (4.295) P E Chapter 4 External Convective Heat and Mass Transfer 391 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press * * an vn =  anb vnb + b + An ( p* − p* ) P N * * at wt* =  anb wnb + b + At ( p* − pT ) P (4.296) (4.297) Since the velocity field obtained from eqs. (4.295) – (4.297) is based on a guessed pressure field, it may not satisfy the continuity equation. An algorithm based on pressure correction (Patankar, 1980) can be used to obtain the correct pressure field and consequently, the velocity field. If the correct pressure can be obtained by p = p* + p ′ (4.298) where p′ is pressure correction, the new velocity field based on the corrected pressure field will be u = u* + u′, v = v* + v′, w = w* + w′ (4.299) where u′, v′, and w′ are velocity corrections. The corrected velocity field satisfies eqs. (4.292) – (4.294), i.e., * * ′ ′ ′ ′ ae (ue + ue ) =  anb (unb + unb ) + b + Ae [( p* + pP ) − ( p* + pE )] P E * * ′ ′ ′ ′ an (vn + vn ) =  anb (vnb + vnb ) + b + An [( p* + pP ) − ( p* + pN )] P N * * ′ ′ ′ at (wt* + wt′ ) =  anb (wnb + wnb ) + b + At [( p* + pP ) − ( pT + pT )] P Subtracting eqs. (4.295) – (4.297) from the above three equations yields the following equations for the velocity corrections: ′ ′ ′ ′ ae ue =  anb unb + Ae ( pP − pE ) (4.300) ′ ′ ′ an vn =  anb vnb + An ( pP − p′ ) N ′ ′ ′ at wt′ =  anb wnb + At ( pP − pT ) (4.301) (4.302) It can be seen from eqs. (4.300) – (4.302) that the correction of the velocity can be from (1) the difference on pressure correction at two neighbor points in the main grid (the second terms), and (2) velocity correction from all neighbor points (the first terms). The solution of the velocity correction from eqs. (4.300) – (4.302) will be very complicated because the velocity correction at neighbor points will be related to pressure corrections at more grid points, and eventually velocity correction at any point will be related to the pressure corrections in the entire computational domain. Assuming the contribution of the latter is negligible, which is equivalent to assuming anb = 0 in eqs. (4.300) – (4.302), the velocity correction becomes ′ ′ ′ ′ ′ ′ ′ ue = de ( pP − pE ), vn = dn ( pP − p′ ), wt′ = dt ( pP − pT ) (4.303) N where de = Ae / ae , dn = An / an , dt = At / at (4.304) Substituting eq. (4.303) into eq. (4.299), we have 392 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press * ′ ′ ue = ue + de ( pP − pE ) * ′ vn = vn + dn ( pP − p′ ) N (4.305) (4.306) (4.307) ′ ′ wt = wt* + dt ( pP − pT ) which can be used to obtain the corrected velocity field from the pressure correction. The assumption of neglecting the first terms in eqs. (4.300) – (4.302) significantly simplifies the equations in order to obtain the corrected velocity field. If, however, this simplification results in unacceptable error in the final velocity field, it may not be used. As the iteration progresses, v∗ , gradually approaches the true velocity, the velocity correction at each point approaches zero. After the converged velocity is obtained, the velocity correction from all the neighbor points will be zero. Therefore, neglecting the first terms in eqs. (4.300) – (4.302) will not affect the ultimate velocity field. It will, however, affect how quickly a converged solution can be obtained. We now turn our attention to the algebraic equation for the pressure correction. Since the starred velocity field (u*, v*, and w*) does not satisfy the continuity equation, we hope that the corrected velocity will satisfy the continuity equation. For a three-dimensional flow, the continuity equation is ∂ρ ∂ ( ρ u) ∂ ( ρ v) ∂ ( ρ w) (4.308) + + + =0 ∂t ∂x ∂y ∂z Integrating eq. (4.308) over the main control volume [see Fig. 4.23(a) for the two-dimensional projection of the three-dimensional control volume] and using a fully-implicit scheme, we obtain 0 ρP − ρP Δx ΔyΔz + [( ρ u)e − ( ρ u) w ]ΔyΔz Δt +[( ρ v)n − ( ρ v) s ]Δx Δz + [( ρ w)t − ( ρ w)b ]Δx Δy = 0 Substituting eqs. (4.305) – (4.307) into the above equation, the following algebraic equation for pressure correction is obtained: ′ ′ ′ ′ ′ ′ ′ aP pP = aE pE + aW pW + aN pN + aS pS + aT pT + aB pB + b (4.309) where aE = ρ e de ΔyΔz , aW = ρ w dw ΔyΔz aN = ρn dn Δz Δx , aS = ρ s ds Δz Δx aT = ρt dt Δx Δy, aB = ρb db Δx Δy aP = aE + aW + aN + aS + aT + aB b= 0 ρP − ρP (4.310) (4.311) (4.312) (4.313) Δt Δx ΔyΔz + [( ρ u* ) w − ( ρ u* )e ]ΔyΔz +[( ρ v* ) s − ( ρ v* )n ]Δx Δz + [( ρ w* )b − ( ρ w* )t ]Δx Δy (4.314) Chapter 4 External Convective Heat and Mass Transfer 393 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press N Boundary vn ue (given) uw W P vs S E Figure 4.24 Control volume for continuity equation at boundary. If the starred velocity terms (u*, v*, and w*) already satisfy the continuity equation, we have b = 0 from eq. (4.314). Therefore, b defined in eq. (4.314) can be viewed as the negative value of the “residual” mass for grid point P. If the continuity equation is satisfied, such residual mass should be zero. Practically, the maximum absolute value of b in the entire computational domain as well as the absolute value of the summation of b in the entire domain are used as criteria for the convergence. In order to apply eq. (4.309) to obtain the pressure correction, the boundary conditions for pressure correction have to be specified. If the pressure at the boundary is given, no correction for the pressure at the boundary is needed so p′ = 0 at the boundary. If the boundary velocity is given (see Fig. 4.24), the velocity correction for the boundary velocity (ue in Fig. 4.24) will be zero. It follows from eq. (4.303) that de for the control volume shown in Fig. 4.24 has to be zero. Therefore, one can set aE = 0 [see eq. (4.310)] if ue in Fig. 4.24 is given. The SIMPLE Algorithm The above approach for the solution of the incompressible flow field was named SIMPLE, which stands for Semi-Implicit Method for Pressure-Linked Equations. It was originally proposed by Patankar and Spalding (1972) and summarized in Patankar (1980). This algorithm is a semi-implicit method because the first terms on the right-hand side of eqs. (4.300) – (4.302) are neglected. If these terms are retained, we will have to solve for the velocity ′′ corrections ( ue , vn , and wt′ ) for the entire flow field simultaneously and the algorithm will become fully-implicit. The procedure for SIMPLE algorithm can be summarized as following: 1. Guess a pressure field, p*. 2. Obtain the starred velocity field, (u*, v*, and w*) from eqs. (4.295) – (4.297). 3. Solve the pressure correction equation (4.309) to get p′ . 4. Obtain a new pressure field from eq. (4.298). 5. Improve the velocity field using eqs. (4.305) – (4.307). 394 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 6. Obtain solutions for other φ ’s from eq. (4.274). If the flow field is not affected by a particular φ , it should be solved after a converged solution for flow field is obtained. 7. Treat the corrected pressure in Step 4 as a new starred pressure and go back to step 2. The iteration procedure is repeated until a converged solution is obtained. Since the invention of the SIMPLE algorithm in 1972, it has evolved into the classic approach for computational fluid dynamics and heat transfer. Although there are several improved versions in the literature, SIMPLE still remains one of basic the most powerful algorithms. The objective of this section is to introduce the readers to the world of computational fluid dynamics and heat transfer without the overwhelming numerical details. Additional information can be found in the seminal books of Patankar (1980) and Tao (2001), as well as the Handbook of Numerical Heat Transfer (Minkowycz et al., 2006). 4.8.4 Numerical Simulation of Interfaces and Free Surfaces Many engineering applications involve interfacial phenomena because of the high heat transfer rates that can be achieved. There are many complexities that need to be addressed when modeling an interface, since an interface is generally irregular, involves mass transfer, is three-dimensional, and is not at a fixed location. There are two different approaches when considering an interface: a continuum and a noncontinuum approach. In many problems, the scale of the entire computational domain is much larger than the scale of the interfacial thickness; therefore the interface can be considered as a planar surface and the continuum approach will give highly accurate results. With the development of nanotechnology, the scales of systems are becoming very small. Therefore, neglecting the thickness of the interface is no longer a valid assumption, and a noncontinuum approach is needed. This approach can be achieved through molecular dynamic simulations. A detailed discussion on tracking an interface using the various approaches is given by Faghri and Zhang (2006). Interface Tracking Techniques There are many engineering problems that involve separated phases in which the interface between phases is clearly defined, but is not at a fixed location. The interface can be assumed to be a planar surface in these situations, and therefore, a continuum approach is valid. In these problems, it is necessary to solve for both the phases as well as the interfacial location. These problems can be solved numerically on an Eulerian or Lagrangian mesh. An Eulerian mesh is stationary and defined prior to the start of a solution. When using an Eulerian mesh, the interface is tracked by solving an additional scalar equation. In the Lagrangian Chapter 4 External Convective Heat and Mass Transfer 395 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Actual Interface Interface Modeled as a Boundary (a) Actual interface Interface approximated by piecewise interpolation from cell volume fraction (c) Phase interface fitted grid Interface captured by a separate interfacial grid (b) Volume of Fluid model (d) Front tracking/capturing Figure 4.25 Interfacial representations for different interface tracking techniques. approach, a boundary of the mesh is aligned with the interface, and this boundary moves with the interface. When thinking of a multiphase system from a continuum approach, in the bulk region a phase is continuous and is discontinuous at an interface between different phases. In general, the interface is free to deform based on the nature of the flow. Therefore, it is difficult to efficiently capture an interface between phases with just one model. Consequently, there have been strong efforts to track an interface based on several different techniques, each with its own pros and cons. An actual interface, represented by the different numerical techniques, is presented in Fig. 4.25. The techniques are as follows: 396 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1. Stationary Grid Approach: standard CFD modeling in cells that contain only a single phase; special consideration taken in cells in the vicinity of an interface. 2. Lagrangian Techniques: grid and fluid move together, interface is directly captured. 3. Phase Interface Fitted Grid: interface directly tracked as boundary, rest of grid moves as a function of interfacial movement, “semi-Lagrangian”. 4. Front Tracking Methods: stationary grid is used and modified near the interface, so that the grid is aligned with the interface; combination of 2 and 3. The above methods are discussed in detail in Faghri and Zhang (2006). The stationary grid approach, which is used widely for solving single phase problems, is presented below. Stationary Grid Approach The stationary grid approach is based on a stationary grid where the fluid interface is captured directly. The first of such approaches is the marker-and-cell approach originated by Harlow and Welch (1965). In this approach, massless particles are introduced into the flow field, and the locations are projected from their interpolated velocities. Cells with a particle are considered to have one phase in them, and cells without a particle do not have that phase in them. An interface is considered where cells with particles are neighbored with cells without particles. For efficiency, this method was extended to only track particles on the surface (Nichols and Hirt, 1975). Further development of this class of models led to the Volume of Fluid method (VOF) by Hirt and Nichols (1981) in which they used a donor-acceptor method to effectively eliminate numerical diffusion at an interface. The VOF method is one of the most widely used routines to solve a free surface problem on an Eulerian mesh. In this method, there is one velocity, pressure, and temperature field, which is shared by all of the phases modeled. The interface between the phases is tracked by the volume fraction of phase k, ε k . The volume fraction equation is the continuity equation of phase k divided by the density of that phase. 1 ∂  ( ε k ρk ) + ∇ ⋅ ( ε k ρk V ) = ρk  ∂t    jk m ′′′  j =1( j ≠ k )    Π (4.315) where Π is the number of phases. When the volume fraction is between 0 and 1 in a computational cell, the cell is considered to be an interfacial cell. For a two-phase system, the phases are k and j. When the volume fraction is 1, that cell is occupied by only phase k, and when it is 0, that cell is occupied by only phase j. The sum of the volume fractions is unity, Chapter 4 External Convective Heat and Mass Transfer 397 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ε k =1 Π k =1 (4.316) Therefore, ( Π − 1 ) volume fraction equations need to be solved. The fluid properties, such as density, viscosity, and thermal conductivity, are calculated by their volume weighted average. Φ eff =  ε k Φ k k =1 Π (4.317) The overall continuity equation in conjunction with the VOF method is: ∂ ( ρeff ) + ∇ ⋅ ( ρeff V ) = 0 ∂t (4.318) The continuity equation is the same as eq. (3.38) except that it uses the effective density and the velocity field is shared by both phases; therefore, subscript k is dropped. The momentum equation is Π Π ∂  ′ ( ρeff V ) + ∇ ⋅ ( ρeff VV ) = ∇ ⋅τ eff + ρeff X +   m′′′Vjk ∂t k =1 j =1( j ≠ k ) (4.319) The momentum equation is also the same as eq. (3.52), with the same exceptions as the continuity equation and the momentum transfer due to phase change. In the energy equation, the enthalpy is mass averaged instead of volume averaged. heff = 1 ρeff ε k =1 Π k ρk hk (4.320) Therefore the energy equation, neglecting pressure effects and viscous dissipation, is: Π Π ∂  jk  ′′′ ( ρeff heff ) + ∇ ⋅ ( ρeff Vheff ) = −∇ ⋅ q′′ +   m′′′ hk + qeff eff ∂t k =1 j =1( j ≠ k ) (4.321) The energy equation has the same form as eq. (3.79) except effective properties are incorporated and a latent heat term due to phase change is added. The continuity, momentum, and energy equations can be solved by standard solution procedures, such as the SIMPLE class of algorithms in Patankar (1980) and Patankar (1991). However, if the volume fraction equation is solved using a standard implicit or explicit scheme, the interface will quickly lose resolution due to numerical diffusion. The numerical diffusion will lead to inaccurate solutions and/or a solution that will not converge. Therefore, special consideration must be taken on interfacial cells to construct the interface so that the advection transport of fluid is representative of the physical problem. An actual interface between two phases, as well as the corresponding interface represented by two special interpolation schemes, is presented in Fig. 4.26. As noted before, the first of such methods is the donor-acceptor scheme proposed by Hirt and Nichols (1981). If a cell is an interfacial cell, 0 < ε k < 1 , the fluid will be rearranged in the cell to be on one face, as shown in Fig. 4.26(b). The face on which the fluid will be rearranged will depend on the normal direction of the interface. The normal direction of the interface, with respect to phase k, can be calculated by the gradient of the volume fraction: 398 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press n (a) (b) (c) Figure 4.26 (a) An actual interface between two phases, (b) an interfacial representation using the Donor-Acceptor scheme, and (c) a piecewise linear reconstruction scheme with the VOF method. nk = − ∇ε k ∇ε k (4.322) The component in which the normal is the greatest will occur where the fluid is perpendicular to the interface while the interface is either horizontal or vertical. Once this is done, one cell is designated as a donor, while its neighbor is designated as an acceptor. The amount of fluid leaving the donor cell is exactly equal to the amount of fluid entering the acceptor through each computational face. Also, the maximum amount of fluid that can leave a cell is equal to either the amount of fluid in that cell or the amount that would make another cell fill with fluid. A more refined interface interpolation scheme was developed by Youngs (1982), in which the interface was approximated as piecewise linear in each cell. The normal of the reconstructed interface in each cell is the same as the normal calculated in eq. (4.322), as shown in Fig. 4.26(c). The advection of fluid through each face is calculated in a manner similar to the donor-acceptor scheme. Both the donor-acceptor scheme and piecewise linear interpolation can only be run in transient simulations, and the time step must be kept small enough so that the fluid near the interface will only advance by one cell at a time. Also, the Chapter 4 External Convective Heat and Mass Transfer 399 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press resolution of the interface is limited to the grid spacing of the computational mesh. Any surface waves that are smaller than the spacing of the mesh will be smoothed out, as seen in Fig. 4.26 on the left side of the interface. Surface tension effects are important in many free surface problems. The surface tension effects can be applied as a body force in the momentum equations. This method is called the continuum surface force (CSF) model and was proposed by Brackbill et al. (1992). The curvature is defined as the divergence of the surface normal vector. K k = ∇ ⋅ nk (4.323) The volume force due to surface tension, Fσ , is Π k −1 ε j ρ j K k ∇ε j + ε k ρ k K j ∇ε k Fσ =  −σ jk (4.324) 1 k =1 j =1 ( ρ j + ρk ) 2 If only two phases are present, the body force reduces to ρ eff K k ∇ε j Fσ = σ jk 1 ( ρk + ρ j ) 2 (4.325) It can be seen that the body force is proportional to the cell density. It is important to note that when surface tension forces are large compared to other flow characteristics, numerical inaccuracies of the surface tension forces in the CSF model create artificial currents called parasitic currents. Much work has been done to eliminate parasitic currents, such as the second gradient method proposed by Jamet et al. (2002). Despite the adverse effects of parasitic currents, the VOF method has been widely used and can give reasonably accurate results for a wide range of applications. It is robust in its handling of free surface problems with large interface distortion and can easily handle problems in which the free surface breaks apart, such as droplet formation. One drawback of the VOF method, however, is that the interface resolution is limited to the grid spacing. Therefore, a refined mesh is needed anywhere the interface is going to travel. This refinement can lead to many computational cells in regions of the mesh resided in by the interface for a short period of time, which will increase the total computational time of the solution. Therefore, advanced remeshing algorithms are needed for problems of this type. 4.9 Application of Computational Methods The computational methods presented in the previous sections provide the solution techniques for solving the governing equations and boundary conditions. However, it is equally important to mathematically set up the physical problem before applying the numerical techniques. In many practical problems for external flow one needs to deal with conjugate effects, free surface, and nonconventional geometry including multiple domain regions with non-uniform 400 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press boundary conditions. In order to show these effects, and the approach to deal with them, an application is presented below that includes all these nonconventional effects. The physical problem that will be presented in this section is heat and mass transfer from a controlled liquid impinging jet over a rotating disk, Figure 4.27, including the conjugate wall effect for both heating with and without evaporation from the free surface. Figure 4.27 Controlled liquid impinging jet onto a rotating disk (Rice et. al. 2005). Figure 4.28 Computational domain for a controlled liquid impinging jet onto a rotating disk (Rice et. al. 2005). Chapter 4 External Convective Heat and Mass Transfer 401 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press First we need to develop the conservative governing equations and boundary conditions for fluid mechanics, heat transfer, and volume of fluid (VOF) for a controlled liquid impinging jet on a rotating disk with a free surface. Figure 4.28 shows a computational domain for the three regions: disk, liquid, and vapor, with an extended domain to include the edge effects. The basic assumptions are two dimensional, constant properties, incompressible, and laminar flow. For evaporation, the effect of mass transfer is neglected. It is also assumed that there is constant heat flux along the heated section. Flow enters the disk between two circular plates, one being the collar, and the other being the disk. The space between the collar and the disk is δin. After the flow leaves the entrance region between the collar and the disk at rhin, the flow turns from an internal fully developed flow to a free surface flow. The heater provides a constant heat flux at the outside wall, from the bottom of the aluminum disk. The heat is conducted through the disk to the fluid. The basic assumptions of the problem at hand are that the flow field is incompressible, and the fluid is considered to be in the laminar flow regime while all of the fluid properties are constant. When evaporation is being modeled, no mass transfer is modeled and the effects are only considered in the energy equation, because the evaporation rate of the liquid is much less than the inlet mass flow rate of the liquid. The Navier-Stokes equations are solved to compute the fluid flow. The continuity and momentum equations are as follows: ∇⋅V = 0 (4.326) ρ DV = −∇p + ∇ ⋅ ( μ∇V ) + ρ g + F Dt (4.327) Since the flow field is assumed to be independent of the temperature field, a steady-state energy equation is solved in the fluid region once the flow field has been resolved. ∇ ⋅ ( Vρ h ) = ∇ ⋅ ( k∇T ) (4.328) In the solid region, only conduction is solved. ∇ 2T = 0 (4.329) The free surface is tracked by the Volume of Fluid (VOF) method, where the volume fraction of the fluid, ε, is tracked through each computational cell. The VOF equation is: ∂ε + V ⋅∇ε = 0 (4.330) ∂t The interface between fluids is represented by a piecewise linear approach, similar to the work of Youngs (1982), in order to greatly limit numerical diffusion of the interface. Surface tension effects are modeled in the numerical simulations. The surface tension forces are represented by the F term in eq. (4.327). A continuum surface force method proposed by Brackbill et al. (1992) is used to model surface tension. ρ K ∇ε F =σ (4.331) 0.5 ( ρ  + ρ v ) 402 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The curvature, K, is defined as: K =− ∇ 2ε ∇ε (4.332) The fluid properties are calculated by the volume weighted average: φ = εφ + (1 − ε ) φv (4.333) The general fluid property, φ , represents either density, viscosity or thermal conductivity. The enthalpy is calculated using a mass-weighted average instead of a volume-weighted average. ερ c + (1 − ε ) ρ v cv  ( T − Tref )   (4.334) h= ρ The boundary conditions are as follows: At the inlet ( r = rin ) ( 0 < x < δ in ) , v = vin , T = Tin (4.335) (4.336) ( −dd < x < 0 ) , At the collar ( x = δ in , rin < r < rhin ) ∂T =0 ∂r ∂T =0 ∂x u = v = 0, w = ωr , (4.337) At the disk surface ( x = 0, rin < r < rout ) u = v = 0, w = ωr , k ∂T ∂x ∂T ∂x ∂T ∂x =k  ∂T ∂x (4.338) S At the disk bottom ( x = −dd ) ( rin < r < rhin ) : k k =0 S (4.339) (4.340) (4.341) ( rhin < r < rhout ) : ( rhout ′′ = qw S < r < rout ) : k ∂T ∂x =0 S At the outer boundaries ( r = rhin , δ in < x < 100δ in ) , ( rhin < r < rout , x = 100δ in ) , ( r = rout , δ in < x < 100δ in ) p = pref (4.342) if ( V ⋅ n > 0 ) then ∂T  = 0, else T = Tref ∂n (4.343) At the liquid-vapor interface p − pv − μ ∂V ⋅ n ∂V ⋅ n + μv =σK ∂n  ∂n v (4.344) Chapter 4 External Convective Heat and Mass Transfer 403 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press μ ∂V ⋅ t ∂V ⋅ t = μv ∂n  ∂n k ∂T ∂n = kv  (4.345) v heating: ∂T ∂n (4.346) v evaporation: (4.347) Rice et al. (2005) solved this problem numerically with the result presented below. Air was used as the vapor for all of the simulations. The viscosity and thermal conductivity of air were lessened by an order of magnitude for selective cases, and were found not to change the results. The boundary condition at the interface is therefore essentially a zero shear stress for all cases, an insulated boundary for heating cases, and Tsat for evaporation. As was noted before, the flow field is considered to be independent of the thermal field; therefore the flow field is first solved as a time-dependent solution, until a steady state solution has been reached. Once a steady-state solution has been reached, the energy equation is solved for both the fluid and solid regions, as a steady-state solution. The criterion used to determine a steady-state solution was a mass flow rate of liquid leaving the computational domain within 0.05% of the mass flow rate of liquid entering the domain for greater than 0.05 seconds. The flow field was solved using the following procedure: 1. Solve momentum equations 2. Solve continuity (pressure correction) equation, and update pressure and face mass flow rate 3. Repeat steps 1 and 2 until converged, for each time step A co-located finite volume computational scheme as described in Section 4.8, where both the flow-field variables and the pressure are stored in the cell centers, was used to solve the governing equations. The pressure was discretized in a manner similar to a staggered-grid scheme, while the pressure and velocity were coupled using the SIMPLE algorithm described in the previous section. All of the convective terms in the governing equations were discretized using a second-order upwind scheme. The grid used for the numerical simulations consisted of rectangular-shaped cells, which were produced in four axial layers in the liquid regions, and one axial layer in the solid region. In the first layer of the fluid region, spanning from the disk surface to δin, there were 25 cell rows, with a growth rate of 1.01. The second layer, spanning from δin to 3δin, consisted of 30 cell rows with a growth rate of 1.02. The third layer, spanning from 3δin to 10δin, consisted of 25 cell rows, and had a growth rate of 1.05. In the final layer, spanning from 10δin to 100δin, there were 20 cell rows, with a growth rate of 1.2. In the solid region, there were 15 evenly spaced cell rows. In the radial direction, there were 15 evenly spaced cell columns between rin and rhin, and 75 evenly spaced cell columns between rhin and rout. The grid spaned such a great distance in the axial T = Tsat 404 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press direction (100δin), because convergence is greatly increased with the added cells, so the computations ran more rapidly. 2.5 2 δ / δin 1.5 1 0.5 0 1 2 100 rpm (present) 200 rpm (present) r / rhin 3 4 Figure 4.29 Film thickness vs. radial distance at an inlet temperature of 20ºC at δin = 0.254 mm and a liquid flow rate of 7 RPM (Rice et. al. 2005) 500 200 rpm 100 rpm 50 rpm 400 Nu 300 200 100 1 1.5 2 r / rhin 2.5 3 3.5 Figure 4.30 Local Nusselt Numbers vs. radial distance for δin = 0.204 mm, inlet temperature of 40ºC and liquid flow rate of 7 RPM (Rice et. al. 2005). Chapter 4 External Convective Heat and Mass Transfer 405 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Typical results for local film thickness δ and Nussult number, Nu = q′′rin / [k (Tw − Tin )] , as a function of radial distance are presented in Fig. 4.29 and 4.30 respectively. The heating boundary condition was used for all cases with an inlet temperature of 20ºC to 40ºC. The evaporation boundary condition was used for all cases with an outlet temperature of 100ºC. 4.10 Analogies and Differences in Different Transport Phenomena Most early work in theoretically predicting the heat and/or mass transfer in both laminar and turbulent flow cases was done using the analogy between momentum, heat, and mass and by predicting the approximate results for the heat and/or mass transfer coefficient from the momentum transfer or skin friction coefficient. Clearly there are severe limitations in using this simple approach; however, it is beneficial to understand the advantages and similarities for physical and mathematical modeling as well as the constraints involving this approach. We present this analogy for the classic problem of heat and mass transfer over a flat plate in this section. Its application to more complicated geometries and boundary conditions as well as turbulent flow is not proven and caution should be taken in applying this approach to other cases. As presented before, a flat plate at constant wall temperature Tw is exposed to a free stream of constant velocity U∞, temperature T∞ and mass fraction ω1,∞. Due to binary diffusion it can be presented with the following dimensionless conservation boundary layer equations and boundary conditions corresponding to Fig. 4.31. y, v U∞ T∞ ω1, ∞ δc δ δT x, u Tw and ω1, w Figure 4.31 Mass, momentum, and heat transfer in a laminar boundary layer. 406 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Continuity ∂u + ∂v+ + =0 ∂x + ∂y + (4.348) Momentum u+ ∂u + ∂u + 1 ∂ 2u+ + v+ + = Re ∂y + 2 ∂x + ∂y (4.349) Energy u+ ∂θ ∂θ 1 ∂ 2θ + v+ + = + ∂x ∂y Pr Re ∂y +2 ∂φ ∂φ 1 ∂ 2φ + v+ + = + ∂x ∂y Sc Re ∂y +2 y = 0, v = v + + + + w (4.350) Species u+ (4.351) (4.352) (4.353) (4.354) at y + = 0, u + = 0, θ =φ =0 y = ∞, u + = 1, θ = φ = 1 where the dimensionless variables are defined as: ω − ω1, w T − Tw u u+ = , θ= , φ= 1 , ω1,∞ − ω1, w U∞ T∞ − Tw v y x v= , y+ = , x + = , U∞ L L + ν Let’s first consider the analogy between momentum and heat transfer. Equations (4.349) and (4.350) and appropriate boundary conditions are the same if Pr = 1. The solutions for u+ and θ are, therefore, exactly the same if Pr = 1 and one expects to have a relationship between the skin friction coefficient, cf, and heat transfer coefficient, h. cf 2 = Re = U∞L (4.355) τw y =0 = = ρU ∞ 2 ρU ∞ 2 −k ∂T ∂y μ ∂u ∂y ν ∂u + ∂y + y+ = 0 U∞L = ∂u + ∂y + y+ = 0 Re (4.356) L hL ∂θ y=0 Nu = = =+ k k ( Tω − T∞ ) ∂y (4.357) y+ = 0 Since θ = u+ for Pr = 1, one also concludes that ∂θ ∂u + ∂y + = y+ = 0 ∂y + Nu Re (4.358) y+ = 0 Therefore, combining eqs. (4.356) and (4.357) and using eq. (4.358) gives cf 2 = (4.359) Chapter 4 External Convective Heat and Mass Transfer 407 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press This relationship between the friction coefficient and Nu is referred to as Reynolds analogy and is appropriate when Pr = 1. If Pr ≠ 1, we already concluded that δT = Pr −1/3 for 0.5 ≤ Pr ≤ 10 δ from the similarity solution presented in Section 4.6. Using this information, one can generalize the result of the Reynolds analogy to Pr ≠ 1 by cf 2 = Nu −1/3 Pr Re (4.360) Now, we focus on similarities between the heat and mass transfer for comparison of eqs. (4.350) and (4.351) with their appropriate boundary conditions. It is clear that the solution for differential equations (4.350) and (4.351) for θ and φ are the same if Sc and Pr are interchanged appropriately. We already know the solution for eq. (4.350) for vw = 0 from the similarity solution for 0.5 ≤ Pr ≤ 10 Nu x = 0.332 Re1/ 2 Pr1/3 (4.361) x Therefore, it can also be assumed that the solution of eq. (4.351) for vw = 0 is Sh x = 0.332 Re1/ 2 Sc1/3 (4.362) x Combining eqs. (4.361) and (4.362) gives Nu Pr1/3 = Sh Sc1/3 (4.363) It should be noted that the convective effect due to vertical velocity at the surface was neglected in predicting h and hm in the above analysis and therefore the analogy presented in (4.363) is for very low mass transfer at the wall. The effect is that the contribution of vw in heat and mass flux was primarily due to diffusion. ′′ qw = ρ c p vw ( Tw − T∞ ) − k ∂T ∂y (4.364) y=0  ∂ω  ′′ mw = ρ  ω1,ω vw − D12 1  ∂y     y=0  (4.365) The analogy in momentum, heat, and mass transfer can also be applied to complex and/or coupled transport phenomena problems including phase change and chemical reactions. Obviously, in these circumstances, the simple relationship developed in eq. (4.363) is not applicable. To show the usefulness of this analogy to more complex transport phenomena problems we will apply it to sublimation with a chemical reaction for forced convective flow over a flat plate. During combustion involving a solid fuel, the solid fuel may either burn directly or be sublimated before combustion. In the latter case – which will be discussed in this subsection – gaseous fuel diffuses away from the solid-vapor interface. Meanwhile, the gaseous oxidant diffuses toward the solid-vapor interface. Under the right conditions, the mass flux of vapor fuel and the gaseous 408 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press oxidant meet and the chemical reaction occurs at a certain zone known as the flame. The flame is usually a very thin region with a color dictated by the temperature of combustion. Figure 4.32 shows the physical model of the problem under consideration. The concentration of the fuel is highest at the solid fuel surface, and decreases as the location of the flame is approached. The gaseous fuel diffuses away from the solid fuel surface and meets the oxidant as it flows parallel to the solid fuel surface. Combustion occurs in a thin reaction zone where the temperature is the highest, and the latent heat of sublimation is supplied by combustion. The combustion of solid fuel through sublimation can be modeled as a steady-state boundary layer type flow with sublimation and chemical reaction. To model the problem, the following assumptions are made: 1. The fuel is supplied by sublimation at a steady rate. 2. The Lewis number is unity, so the thermal and concentration boundary layers have the same thickness. 3. The buoyancy force is negligible. The conservation equations for mass, momentum, energy, and species in the boundary layer are ∂ ( ρ u ) ∂ ( ρ v) (4.366) + =0 ∂x ∂y ∂u ∂u ∂  ∂u  +v = ν  ∂x ∂y ∂y  ∂y  ∂ ∂ ∂  ∂T   ′′′ ( ρ c p uT ) + ( ρ c p vT ) = k  + mo hc ,o ∂x ∂y ∂y  ∂y  u ∂ωo  ∂ ∂ ∂  ′′′ ( ρ uωo ) + ( ρ vωo ) =  ρD  − mo ∂x ∂y ∂y  ∂y  (4.367) (4.368) (4.369) U∞ U∞ y,v Figure 4.32 Sublimation with chemical reaction (Kaviany, 2001). Chapter 4 External Convective Heat and Mass Transfer 409 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  ′′′ where m0 is the rate of oxidant consumption (kg/m3·s), hc,0 is the heat released by combustion per unit mass consumption of the oxidant (J/kg), which is different from the combustion heat, and ω0 is the mass fraction of the oxidant in the gaseous mixture. The corresponding boundary conditions of eqs. (4.366) to (4.369) are u → U ∞ , T → T∞ , ωo → ωo ,∞ _ at _ y → ∞ (4.370) (4.371) ρ f where m′′ is the rate of solid fuel sublimation per unit area (kg/m2·s) and ρ is the u = 0, v = , f m ′′ ∂ωo = 0 _ at _ y = 0 ∂y density of the mixture. The shear stress at the solid fuel surface is τw = μ ∂u , y=0 ∂y ∂T , y=0 ∂y (4.372) The heat flux at the solid fuel surface is ′′ qw = −k (4.373) ∂y  Considering the assumption that the Lewis number is unity, i.e. Le = ν / D = 1 , eq. The exact solution of the heat and mass problem described by eqs. (4.366) to (4.369) can be obtained using numerical simulation. It is useful here to introduce the results obtained by Kaviany (2001) using an analogy between momentum, heat, and mass transfer. Multiplying eq. (4.369) by hc,0 and adding the result to eq. (4.368), one obtains ∂ω  ∂ ∂ ∂  ∂T  ρ u(c p T + ωo hc ,o )  +  ρ v(c p T + ωo hc ,o )  = + ρ Dhc , o o  (4.374) k     ∂x ∂y ∂y  ∂y (4.374) can be rewritten as ∂ ∂  ρ u(c pT + ωo hc ,o )  +   ∂y  ρ v(c p T + ωo hc , o )    ∂x  ∂ ∂ =  ρα ∂y (c p T + ωo hc , o )  ∂y   (4.375) which can be viewed as an energy equation with quantity ( c pT + ω0 hc ,0 ) as a dependent variable. Since ∂ωo / ∂y = 0 at y = 0, i.e., the solid fuel surface is not permeable for the oxidant, eq. (4.373) can be rewritten as ′′ qw = − ρα ∂ (c p T + ωo hc ,o ) , y = 0 ∂y (4.376) Analogy between surface shear stress and the surface energy flux yields τ ′′ qw = w  (c p T + ωo hc ,o ) w − (c p T + ωo hc ,o )∞   U∞  (4.377) τw  c p (Tw − T∞ ) + hc ,o (ωo, w − ωo, ∞ )   =  U∞ 410 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The energy balance at the surface of the solid fuel is (4.378) where the two terms on the right-hand side of eq. (4.378) represent the latent heat of sublimation, and the sensible heat required to raise the surface temperature of the solid fuel to sublimation temperature and heat loss to the solid fuel. Combining eqs. (4.377) and (4.378) yields the rate of sublimation on the solid fuel surface τ f m ′′ = Z w (4.379) U∞ f ′′ ′′ −qw = m′′ hsv + q where Z is the transfer driving force or transfer number defined as c p (T∞ − Tw ) + hc ,o (ωo ,∞ − ωo , w ) Z= ′′  f hsv + q / m′′ (4.380) Using the friction coefficient gives τw cf = 2 ρU ∞ / 2 eq. (4.379) therefore becomes f m ′′ = cf (4.381) 2 = ρU ∞ Z cf (4.382) The surface blowing velocity of the gaseous fuel is then (4.383) U∞ Z ρ 2 where the friction coefficient, cf, can be obtained from the solution of the boundary layer flow over a flat plate with blowing on the surface. The similarity solution of the boundary layer flow problem exists only if the blowing velocity satisfies vw ∝ x −1/ 2 . In this case, one can define a blowing parameter as ( ρ v) w 1/ 2 Re x (4.384) B= ( ρ u) ∞ Combining eqs. (4.383) and (4.384) yields vw = B= Z 1/ 2 Re x c f 2 f m′′ (4.385) Glassman (1987) recommended an empirical form of eq. (4.385) based on numerical and experimental results: B= ln(1 + Z ) 2.6 Z 0.15 (4.386) Example 4.4: Air with a temperature of 27 °C flows at 1 m/s over a 1-m long solid fuel surface with a temperature of 727 °C. The concentration of the oxidant at the solid fuel surface is 0.1, and the heat released per unit mass of the oxidant consumed is 12000 kJ/kg. The latent heat of sublimation for the solid fuel is 1500 kJ/kg. Neglect the sensible heat required to raise the surface temperature of the solid fuel to sublimation temperature, and Chapter 4 External Convective Heat and Mass Transfer 411 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press heat loss to the solid fuel. Estimate the average blowing velocity due to sublimation on the fuel surface. Solution: The mass fractions of the oxygen at the solid fuel surface and in the incoming air are, respectively, ω0, w = 0.1 and ωo,∞ = 0.21 . The specific heat of gas, approximately taken as the specific heat of air at Tave=(Tw+T∞)/2=377 °C, is cp=1.063 kJ/kg-K. The combustion heat per unit oxidant consumed is hc ,o = 12000kJ / kg . The latent heat of sublimation is hsv = 1500kJ / kg . The density at the wall and the incoming temperatures are respectively ρ w = 0.3482 kg / m3 and ρ∞ = 1.1614 kg / m3 . The viscosity at Tave is ν =60.21×10-6m2/s The transfer driving force can be obtained from eq. (4.380), i.e., c p (T∞ − Tw ) + hc , o (ωo ,∞ − ωo , w ) Z= hsv 1.063 × (27 − 727) + 12000 × (0.21 − 0.1) 1500 = 0.3839 = ln(1 + Z ) ln(1 + 0.3839) = = 0.1443 2.6 Z 0.15 2.6 × 0.38390.15 The blowing parameter obtained from eq. (4.386). is B= The blowing velocity at the surface is obtained from eq. (4.383): ρ ρ 1/ 2 vw = ∞ BU ∞ Re −1/ 2 = ∞ B (U ∞ν ) x −1/ 2 x ρw ρw which can be integrated to yield the average blowing velocity: 2ρ 1/ 2 vw = ∞ B (U ∞ν L ) ρw 1/ 2 2 × 1.1614 × 0.1443 × (1× 60.21× 10−6 × 1) 0.3482 = 0.007469 m/s = 4.11 Turbulence 4.11.1 Turbulent Boundary Layer Equations The generalized governing equations for three-dimensional turbulent flow have been presented in Chapter 2 (see Section 2.5). For two-dimensional steady-state turbulent flow, the governing equations can be simplified to: ∂u ∂v + =0 ∂x ∂y (4.387) 412 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u  ∂ 2 u ∂ 2 u   ∂ u′2 ∂ v′u′  ∂u ∂u 1 ∂p +v =− +ν  2 + 2  −  +    ∂x ∂y ∂y   ∂x ∂y  ρ ∂x  ∂x (4.388) (4.389) (4.390) (4.391) u  ∂ 2 v ∂ 2 v   ∂ v′2 ∂ u′v′  1 ∂p ∂v ∂v +v =− +ν  2 + 2  −  +  ρ ∂y ∂x ∂y ∂y  ∂y   ∂x  ∂x   u u  ∂ 2T ∂ 2T ∂T ∂T +v =α 2 + 2 ∂x ∂y ∂y  ∂x   ∂ u′T ′ ∂ v′T ′  +  −  ∂y    ∂x   ∂ 2ω ∂ 2ω   ∂ u′ω ′ ∂ v′ω ′  ∂ω ∂ω +v = D 2 + 2  − +  ∂x ∂y ∂y   ∂x ∂y   ∂x   To obtain the turbulent boundary layer governing equations, scale analysis can be performed to eqs. (4.388) – (4.391). While the treatments of the time-averaged quantities are similar to the cases of laminar flow, special attention must be paid to the time averaging of the products of the fluctuations. It can be shown through a scale analysis that the first terms in the last parentheses on the right hand side of eqs. (4.388) – (4.391) are negligible compared to the second terms in the parentheses (Oosthuizen and Naylor, 1999). Therefore, eqs. (4.388) – (4.391) can be simplified to: u ∂u ∂u ∂ 2 u ∂ v′u′ 1 dp +v =− +ν 2 − ρ dx ∂x ∂y ∂y ∂y u u ∂T ∂T ∂ 2T ∂ v′T ′ +v =α 2 − ∂x ∂y ∂y ∂y ∂ω ∂ω ∂ 2ω ∂ v′ω ′ +v =D 2 − ∂x ∂y ∂y ∂y (4.392) (4.393) (4.394) where eq. (4.389) became ∂p / ∂y = 0 and the partial derivative of time-averaged pressure in eq. (4.388) became: ∂p / ∂x = dp / dx , which has been reflected in eq. (4.392). When molecules or eddies in a turbulent flow cross a control surface, they will carry momentum with them. Thus, the shear stress in a turbulent flow can be caused by both molecular and eddy level activities. While the molecular level activity is the only mechanism of shear stress, transport of momentum by eddies can only be found in turbulent flow. The time-averaged shear stress tensor can be expressed as τ = τm + τt (4.395) m t where τ is the contribution of the molecular motion, and τ is caused by eddy level activity – referred to as Reynolds stress. For two-dimensional flow, the shear stress can be expressed as τ yx = μ ∂u − ρ u′v′ ∂y (4.396) Chapter 4 External Convective Heat and Mass Transfer 413 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Similarly, the heat flux in the turbulent flow can also be caused by molecular level activity (conduction) and eddy level activity. The time-averaged heat and mass flux can be respectively expressed as ′′ ′′ ′′ qy = q y m + q y t = − k ∂T + ρ c p v′T ′ ∂y (4.397) The mass flux can be obtained by a similar way: ∂ω  ′′  ′′  ′′ my = my m + my t = − ρ D 1 + ρ v′ω1′ ∂y (4.398) 4.11.2 Algebraic Models for Eddy Diffusivity In order to model turbulent flow, one must express the turbulent transports in terms of time-averaged quantities. Such relationships cannot be derived merely from the first principle and therefore must rely on empirical or semi-empirical approaches. As stated in the preceding subsection, the transport quantities for turbulent flow can be expressed as a sum of molecular and eddy effects. The contributions from the molecular level activities (laminar) are proportional to the gradients of the averaged physical quantities. If it is assumed that the contribution from the eddy level activities is also proportional to the gradients of the averaged quantities, one has − ρ u′v′ = μ t ∂u ∂y ∂T ∂y ∂ω1 ∂y (4.399) (4.400) (4.401) ρ c p v′T ′ = −k t ρ v′ω1′ = − ρ Dt where μ t , k t , and Dt are turbulent viscosity, conductivity and mass diffusivity, respectively. Substituting eqs. (4.399) – (4.401) into eqs. (4.396) – (4.398), we have τ yx = ( μ + μ t ) ′′ qy = − k ∂u ∂u = ρ (ν + ε M ) ∂y ∂y (4.402) (4.403) (4.404) ∂T ∂T  ν ε  ∂T − kt = − ρ c p  + Mt  ∂y ∂y  Pr Pr  ∂y ∂ω1  ν ε  ∂ω = − ρ  + Mt  1 ∂y  Sc Sc  ∂y  ′′ my = − ρ ( D + D t ) where εM is the eddy diffusivity for momentum, and Prt and Sct are turbulent Prandtl and Schmidt numbers, which are defined as ε Pr t = M (4.405) εH 414 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press εM (4.406) εD where ε H and ε D are eddy diffusivities for heat and mass, respectively. Substituting eqs. (4.399) – (4.401) into eqs. (4.392) – (4.394), the momentum, energy, and species equations for turbulent boundary layer become: Sct = u ∂u ∂u 1 dp ∂2u +v =− + (ν + ε M ) 2 ρ dx ∂x ∂y ∂y u u ∂T ∂T  ν ε M  ∂ 2T +v = + ∂x ∂y  Pr Pr t  ∂y 2   ∂ω ∂ω  ν ε M +v = + ∂x ∂y  Sc Sct  ∂ ω  ∂y 2  2 (4.407) (4.408) (4.409) which – together with eq. (4.387) – are governing equations for turbulent boundary layer. Appropriate models for eddy diffusivity for momentum, and turbulent Prandtl and Schmidt numbers are needed in order to describe the transport phenomena in the turbulent boundary layer. Until a viable expression of eddy diffusivity for momentum becomes available, the mathematical description of turbulent boundary layer is not complete. As for the turbulent Prandtl and Schmidt numbers, they are often assumed to be constants near unity because the mechanisms of turbulent transport of momentum, heat and mass are the same. Therefore, our attention now is turning to the models of eddy diffusivity for momentum. Mixing Length Model The mixing length model proposed by Prandtl is the simplest turbulent model. The distinctive feature of turbulent flow is the existence of eddies and vortices so that the transport in the turbulent flow is dominated by packets of the molecules, instead of the behavior of individual molecules. Considering a turbulent flow near a flat plat as shown in Fig. 4.33, the mixing length can be defined as the y u+ ∂u l ∂y u = f ( y) u′ v′ l B A u Figure 4.33 Mixing length model u Chapter 4 External Convective Heat and Mass Transfer 415 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press maximum length that a packet can travel vertically while maintaining its time averaged velocity unchanged. The concept of mixing length for turbulent flow is similar to the mean free path for random molecular motion. When a fluid packet located at point A travels to point B by moving upward a distance that equals the mixing length, l, its time-averaged velocity should be kept at u according to the definition of mixing length. On the other hand, the time-averaged velocity at point B is u + (∂u / ∂y)l according to the profile of the time-averaged velocity (see Fig. 4.33). Therefore, the packet must have a negative velocity fluctuation equal to −(∂u / ∂y)l in order to keep its time averaged velocity unchanged. Thus, the fluctuation of the velocity component in the x-direction is:  ∂u  u′ = −l    ∂y  (4.410) When the velocity component in the x-direction has the above negative fluctuation, the velocity component in the y-direction must have a positive fluctuation, v′ , with the same scale, i.e.,  ∂u  v′ = Cl    ∂y  (4.411) where C is a local constant. Thus, the time-average of the product of the velocity fluctuations, u′v′ , must be negative for this case. Similarly, we can also analyze motion of the fluid packet from point B to point A, in which case u′ will be positive and v′ will be negative. Therefore, u′v′ must be negative for any cases. Combining eqs. (4.410) and (4.411) yields  ∂u  −u′v′ = Cl 2    ∂y   ∂u  −u′v′ = l 2    ∂y  2 Since l is still undetermined, it will be beneficial to absorb C into l and yield 2 (4.412) It follows from the definition of the eddy diffusivity that [see eqs. (4.399) and (4.402)]: ε M = l2 ∂u ∂y (4.413) where the absolute value is to ensure a positive eddy diffusivity. While the general rule for determining the mixing length, l, is lacking, the mixing length for turbulent boundary layer cannot exceed the distance to the wall. Therefore, we can assume: l = κy (4.414) where κ is an empirical constant with order of 1, and is referred to as von Kármán’s constant. Equation (4.414) is valid only if κ is really a constant. Substituting eq. (4.414) into eq. (4.413), the eddy diffusivity of momentum becomes 416 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ε M = κ 2 y2 ∂u ∂y (4.415) Substituting eq. (4.415) into eq. (4.402), the shear stress in the twodimensional turbulent flow becomes τ yx = ρ ν + κ 2 y 2   ∂u  ∂u  ∂y  ∂y (4.416) Two-Layer Model The above mixing length theory assumes that the turbulent flow takes place within the entire boundary layer thickness. In reality, the turbulent boundary layer can be further divided into two regions (see Fig. 4.34). The first region is viscous sublayer and it occupies less than 1% of the total turbulent boundary layer thickness. The momentum and heat transfer in this region are dominated by viscous shear and heat conduction, respectively. Outside of the viscous sublayer is a full turbulent region which comprises almost the entirety of the turbulent boundary layer. Figure 4.34(a) shows the instantaneous velocity profiles in a turbulent boundary layer on a flat plate taken at 17 instances. The 17 superimposed profiles as well as the average profile are shown in Fig. 4.34 (b) and (c), respectively. It can be clearly seen from 4.34(c) that the velocity increases sharply in the sublayer near the wall and its change in the second region is fairly insignificant. While the thickness of the turbulent boundary layer increases with increasing x, the viscous sublayer remains at a fairly constant thickness. As the thickness of the turbulent boundary layer increases, the viscous sublayer inhabits a smaller and smaller portion of the entire turbulent boundary Figure 4.34 Velocity profiles in the turbulent boundary layer on a flat plate: (a) velocity distributions at different instances, (b) 17 superimposed profiles, and (c) time-averaged profile (Cebeci, 2004; Reprinted with permission from Elsevier). Chapter 4 External Convective Heat and Mass Transfer 417 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press layer. Therefore, the mixing length model based on the assumption that the entire turbulent boundary layer is full turbulent flow should be improved by considering the two-layer structure of the turbulent boundary layer. In the two-layer structure of the turbulent boundary layer, the velocity is dominated by viscous shear stress in the sublayer and by turbulent mixing in the fully turbulent region. Near the wall, the shear stress is equal to the shear stress at the wall, i.e., τ yx = τ w . Thus, eq. (4.402) near the wall becomes τw du = (ν + ε M ) (4.417) ρ dy where the partial derivative ∂u / ∂y has been approximated as du / dy because u (∂u / ∂x ) ≈ 0 in the region near the wall. For the viscous sublayer, ε M = 0 and eq. (4.417) is simplified to: τ du = w dy (4.418) μ Equation (4.418) can be integrated from the wall to yield u yτ w 0 du = 0 μ dy i.e., τ u= wy (4.419) μ Introducing the wall coordinate y+ = yuτ ν , u+ = u uτ (4.420) where cf τw = U∞ 2 ρ is shear velocity or friction velocity, eq. (4.419) becomes uτ = u+ = y+ (4.421) (4.422) which is the velocity distribution in the viscous sublayer. For the fully turbulent region, ε M  ν and eq. (4.417) is simplified to: τw du = εM (4.423) ρ dy which can be rewritten to the following format using the Prandtl’s mixing length model τ w = ρκ 2 y 2   du    dy  2 (4.424) Substituting eq. (4.420), eq. (4.424) can be nondimensionalized into du + 1 = dy + κ y + (4.425) In order to integrate eq. (4.425), we must know the thickness of the viscous sublayer. For the case without blowing or suction on the wall and zero-pressure 418 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press gradient, we can choose y + = 10.8 as the dimensionless thickness of the viscous sublayer. Therefore, eq. (4.425) can be integrated from y + = 10.8 , i.e.,  u+ 10.8 du + =  y+ 10.8 1 dy + κ y+ (4.426) which can be integrated to yield (4.427) where the von Kármán’s constant is taken as κ = 0.4 . Equation (4.427) is often referred to as the law of the wall. It can also be approximated into the following power law form: u + = 8.75( y + )1/ 7 (4.428) + which fits into eq. (4.427) at least up to y = 1500 . The two-layer model can then be summarized as  y+ y + < 10.8 u+ =  + + 2.44 ln y + 5.0 y > 10.8 u + = 2.44 ln y + + 5.0 (4.429) Figure 4.35 shows the comparison between the above two-layer model and experimental results obtained using water and air. The Reynolds number shown in the figure is the momentum thickness Reynolds number defined as Reδ = U ∞δ 2 / ν , where δ 2 = 0.664 ν x / U ∞ is the momentum thickness. It can be 2 seen that the results from the two-layer model agreed with the experimental results very well except at the outer region of the turbulent boundary layer. In addition, the agreement between eq. (4.429) and the experimental results is also not very good in the region near y + = 10.8 . Some researchers suggested there is a buffer region ( 5 < y + < 30 ) between the sublayer and the fully turbulent region 30 Johnson and Johnston (1989) Johnson (1989) 25 20 u+=2.44ln(y+)+5.0 u+ 15 u+=y+ 10 Reδ2=1500 5 0 1 10 100 y+ 1000 10000 Figure 4.35 Comparison of two-layer model and experimental results. Chapter 4 External Convective Heat and Mass Transfer 419 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press and the velocity profile in the buffer region is (Bejan and Kraus, 2003): u + = 5 + 5ln( y + / 5) (4.430) Van Driest Model While the above two-layer model works well in the viscous sublayer and the fully turbulent region, a discrepancy between the two-layer model and experimental results in the transition region (near y + = 10.8 ) is apparent. The reason for this discrepancy is that the velocities in these two regions were assumed to be dominated by different mechanisms: molecular activities in the viscous sublayer and eddy level activities in the fully turbulent region. The contribution of the eddy diffusion in the viscous sublayer was completely neglected. In fact, relatively large elements of fluid in the viscous sublayer can lift off the wall and immediately be replaced by the other fluid from the fully turbulent region. This phenomenon is referred to as “bursting” (Kays, et al., 2005). Therefore, the eddy diffusivity in the viscous sublayer is not really zero and the only point that it becomes zero is at the wall where y + = 0 . To account for the effect of eddy diffusivity in the viscous sublayer, the Van Driest hypothesis states that the eddy diffusivity decreases as one nears the wall (Van Driest, 1956), instead of two regions where it is “on” or “off.” Therefore, we can use the Prandtl’s mixing length model in the entire turbulent boundary layer by introducing a damping function, i.e., l = κ y(1 − e y / A ) (4.431) 25 Reδ2=1500 Johnson(1989) Johnson and Johnston (1989) 20 u+=2.44ln(y+)+5.0 15 u+ 10 u+=y+ Van Driest model 5 ∂u +)+( k y+ 2 exp(- + 25 2  ∂u+  + 2 + 1 1=(du /dy (0.41y + [1-(1 − e yy / /25)] du /dy) ) ) = + + +   ∂y  + + 2  +  ∂y    0 1 10 y+ 100 1000 Figure 4.36 Comparison of Van Driest model (A+=25) with two-layer model and experimental results. 420 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where A is an effective sublayer thickness that must be determined empirically. Equation (4.431) can also be expressed in the following dimensionless form: l = κ y(1 − e y / A ) (4.432) + where A = Auτ / ν is the dimensionless sublayer thickness. The value of A+ must be determined in such a way that the calculated value of u+ outside of the sublayer matches the result obtained by the law of the wall, eq. (4.427). It is found that A + = 25 can produce the best match. Figure 4.36 shows comparison of the results obtained by the Van Driest model, the two layer model, and experiments. It can be seen that the results produced by the Van Driest model agreed very well with the experimental results in the entire turbulent boundary layer. It should be pointed out that A + = 25 can only produce good results for the cases without blowing or suction and without pressure gradient. When adverse pressure gradient ( dp / dx > 0 ) and wall velocity are present, the following correlation can be used: + + A+ = 25 + + a{vw + b[ p + / (1 + cvw )]} + 1 vw μ dp / dx , p + = 1/ 2 3/ 2 ρ τw uτ (4.433) where + vw = (4.434) are respectively the blowing velocity and the pressure gradient in the wall coordinate. The empirical constants in eq. (4.433) are a = 7.1, b = 4.25, and c = 10.0 . If p + > 0 , b = 2.9 and c = 0 . Example 4.5: Drive the Van Driest Equation shown in Fig. 4.36. Solution: Combining eqs. (4.414) and (4.416) one obtains: τw  ∂u  ∂u = ν + l 2   ∂y  ∂y ρ  (4.435) which can be rewritten as the following form using the wall coordinate given in eq. (4.420): 2 2  ∂u +  y +   ∂u +   1 =  + + l2    +    ∂y  y   ∂y     (4.436) Using the Prandtl’s mixing length from eq. (4.432), the following expression can be written: + y+   ∂u + (κ y + ) 2 1 − e A+ 1=   ∂y +   2   ∂u +  2      ∂y +      (4.437) which is identical to the equation in Fig. 4.36 if seting A+ = 25. Solving for ∂u + / ∂y + from eq. (4.437) and considering only the physically possible solution, one obtains the final non-linear ODE of the Van Driest equation. Chapter 4 External Convective Heat and Mass Transfer 421 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press +2 y /A 2 ) ∂u + −1 + 1 + 4(κ y ) (1 − e = + + + +2 y /A 2 ∂y 2(κ y ) (1 − e ) + + (4.438) The velocity profile of Van Driest equation shown in Fig. 4.36 can be obtained by solving eq. (4.438) using a forth-order Runge–Kutta method. 4.11.3 K-ε Model In the turbulent models presented in the preceding section, the eddy diffusivity, ε M , was obtained by algebraic expressions. While they are very easy to apply, their shortcoming is their lack of universal applicability because all of them work well only near the wall. As evidenced by Figs. 4.35 and 4.36, the agreement between the wall function and the experimental results is not very good for large y+. When the above algebraic models are applied to the internal flow, they will break down near the center of the tube. Thus, it is imperative to develop a turbulent model that is applicable in the region that is sufficiently far away from the wall. An ideal turbulent model should be applicable to any region in the turbulent flow although it may not be able to be represented by algebraic equations. Significant advancement of high-speed computers makes developments and applications of non-algebraic turbulent models possible. Among different non-algebraic turbulent models, the K-ε model is one of the most widely used models and it will be discussed below. The K-ε model is based on the analogy between the motion of the fluid packet in the turbulent flow and the random motion of an ideal gas. The kinetic theory of gases stated that the kinematic viscosity of gas can be obtained by ν = cλ / 3 where c is the mean speed of the molecules and λ is the mean free path. For turbulent flow, the kinetic energy due to the fluid packet random motion is K= and the mean speed of the fluid packet is K1/2. The eddy diffusivity for the momentum can be defined in a similar way to the kinematic viscosity of gas (Kolmogorov, 1942), i.e., ε M = C μ K 1/ 2 L (4.440) where C μ is a dimensionless empirical constant and L plays the same role as the mean free path in the kinematic viscosity of the gas. Therefore, both K and L must be determined in order to use eq. (4.440) to evaluate the eddy diffusivity. For three-dimensional flow, the equation for kinetic energy associated with the velocity fluctuation was derived in Section 3.5.2. While we can obtain the equation for K by simplifying eq. (2.528) for two-dimensional turbulent boundary flow (Oosthuizen and Naylor, 1999), a more direct approach suggested by Bejan (2004) will be introduced here. If we consider a control volume in a 12 (u ′ + v ′2 + w ′2 ) 2 (4.439) 422 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press horizontal slender flow region (turbulent boundary layer), balance of the kinetic energy for velocity fluctuation is Convection of K = Eddy diffusion of K + Rate of K generation − Rate of K destruction (4.441) (4.442) For a two-dimensional steady state flow, the convection of K is Convection of K = u ∂K ∂K +v ∂x ∂y Since eddy diffusion in the x-direction is negligible for a boundary layer type flow, only diffusion in the y-direction needs to be considered: ∂  ε M ∂K  Eddy diffusion of K = (4.443)   ∂y  σ k ∂y  where σ K is a dimensionless empirical constant. The rate of K production can be obtained by the eddy shear stress and the time averaged velocity gradient  ∂u   ∂u  1/ 2  ∂u  Rate of K generation =  ε M ⋅  = Cμ K L   ∂y   ∂y    ∂y  2 (4.444) where eq. (4.440) was used to evaluate the eddy diffusivity. The rate of K destruction, ε , can be evaluated by analyzing a fluid packet with diameter L, oscillating in the turbulent flow field. If the oscillating velocity is K1/2, the drag force acting on the fluid packet will be C D ρ L2 ( K 1/ 2 )2 , where CD is the drag coefficient that is approximately equal to 1. The mechanical power dissipated per unit mass is [C ρ L2 ( K 1/ 2 ) 2 ]K 1/ 2 K 3/ 2 = CD (4.445) ε= D L ρ L3 The K-equation for a boundary layer type flow can be obtained by substituting eqs. (4.442) – (4.445) into eq. (4.441) 2  ∂u  ∂K ∂K ∂  ε M ∂K  (4.446) u +v =   + εM   −ε ∂x ∂y ∂y  σ k ∂y   ∂y  By following a similar procedure, the dissipation equation can be obtained as the following: 2  ∂u  ε ε2 ∂ε ∂ε ∂  ε M ∂ε  (4.447) u +v = − C2   + C1ε M   K ∂x ∂y ∂y  σ ε ∂y   ∂y  K where σ ε , C1 , and C 2 are additional dimensional empirical constants. After K and ε are obtained, the eddy diffusivity can be obtained by eliminating L between eqs. (4.440) and (4.445), i.e. (4.448) ε where CD in eq. (4.445) has been set to 1. Equations (4.446) – (4.448) become the equations for the K-ε model. It is different from the algebraic equations presented in the preceding subsection because additional partial differential equations must be solved. Similar to eq. (4.407), both eqs. (4.446) and (4.447) ε M = Cμ K2 Chapter 4 External Convective Heat and Mass Transfer 423 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press are parabolic because diffusion in the x-direction is neglected. Jones and Launder (1972) suggested the following values for the K-ε model: C μ = 0.09, C1 = 1.44, C 2 = 1.92, σ K = 1, and σ ε = 1.3 Since the effect of molecular viscosity, ν, was neglected in eqs. (4.446) and (4.447), the K-ε model cannot be applied in the viscous sublayer where molecular viscosity dominates the velocity profiles. In the sublayer region, Prandtl’s mixing length theory is valid and eq. (4.413) can be rewritten as 2 2  ∂u  εM (4.449)  = 4  ∂y  l The generation of K represented by the second term on the right-hand side of eq. (4.446) becomes 2  ∂u  ε3 εM   = M 4  ∂y  l Substituting eq. (4.440) into the right-hand side of the above equation, one obtains: εM  3 C μ K 3/ 2 L3  ∂u  =  l4  ∂y  2 Considering eq. (4.445), the above equation can be rewritten as C 3  L 4  ∂u  εM   = μ   ε CD  l   ∂y  2 (4.450) If the length scale is taken as C L = l D  C3 μ     1/ 4 eq. (4.450) becomes  ∂u  εM   = ε  ∂y  2 (4.451) where the left and right hand sides, respectively, represent the rates of production and dissipation of K, i.e. generation of K is balanced by dissipation of K in the viscous sublayer. Therefore, the mixing length theory can be used in the sublayer region and the fully turbulent region can be described using the K-ε model. 4.11.4 Momentum and Heat Transfer for Turbulent Flow over a Flat Plate To obtain the boundary layer thickness for turbulent flow over a flat plate, von Kármán’s momentum integral can be employed. The integral momentum equation (4.167) that was derived for the case of laminar flow is still valid except that the instantaneous velocity should be replaced by the time-averaged velocity. For the flat plate without bellowing or suction ( vw = 0 ) and without pressure gradient, eq. (4.167) is simplified to 424 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press τw d δ = (4.452) u (U ∞ − u ) dy ρ dx 0 While the velocity profile in the laminar boundary layer can be adequately described by a polynomial function, the velocity profile in a turbulent flow is too complicated to be described by a single function for the entire boundary layer. Equation (4.428) states that the velocity profile in the turbulent boundary layer can be approximated as u + = 8.75( y + )1/ 7 , which can be rewritten into the following dimensionless form: u y =  U∞  δ  1/ 7 (4.453) Although eq. (4.453) can reasonably represent velocity profile in most parts of the boundary layer, the velocity gradients at both y = 0 and δ are incorrect: (∂u / ∂y) y = 0 → ∞ and (∂u / ∂y) y =δ ≠ 0 . To get the shear stress at the wall, let us substitute eq. (4.420) and (4.421) into eq. (4.428):  δ τw / ρ = 8.75   ν τw / ρ  U∞     1/ 7 (4.454) Solving for τ w yields 2 τ w = 0.0225ρU ∞   δU∞   ν  −1/ 4 (4.455) which is referred to as the Blasius relation and which is valid for Re x < 107 . Substituting eqs. (4.453) and (4.455) into eq. (4.452), one obtains −1/ 4 1/ 7 2/ 7 d δ 2  y  y  2  δU ∞  (4.456) 0.0225U ∞  = U ∞    −    dy  dx 0 δ   ν   δ    Performing the integration and differentiation on the right-hand side of eq. (4.456) yields the following differential equation for the boundary layer thickness: −1/ 4 7 dδ 2  δU ∞  (4.457) = 0.0225U ∞  ν 72 dx   which can be integrated to obtain: ν     U∞  1/ 4 x = 3.45δ 5/ 4 + C (4.458) where C is the unspecified integration constant. If we assume that the turbulent boundary layer starts from the edge of the flat plate – which is not a good assumption (see Fig. 4.3) – the integration constant C becomes zero, and eq. (4.458) becomes δ 0.376 = (4.459) 1/5 x Re x Chapter 4 External Convective Heat and Mass Transfer 425 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δturbulent ≈ 4δlaminar 1.2 δlaminar 1 0.8 u/U∞ 0.6 Turbulent flow 0.4 Laminar flow 0.2 0 0 0.5 1 1.5 2 2.5 3 3.5 4 y/δlaminar Figure 4.37 Comparison of velocity profiles in laminar and turbulent boundary layers. The local friction coefficient can be found as τw 0.0576 cf = = (4.460) 2 ρU ∞ / 2 Re1/5 x It should be pointed out that eqs. (4.459) and (4.460) are valid for the case that turbulent boundary layer starts from the leading edge of the flat plate and Re x < 107 – beyond which the Blasius relation becomes invalid. Figure 4.37 shows a comparison of the velocity profiles of laminar and turbulent boundary layers at Re x = 5 × 106 (Welty et al., 2000). It can be seen that the turbulent boundary layer is much thicker than the laminar boundary layer at the same Reynolds number. The mean velocity in the turbulent boundary layer is much higher than that of the laminar boundary layer. The large mean velocity of the turbulent boundary layer allows for a much stronger momentum, heat and mass transfer. Another advantage of the turbulent boundary layer is that it can resist separation better than a laminar boundary layer. Due to its strong ability on momentum, heat and mass transfer and resisting separation, a turbulent boundary layer is desirable in many engineering applications. Example 4.6: Obtain Nusselt number for turbulent boundary layer based on analogy between momentum and heat transfer. The Prandtl number and turbulent Prandtl number can be assumed to be equal to 1. Solution: The shear stress and heat flux for turbulent boundary layer can be obtained from eqs. (4.402) and (4.403). These two equations can be rewritten as 426 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂u τ = ∂y ρ (ν + ε M ) ∂T q′′ =− ∂y ρ c p (ν / Pr + ε M / Pr t ) (4.461) (4.462) where the subscripts and bar for shear stress and heat flux have been dropped to simplify the notation. Equations (4.461) and (4.462) can be integrated from the wall (y = 0) to a location y =  in the free stream where u = U ∞ and T = T∞ :  τ U∞ =  dy (4.463) 0 ρ( + ε ) νM Tw − T∞ =   0 ρ c p (ν / Pr + ε M / Pr t ) q′′ dy (4.464) When Pr = Pr t = 1 , the momentum and thermal boundary layers have the same thickness and the distributions for τ and q′′ are similar, i.e.: τ q′′ (4.465) = ′′ τ w qw Therefore, eq. (4.464) becomes: q′′  τ (4.466) Tw − T∞ = w  dy 0 ρ ν +ε τ wcp ( M) Combining eqs. (4.463) and (4.466) to eliminate the integral yields: Tw − T∞ = ′′ qwU ∞ τ wcp The local heat transfer coefficient becomes: τ wcp c f ′′ qw ρ c pU ∞ = = hx = Tw − T∞ U∞ 2 The local Nusselt number is therefore h x c f  U ∞ x   μcp  c f Nu x = x = Re x Pr  = k 2  ν  k  2   Substituting eq. (4.460) into the above equation, and considering Pr =1 yields: Nu x = 0.0288 Re0.8 (4.467) x Equation (4.467) is valid only when Pr = Pr t = 1 . The shear stress for the entire turbulent boundary layer can be obtained from eq. (4.402), which can be non-dimensionalized using eqs. (4.420) and (4.421), i.e., τ  ε M  ∂u + (4.468) = 1+ τw  ν  ∂y +   Defining the dimensionless temperature in the wall coordinate Chapter 4 External Convective Heat and Mass Transfer 427 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press T+ = (Tw − T ) τ w / ρ ′′ qw / ( ρ c p ) (4.469) the heat flux in turbulent boundary layer, eq. (4.403), can be nondimensionalized to q′′  1 ε M / ν  ∂T + (4.470) = + + t ′′ qw  Pr Pr  ∂y To consider heat transfer in a turbulent boundary layer, the entire turbulent boundary layer is divided into an inner region ( y + < 5 ), a buffer region ( 5 ≤ y + ≤ 30 ), and an outer region ( y + > 30 ). In the inner region, ′′ ε M = 0 and q′′ = qw , and eq. (4.470) becomes 1= 1 ∂T + Pr ∂y + (4.471) which can be integrated to yield: T + = Pr y + (4.472) ′′ For the buffer region τ = τ w , q′′ = qw , and the velocity can be obtained by (4.430). The velocity gradient in the buffer region is ∂u + 1 =+ + ∂y y (4.473) Substituting eq. (4.473) into eq. (4.468) and considering τ = τ w , an expression for the eddy diffusivity in the buffer region is obtained: ε M y+ = −1 (4.474) ν 5 ′′ Substituting eq. (4.474) into eq. (4.470) and considering q′′ = qw , the following expression is obtained:  1 y + − 5  ∂T + 1=  + + t  Pr 5 Pr  ∂y (4.475) Integrating eq. (4.475) over the buffer region yields  i.e., T+ 5 Pr dT + =  y+ 5 dy + 1/ Pr + ( y − 5) / (5 Pr t ) + (4.476)   Pr  y + − 1  , 5 < y + < 30 (4.477) T + − 5 Pr = 5 Pr t ln 1 + t    Pr  5 At the top of the buffer region where y + = 30 , eq. (4.477) becomes  5 Pr  Tb+ = 5 Pr + 5 Pr t ln  1 + t  (4.478) Pr   For the outer region, ε M  ν and ε H  α , so eqs. (4.468) and (4.470)become τ ε M ∂u + (4.479) = τ w ν ∂y + 428 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press q′′ ε M / ν ∂T + = ′′ Pr t ∂y + qw (4.480) Since the distributions of shear stress and the heat flux in the outer layer are similar, one can assume: τ q′′ (4.481) = ′′ τ b qb ′′ where τ b and qb are shear stress and heat flux at the bottom of the outer region + ( y = 30 ). Both shear stress and heat flux across the sublayer and buffer region ′′ ′′ are constant, i.e., τ b = τ w and qb = qw , and eq. (4.481) becomes τ q′′ (4.482) = ′′ τ w qw Substituting eqs. (4.479) and (4.480) into eq. (4.482), we have ∂T + ∂u + = Pr t + ∂y + ∂y (4.483) (4.484) (4.485) Integrating eq. (4.483) across the outer layer yields + + + T∞ − Tb+ = Pr t (u∞ − ub ) where u can be obtained from eq. (4.430), i.e., + ub = 5(1 + ln 6) + b Substituting eq. (4.485) into eq. (4.484), we get: (4.486) The dimensionless temperature at the free stream can then be obtained by substituting eq. (4.478) into eq. (4.486):  Pr t + 5 Pr  + + t T∞ = 5 Pr + 5 Pr t ln   + Pr (u∞ − 5) 6 Pr t   + + T∞ − Tb+ = Pr t [u∞ − 5(1 + ln 6)] (4.487) According to eq. (4.469), the dimensionless temperature at free stream is 1/ 2 (T − T∞ ) τ w / ρ Pr Re x  c f  + = T∞ = w (4.488)  ′′ Nu x  2  qw / ( ρ c p ) Combining eqs. (4.487) and (4.488) yields the local Nusselt number as following: Nu x = Re x (c f / 2)1/ 2 Pr t  Pr t + 5 Pr  Pr t + 5+5 ln  (u∞ − 5) + Pr  6 Pr t  Pr (4.489) The dimensionless free stream velocity can be obtained from eq. (4.420) and (4.421): + u∞ = U∞ U∞ 2 = = uτ cf τw / ρ (4.490) Substituting eq. (4.490) into eq. (4.489), the following expression for the local Nusselt number is obtained: Chapter 4 External Convective Heat and Mass Transfer 429 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Nu x = Re x Pr cf / 2 cf    Pr t + 5 Pr  t t t 5 Pr + 5 Pr ln   − 5 Pr  + Pr t 2  6 Pr   (4.491) Assuming Prt = 1 and evaluating the friction coefficient using eq. (4.460), the Nusselt number for turbulent boundary layer becomes Nu x = 0.288 Re0.8 G (4.492) x where G= Pr 0.0288   1 + 5 Pr   5 Pr + 5ln   − 5 + 1 Re0.2  6  x (4.493) References Bejan, A., 2004, Convection Heat Transfer, 3rd ed., John Wiley & Sons, Hoboken, NJ. Bejan, A., and Kraus, A.D., 2003, Heat Transfer Handbook, John Wiley & Sons, Hoboken, NJ. Blasius, H., 1908, Grengschichten in Flüssigkeiten mit kleiner Reibung, Z. Math. Phy., Vol. 56 p.4, Also NACA TM 1256 Brackbill, J.U., Kothe, D.B., and Zemach, C., 1992, “A Continuum Method for Modeling Surface Tension,” Journal of Computational Physics, Vol. 100, pp. 335-354. Burmeister, L.C., 1993, Convective Heat Transfer, 2nd ed., John Wiley & Sons, Hoboken, NJ. Cebeci, T., 2004, Analysis of Turbulent Flows, 2nd Revised and Expanded Edition, Elsevier, Oxford, UK. Eckert, E.R.G., and Goldstein, R.J., 1976, Measurement in Heat Transfer, McGraw-Hill, New York, NY Eckert, E.R.G., 1942, Die Berchnung des Wärmeüberganges in der Laminaren Grenzschicht um strömter körper, VDI – Forschungsheft, vol. 416 pp. 1-24 Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems, Elsevier, Burlington, MA. Falkner, V.M., and Skan, S.W., 1931, Some Appropriate Solutions of the Boundary Layer Equations, Phil. May 12, pp. 865-896 Glassman, I., 1987, Combustion, 2nd ed. Academic Press, Orlando, FL Harlow, F.H., and Welch, J.E., 1965, “Numerical Calculation of Time-Dependent Viscous Incompressible Flow,” Physics of Fluids, Vol. 8, pp. 2182-2189. 430 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Hirt, C.W. and Nichols, B.D., 1981, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” Journal of Computational Physics, Vol. 39, pp. 201-225. Incropera, F.P., Dewitt, D.P., Bergman, T.L., and Lavine, A.S., 2006, Fundamentals of Heat and Mass Transfer, John Wiley & Sons, Hoboken, NJ. Jamet, D., Torres, D. and Brackbill, J.U., 2002, “On the Theory and Computation of Surface Tension: The Elimination of Parasitic Currents through Energy Conservation in the Second-Gradient Method,” Journal of Computational Physics, Vol. 182, pp. 262-276. Jones, W.P., and Launder, B.E., 1972, “The Prediction of Laminarization with a Two-Equation Model of Turbulence,” Int. J. Heat Mass Transfer, Vol. 15, pp. 301-314. Johnson, P.L., and Johnston, J.P., 1989, “The Effects of Grid-Generated Turbulence on Flat and Concave Turbulent Boundary Layers,” Seventh Symposium on Turbulent Shear Flows, Stanford University, pp. 20.2.1-20.2.6. Kaviany, M., 2001, Principles of Convective Heat Transfer, 2nd ed., SpringerVerlag, New York. Kays, W.M., Crawford, M.E., and Weigand, B., 2005, Convective Heat Transfer, 4th ed., McGraw-Hill, New York, NY Kolmogorov, A.N., 1942, “Equations of Turbulent Motion of an Incompressible Turbulent Fluid,” Izv. Akad. Nauk SSSR, Ser. Fiz., Vol. VI, pp. 56-58. Kurosaki, Y., 1974, “Coupled Heat and Mass Transfer of a Flat Plate with Uniform Heat Flux in a Laminar Parallel Flow,” Journal of the Japan Society of Mechanical Engineers, part B, Vol. 40, pp. 1066-1072 Leonard, B.P., 1979, “A Stable and Accurate Convective Modeling Procedure based on Quadratic Upstream Interpolation,” Computer Methods in Applied Mechanics and Engineering, Vol. 29, pp. 59-98. Leonard 1981, B.P., “A Survey of Finite Differences with Upwinding for Numerical Modeling of the Incompressible Convection Diffusion Equation,” Computational Techniques in Transient and Turbulent Flows, Taylor, C., and Morgan, K., Eds., Pineridge Press, Swansea, pp. 1-35. Minkowycz, M.J., Sparrow, E.M., and Murthy, J.Y., eds., 2006, Handbook of Numerical Heat Transfer, 2nd ed. John Wiley & Sons, Hoboken, NJ. Nichols, B.D., and Hirt, C.W., 1975, “Methods for Calculating Multidimensional, Transient Free Surface Flows Past Bodies,” Proc. of the First International Conf. On Num. Ship Hydrodynamics, Gaithersburg, MD, Oct. 2023. Oosthuizen, P.H., and Naylor, D., 1999, Introduction to Convective Heat Transfer Analysis, WCB/McGraw-Hill, New York. Chapter 4 External Convective Heat and Mass Transfer 431 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Patankar, S.V., and Spalding, D.B., 1972, “A Calculation Procedure for Heat, and Mass and Momentum Transfer in Three-Dimensional Parabolic Flows,” International Journal of Heat and Mass Transfer, Vol. 17, pp. 1787-1806. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC. Patankar, S.V., 1981, “A Calculation Procedure for Two-Dimensional Elliptic Situations,” Numerical Heat Transfer, Vol. 4, pp. 409-425. Patankar, S.V., 1991, Computation of Conduction and Duct Flow Heat Transfer, Innovative Research. Prandtl, L., 1904, Flussigkeitsbewegung bei sehr Reibung, Proceedings of the 3rd International Mathematics Congress Heidelberg, Germany also NACA TM 452, 1928 Rice, J., Faghri, A., and Cetegen, B.M., 2005, “Analysis of a Free Surface Film from a Controlled Liquid Impinging Jet over a Rotating Disk including Conjugate Effects, with and without Evaporation,” Int. Journal of Heat and Mass Transfer, Vol. 48, pp.5192-5204. Schlichting, H. and Gersteu, K., 2000, Boundary layer theory, 8th enlarged and revised ed., Springer-Verlag New York. Spalding, D.B., 1972, “A Novel Finite-difference Formulation for Differential Expressions Involving both First and Second Derivatives,” Int. J. Num. Methods Eng., Vol. 4, pp. 551-559. Tao 2001, W.Q., Numerical Heat Transfer, 2nd Ed., Xi’an Jiaotong University Press, Xi’an, China (in Chinese). Van Driest, E.R., 1956, “On Turbulent Flow near a Wall,” J. Aero. Sci., Vol. 23, pp. 1007-1011. Welty, J.R., Wicks, C.E., Wilson, R.E., Rorrer, G., 2000, Fundamentals of Momentum, Heat, and Mass Transfer, 4th ed., Wiley, New York. White, F.M., 2005, Viscous Fluid Flow, 3rd Ed., McGraw-Hill, New York, NY. Youngs, D.L., 1982, “Time-Dependent Multimaterial Flow with large Fluid Distortion,” numerical Methods for Fluid Dynamics, K.W. Morton and M.J. Baines, eds. Academic Press, pp. 273-285. Zhang, Y., Chen, Z.Q., and Chen, M., 1996, “Local Non-Similarity Solutions of Coupled Heat and Mass Transfer of a Flat Plate with Uniform Heat Flux in a Laminar Parallel Flow,” Journal of Thermal Science, Vol. 5, pp. 112-116. 432 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Problems 4.1. Use boundary layer approximation as well as scaling and order of magnitude analysis to show that the laminar thermal boundary equation, including viscous dissipation effect, in a Cartesian coordinate system, for a Newtonian fluid, is ρcp   ∂T ∂T ∂T  ∂  ∂T  +u +v = k + ∂x ∂y  ∂y  ∂y   ∂t μ  ∂u  ∂P   ∂P +u  + βT   ∂y  ∂t ∂x    2 where β is the coefficient of thermal expansion. 4.2. Using order of magnitude analysis and scaling, show that the general, two dimensional, unsteady, Laminar boundary layer equations for a nonnewtonian fluid with variable properties in a Cartesian coordinate system are  ∂u ∂u ∂u  ∂p ∂τ momentum: ρ  + u + v  = X − + yx energy: dy  dx ∂y ∂y  ∂t  ∂T ∂T ∂T  ∂  ∂T  ∂u ∂p   ∂p ρ cp  +u +v + βT  +u  = k  + τ yx ∂t ∂x ∂y  ∂y  ∂y  ∂y ∂t ∂x    4.3. Show that eq. (4.68) is the approximate physical boundary conditions for the Blasius equation (4.66) for the case of no blowing or suction at the wall. 4.4. Show that the Blasius equation can also be developed by introducing a stream function, ψ , where u= ∂ψ ∂ψ and v = − ∂x ∂y y where ψ = ν xU ∞ ς (η ) and η = ν x / U∞ 4.5. Use the shooting method to solve eq. (4.66) and the boundary condition eq. (4.68). In this method, a value of ς ′′ (0) is guessed and the function ς is predicted up to a larger η, such as 15. If ς ′′ (15) is different than 1, an adjustment is made to find a new ς ′′ (0). This process is continued until the requirement at ς ′′ (15) is approximately satisfied. 4.6. Assume steady parallel flow of air over both sides of a flat plate, with length and width of 1 meter, at atmospheric pressure, with free stream temperature and velocity of 80 °F and 8 ft/s, respectively. Assuming the wall temperature is kept at 120 °F, a) determine the total drag force on both Chapter 4 External Convective Heat and Mass Transfer 433 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press sides of the plate, and b) determine the total heat transfer rate from the plate to the air. 4.7. Show that similarity solutions for flow over a wedge exist only if the vertical velocity at the wall varies according to vw ∝ x ( m −1)/ 2 , using the definition of the stream function in Problem 4.4. 4.8. A gas with a velocity, U∞, a concentration of a sublimable substance, w∞, and a temperature, T∞, flows parallel to a flat plate coated with sublimable materials; the back of the flat plate is adiabatic (see Fig. P4.1). Specify the governing conservation equations and the corresponding boundary conditions of the sublimation problem U∞ Figure P4.1 4.9. Suppose the blowing velocity on the surface of the flat plate in Problem 4.8 satisfies vw ∝ x −1/ 2 . Introduce appropriate similarity variables to the governing equations and reduce them into a set of ordinary differential equations. 4.10. Write a computer program to solve for the ordinary differential equations of Problem 4.9. Obtain the local Nusselt number and the Sherwood number. 4.11. Air with a temperature of 27 °C flows at 1 m/s over a 1-m-long solid fuel surface at a temperature of 527 °C. The average blowing velocity due to sublimation of the solid fuel is 0.01 m/s, and the heat released per unit mass of the oxidant consumed is 10,000 kJ/kg. The latent heat of sublimation for the solid fuel is 1350 kJ/kg. The sensible heat required to raise the surface temperature of the solid fuel to sublimation temperature, and heat loss of the solid fuel, can be neglected. Estimate the mass fraction of the oxidant at the solid fuel surface. 4.12. A water film wets a horizontal surface over which flows a free stream gas of ambient temperature T∞, mass fraction of water, ω∞, and velocity, U∞ (Figure P4.2). a) Develop the governing conservation equations with appropriate boundary conditions. b) Obtain the similarity solution for mean Nusselt number Nu , mean Sherwood number Sh , and mass fraction at the liquid surface wδ. 434 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Figure P4.2 4.13. Air at 1 atm and 30 °C with 55% humidity flows over a 0.8 m-long wet flat plate at a velocity of 1 m/sec. Estimate the rate of heat transfer to and evaporation from the wet flat plate. Estimate the evaporation-induced blowing velocity on the surface, and compare it with the incoming velocity of the humid air. 4.14. Obtain the Nussult number, Nu, and the skin friction coefficient, cf, for the case of flow over a flat plate, including the effect of viscous dissipation, using the similarity solution. 4.15. Determine whether any of the following transport phenomena analogies for flow over a wedge are valid: a) b) c) Nu = ( C f / 2 ) Re Pr1/3 Nu = ( C f / 2 ) Re Pr1/3 Nu Pr1/3 = Sh Sc1/3 Nu Pr1/3 = 1/3 d) Sh Sc 4.16. Use the integral method to find the surface temperature variation along the wall and local Nusselt number for steady, laminar flow, with constant properties over a flat plate with a constant wall heat flux at the wall. 4.17. Determine δ, δT, cf, and Nu using the integral method for steady laminar boundary layer flow over a flat plate by using the following profiles for velocity and temperature within the boundary layer: u = c1 + c2y and T – Tw = c3 + c4y where c1, c2, c3, and c4 are constants. 4.18. Repeat Problem 4.17 assuming a velocity and temperature profile of the forms Chapter 4 External Convective Heat and Mass Transfer 435 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u = c1sin(c2y) and T – Tw= c3sin(c4y) 4.19. A subcooled solid is exposed to its superheated vapor as shown in Fig. 4.12(a). The temperature at the left surface of the solid is T0 , which is below the interfacial temperature. Depending on the direction of the overall heat flux at the interface, both sublimation and deposition are possible. Derive the criteria for sublimation and deposition. 4.20. Superheated vapor is brought into contact with a cold surface at a temperature of T0, and deposition takes place on the cold surface. Find the deposition rate by solving transient conduction in the deposited solid phase. 4.21. Manned spacecraft and spacesuits reject excess thermal energy by sublimating water into the vacuum of space. The sublimator consists of a porous plate exposed to vacuum on one side and feed water on the other side. The feed water seeps into the porous plate, where it then freezes. In Fig. P4.3, the sublimator also has a separate coolant heat exchanger that interfaces with the feed water. Describe in detail how the process of sublimation keeps the astronaut or spacecraft cool. Coolant Vacuum Coolant Heat Exchange Forced Water Porous Plate Figure P4.3 4.22. Obtain the value of φP using the central difference scheme, eq. (4.218), for φW = 200 and φE = 100 with different Peclet numbers: -50, -10, -5, -1, 0, 1, 2, 5, 10, and 50. Compare your results with those obtained from the exact solution, eq. (4.210) and comment on the stability of the central difference scheme. 436 Advanced Heat and Mass Transfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 4.23. Repeat the above problem for the hybrid scheme and discuss under which grid Peclet number the largest error occurs. 4.24. If the function A(PeΔ) in eq. (4.244) for the case that PeΔ > 0 is known, show that the expression for B for any PeΔ expressed as eq. (4.251) is correct. 4.25. Integrate the continuity equation for transient two-dimensional compressible flow with respect to t in the interval of (t, t+Δt) and over the control volume, P shown in Fig. 4.21. Show that eq. (4.271) can be reduced to eq. (4.273). 4.26. The control volume shown in Fig. 4.23(a) is square in shape ( Δx = Δy ), and the pressures at E and N are 20 and 30, respectively. The velocity at the west and south faces of the control volume are uw = 16 and vs = 25 , respectively. If the velocities at the east and north faces of the control volume satisfy ue = 0.6( pP − pE ) and vn = 0.8( pP − pN ) , respectively, obtain pP, ue, and vn using the SIMPLE algorithm. 4.27. Obtain the integral solution for heat transfer in a turbulent boundary layer with an unheated starting length. The Prandtl number and turbulent number can be assumed to be equal to one for simplicity. U∞ T∞ y x x0 Figure P4.4 δT δ Tw Chapter 4 External Convective Heat and Mass Transfer 437 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press