Back to Advanced Heat and Mass Transfer Home
Table of Contents
2 GE ERA I ED G VER I G EQUATI S 2 1 I tr ducti The heat and mass transfer systems can either possess one phase or multiple phases: the former is referred to as single-phase system while the latter is referred to as multiphase system. A multiphase system, which is distinguished from a single-phase system by the presence of one or more interfaces separating the phases, can be considered as a field that is divided into single-phase regions by those interfaces, or moving boundaries, between phases. A single-phase system can be described using the standard local instance governing equations with appropriate boundary conditions. Two common techniques are used for describing fluid flow: Lagrangian and Eulerian. The Lagrangian approach requires that the properties of a particular element of the fluid particles be tracked as it traverses the flow. This approach is similar to what we used in particle and rigid-body dynamics. The location of this fluid element is described by its coordinates (x, y, z), which are functions of time. The fluid element can be identified by tracking it from its initial location (x0, y0, z0) at time t = 0, and the velocity of this element at an arbitrary time t is expressed as V = V ( x0 , y0 , z0 , t ) . In order to describe a fluid flow using the Lagrangian approach, the sensors that monitor fluid properties would have to move at the same velocity as the fluid element; this is an impractical and often impossible requirement to meet, especially for such complex cases as three-dimensional transient flow. Therefore, the Lagrangian approach is rarely used in description of fluid flow. The Eulerian approach, on the other hand, observes the flow properties from a fixed location relative to a reference frame, which can be stationary or more generally moves at its own velocity. The Eulerian approach gives the values of the fluid variable at a given point (x, y, z) at a given time t. For example, the velocity can be expressed as V = V ( x , y, z , t ) , where x, y, and z are independent of t. The Eulerian approach requires that the fluid properties be measured at spatial locations that are fixed relative to the reference frame in the fluid field, so Chapter 2 Ge era i ed G ver i g Equati s 89 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the sensors are not required to move with the individual particles. Since the Eulerian approach is consistent with conventional experimental observation techniques, it is widely used in fluid mechanics and will also be adopted in this textbook. The conservation laws for flow and heat transfer can be expressed for a fixed-mass or for a control volume. While the former is a fixed collection of particles with constant mass (fixed-mass), the latter is a defined region relative to the reference frame in space (fixed-volume). From a thermodynamic point of view, the fixed-mass and the control volume can be considered as closed and open systems, respectively. Writing the governing equations for a fixed-mass requires tracking the motion of the particles, i.e., the Lagrangian approach. As we have already mentioned, although the Lagrangian description is applicable to some fluid mechanics problems, it is not a very practical way to describe multiphase systems, with a few exceptions such as liquid droplet tracking in a thermal plasma. The governing equations written for a control volume, on the other hand, express the relationship between the change of properties inside the control volume and the property of the flow into or out of the control volume. The governing equations obtained by writing the conservation laws for a control volume are consistent with the Eulerian approach. However, all of the fundamental laws of mechanics, including conservation of mass, momentum, and energy, are formulated for a collection of particles with fixed identity; that is to say, they are Lagrangian in nature. Therefore, it is necessary to apply the fundamental laws to the fixed-mass first, and then to convert them into expressions for the control volume. Section 2.2 presents local instance macroscopic (integral) formulations of the governing equations of the heat and mass transfer in single-phase systems. The macroscopic (integral) formulation is obtained by performing mass, momentum, energy, entropy and species mass balances over a control volume that includes different phases as well as interfaces. The local instance microscopic (differential) formulation is presented in Section 2.3. The local instance microscopic (differential) formulations are obtained by simplifying the integral formulations for control volumes with only one phase. The local instance differential equations must be supplemented by the jump conditions at the interface. The classification of PDEs and boundary conditions, as well as a rarefied vapor self-diffusion model, are also discussed. Section 2.3 is closed by discussion of a rarefied vapor self-diffusion model and application of the differential formulations to combustion. Section 2.4 presents a thorough overview of various averaging methods used to describe multiphase flow and heat and mass transfer problems; this is followed by the governing equations for the multifluid and homogeneous models. Section 2.4 is closed by single- and multiphase transport phenomena in porous media. 90 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 22 acr sc pic (I tegra ) ca I sta ce F r u ati A fixed-mass system describes an amount of matter that can move, flow and interact with the surroundings, but the control volume approach depicts a region or volume of interest in a flow field, which is not unique and depends on the user. So conservation laws for a fixed-mass system need to be transformed to apply to a control volume. A relation that allows one to mathematically and physically link the conservation laws for a control volume with that of a fixedmass system will be derived. Figure 2.1 shows the flow field under consideration. At time t, the control volume shown by the solid line coincides with a single-phase fixed-mass system depicted by the dashed line. At time t+dt, a portion of the fixed-mass system moves outside of the boundaries of the control volume. It can be seen from Fig. 2.1 that region I is occupied by the system at time t only, region II is common to the system at both t and t+dt, and region III is occupied by the system at t+dt only. For a system with a fixed-mass, the change of the general property Φ, which has a specific value (per unit mass) φ, can be written as (Welty, 1978): dΦ dt = system Φ t + dt −Φ t dt (2.1) Considering that the control mass occupies regions I and II at time t, one obtains Φ t = Φ I + Φ II . At time t + dt , the fixed-mass system occupies both regions II and III, i.e., rewritten as dΦ dt = system Φ t + dt = Φ II + Φ III . Therefore, eq. (2.1) can be t + dt −Φ II t Φ II dt + Φ III t + dt dt − ΦI dt t (2.2) Figure 2 1 Relation between a fixed mass system and a control volume. Chapter 2 Ge era i ed G ver i g Equati s 91 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The first term on the right-hand side of eq. (2.2) may be written as Φ II t + dt −Φ II t dt = dΦ II dt and represents the rate of change of property Φ within the control volume, because region II becomes coincident with the control volume as dt → 0 . For the general case of variable Φ within a control volume, it is appropriate to write the time derivative of Φ as d Φ |II  d Φ  ∂ = = dt  dt CV ∂t  V ρφ dV where V is the volume of the control volume. The second and third terms on the right-hand side of eq. (2.2), respectively, represent the property Φ leaving and entering the control volume due to mass flow across its boundary. If the absolute velocity is V, and the reference frame moves with a constant velocity Vref, the relative velocity is Vrel = V − Vref . For the control volume’s entire surface area, A, the rate of movement of property Φ due to mass flow may be written as (see Fig. 2.2) Φ III t +Δt dt − ΦI dt t =  A ρ (Vrel ⋅ n)φ dA where n is the normal direction of the control volume. Therefore, eq. (2.2) can be rewritten as dΦ dt = system ∂ ∂t  V ρφ dV +  A ρ (Vrel ⋅ n)φ dA (2.3) Vrel V V = Vref + Vrel z n Vrel 0 x Vrel y Figure 2 2 Control volume in a flow field. 92 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press which is the final form of the transformation formula that relates the change of property for a fixed-mass system to that of the control volume. It states that the rate of change of a property Φ for a fixed-mass system is equal to the rate of change of Φ in the control volume (the first term on the right-hand side) plus the net rate of efflux of Φ by mass flow into or out of the control volume (the second term on the right-hand side). It should be pointed out that the control volume moves with the reference frame, which moves with a constant velocity, Vref. The coordinate system is attached to and moves with the reference frame. In other words, the coordinate system is stationary relative to the reference frame. The reference frame is inertial, so Newton’s second law is valid in the coordinate system that moves with the reference frame. Equation (2.3) will be used to obtain the macroscopic (integral) formulation of the basic laws for a control volume. 221C servati f ass The law of the conservation of mass dictates that mass may be neither created nor destroyed. For a control volume that contains only one phase, conservation of mass can be obtained by setting the general and specific property forms to Φ = m and φ = 1 in eq. (2.3), i.e., dm ∂ = dt system ∂t  V ρ dV +  A ρ (Vrel ⋅ n)dA (2.4) where the first term on the right hand side represents the time rate of change of the mass of the contents of the control volume, and the second term on the right hand side represents the net rate of the mass flow through the control surface. The term (Vrel ⋅ n)dA in the mass flow integral represents the product of the velocity component perpendicular to the control surface and differential area. This term is the volume flowrate through dA, and becomes the mass flowrate when multiplied by density, ρ. Since the mass of a fixed-mass system is constant by definition, and the fixed-mass system contains only one phase, the resulting formulation of conservation of mass is ∂ ∂t  V ρ dV +  A ρ (Vrel ⋅ n)dA = 0 (2.5) which shows that the time rate of change of the mass of the contents of the control volume plus the net rate of mass flow through the control surface must be equal to zero. In other words, the sum of the mass flow rate into and out of the control volume must be equal to the accumulation and depletion of the mass within the control volume. For a control volume containing multiple phases separated by interfaces, the conservation of mass can be similarly obtained (Faghri and Zhang, 2006): Chapter 2 Ge era i ed G ver i g Equati s 93 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press   ∂t    k =1 Π ∂ Vk (t ) ρk dV +  Ak (t ) ρk (Vk ,rel ⋅ n k )dA  = 0   (2.6) where k denotes the kth phase in the multiphase system, and Π is the number of phases. 222 e tu Equati Newton’s second law states that, in an inertial reference frame, the time rate of momentum change of a fixed mass system is equal to the net force acting on the system, and it takes place in the direction of the net force. Mathematically, Newton’s second law of motion for fixed-mass in a reference frame that moves at a constant velocity Vref (see Fig. 2.2) is written as F = d (mV )rel dt (2.7) where the left-hand side is the net force vector acting on the fixed-mass system, and the right-hand side is the rate of momentum change. For control volumes that contain only one phase, the integral form of Newton’s second law can be obtained by using eq. (2.3). With the applicable value of Φ and φ in eq. (2.3) defined as Φ = mVrel and φ = Vrel , one obtains  F = ∂t  ∂ V ρ Vrel dV +  A ρ (Vrel ⋅ n)Vrel dA (2.8) which is in the vector form and is valid in all three directions. Forces acting on the control volume include body forces and contact forces that act on its surface. For example, for a multicomponent system that contains N components, if the body force per unit volume acting on the ith species is Xi, the total body force acting on the control volume is  the mass concentration (kg/m3) of the ith species. If the body force per unit mass is the same for different species (as is the case with gravity), the body force term ρ XdV , where ρ is the density of the mixture. The stress is reduced to N   ρi Xi dV , where ρi is V   i =1   V tensor acts on the surface of a fluid control volume, and includes both normal and shear stresses. The net force may be written as  F=  N   ρi Xi dV + V   i =1   A ′ τ rel ⋅ ndA (2.9) where τ ′ is the total stress tensor and n is the local normal unit vector on surface A. The dot product of a tensor of rank two, τ ′ , and a vector, n , is a vector that represents the force acting on the surface of the control volume per unit area. Combining eqs. (2.8) and (2.9), we obtain the momentum equation for the control volume of a single-phase system: 94 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  ∂ = ∂t N   ρi Xi dV + V   i =1    A ′ τ rel ⋅ ndA (2.10)  V ρ Vrel dV + A ρ (Vrel ⋅ n)Vrel dA where the two terms on the left-hand side represent, respectively, the body force and stress on the control volume, and the two terms on the right-hand side represent, respectively, the rate of momentum change in the control volume and the rate of the momentum flow into or out of the control volume. Vrel is the bulk velocity of the mixture that contains N components. When the control volume includes multiple phases and multicomponent, integrations must be performed for each subvolume. In that case, the momentum equation becomes (Faghri and Zhang, 2006)    Π N   ρk ,i Xk ,i dV + V (t )   k =1  k  i =1  Π     τ ′ ,rel ⋅ nk dA  k Ak (t )   (2.11) ∂ =  ∂t k =1    ρk Vk ,rel dV + ρk (Vk ,rel ⋅ n k )Vk ,rel dA  Vk (t ) Ak (t )   Equations (2.10) and (2.11) are momentum equations in a coordinate system that is attached to and moves with an inertial reference frame. For a fixed coordinate system that does not move with the reference frame (while the control volume still moves with the reference frame at velocity Vref), one can substitute the general variables Φ = mV and φ = V into eq. (2.3) to obtain the momentum equation. 2 2 3 E ergy Equati The first law of thermodynamics for a fixed-mass system (closed system) states that ˆ dE δ Q δW = − (2.12) dt system dt dt where Q is positive if heat is transferred into the system, and W is positive if the work is done by the system to the surrounding. The mass of system can store energy internally in a number of different forms. Therefore, the total energy Ê is due to internal energy, E, kinetic, potential, electromagnetic, surface tension and other forms. For the development of equations in this chapter, the contributions of internal and kinetic energies are considered. Other contributions can also be added as a source term or boundary condition based on the physical model. For a control volume including only one phase, the energy equation for the control volume can be obtained by setting the general property as 2 2 Φ = E + mVrel / 2 + mgz and φ = e + Vrel / 2 + gz in eq. (2.3), i.e., Chapter 2 Ge era i ed G ver i g Equati s 95 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ˆ dE dt = system ∂ ∂t  V ρ e +    2  Vrel + gz  dV +  2   A ρ (Vrel ⋅ n)  e +    2  Vrel + gz  dA  2  (2.13) ˆ where E and e are, respectively, total internal energy and specific internal 2 energy, and Vrel = Vrel ⋅ Vrel . Substituting eq. (2.13) into eq. (2.12), an integral form of the energy equation may be written as     V2 V2 δ Q δW ∂ − = ρ  e + rel + gz  dV + ρ (Vrel ⋅ n)  e + rel + gz  dA (2.14)     dt dt ∂t  V  2   A  2  The heat flow, δ Q / dt , is attributed to conduction (and radiation, for the case of participating medium) across the boundary and/or to internal generation, i.e., δQ dt =  A −q′′ ⋅ ndA +  V q′′′dV where q′′ is the heat flux vector at the control volume surface, which can be caused by temperature or concentration gradients, as indicated in eq. (1.70). The dot product of the heat flux vector, q′′ , with the unit normal vector, n , gives heat conducted out of the control volume. The term q′′′ is the internal heat generation per unit volume, and it can be caused by chemical reaction, electrical heating, etc. There are four types of work included in the work rate term. The first is shaft work, δ Wsh / dt , which is done by the control volume to its surroundings that could cause a shaft to rotate or raise a weight through a distance. The second type is the work due to body force, δ Wb / dt . The third term is associated with pressue change, δ Wp / dt , and the remaining part of work is necessary to overcome viscous force due to normal and shear stress, δ Ws / dt . δ W δ Wsh δ Wb δ Wp δ Ws dt = + + + dt dt dt dt (2.15) The work done by the body force is N  δ Wb = −  ρi Xi  ⋅ Vrel dV dt    V i =1   (2.16) The work associated with pressure can be written as δ Wp dt =−  A p(Vrel ⋅ n) ⋅ dA (2.17) The viscous work is evaluated by taking the scalar product of the force acting on the surface of the control volume per unit area, nk ⋅ τ k , and the velocity, Vrel , over the entire surface area of the control volume: δ Ws = − (n ⋅ τ ) ⋅ Vrel dA (2.18) dt  A Substituting eqs. (2.15) and (2.17) into eq. (2.14), one obtains 96 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press δ Q δ Wsh dt = − dt   V − δ Wb dt 2 Vrel − δ Ws dt ∂ ∂t  ρ e +   + gz  dV +  2   A ρ (Vrel ⋅ n)  e +    2  Vrel p + + gz  dA  2 ρ  (2.19) For a single phase system neglecting the shaft work and including pressure work as part of τ ′ , the integral form of the energy equation (2.14) becomes ∂ ∂t  V ρ e +    2 Vrel 2   dV +    A ρ (Vrel ⋅ n)  e +    2 Vrel 2   dA   =−  A q′′ ⋅ ndA +  V q′′′dV +  A ′ (n ⋅ τ rel ) ⋅ Vrel dA +  N   ρi Xi  ⋅ Vrel dV V   i =1  (2.20) where the left-hand side represents the rate of total energy change in the control volume and the rate of total energy flow into or out of the control volume. The four terms in the right-hand side represent heat transfer across the boundary of the control volume, internal heat generation, work done by the stress (including normal, shear stress, and pressure at the boundary of control volume), and the work done by the body force. For control volumes including multiple phases, the energy equation becomes (Faghri and Zhang, 2006)   ∂t   k =1  Π Π ∂ Vk (t ) ρk  ek +    2 Vk ,rel   dV + 2   Ak (t ) ρk (Vk ,rel ⋅ n k )  ek +    2 Vk ,rel    dA  2   =   −  k =1 Ak (t ) q′′ ⋅ n k dA + k  Vk (t ) ′′′ qk dV +  Ak (t ) (nk ⋅ τ ′ ,rel ) ⋅ Vk ,rel dA k (2.21) +   N   ρk ,i Xk ,i  ⋅ Vk ,rel dV  Vk (t )     i =1    where the second integral in the bracket on the left-hand side is advection of energy due to mass flow. Example 2.1 An incompressible fluid flows in a circular tube of cross-sectional area A1, which dumps into a large tube of cross-sectional area A2 as shown in Fig. 2.3. Find the change in internal energy between (1) and (2) for steady condition. Neglect shear stress at the wall. Solution: The control volume is chosen as section 2 with the length away from the inlet. For a steady state flow problem, the continuity equation (2.5) becomes ρ (Vrel ⋅ n)dA = 0 (2.22)  A i.e., − ρ u1 A1 + ρ u2 A2 = 0 (2.23) Chapter 2 Ge era i ed G ver i g Equati s 97 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Control volume p2 A2 ρ u2 p1 A1 ρ u1 Figure 2 3 Flow through a sudden enlargement Therefore, the velocity in section 2 is u2 = A1 u1 A2 (2.24) Neglecting the body force and considering steady state condition, the momentum equation (2.10) becomes ′ τ rel ⋅ ndA = ρ (Vrel ⋅ n)Vrel dA (2.25)  A  A For control volume shown in Fig. 2.3, we have, 2 2 p1 A2 − p2 A2 = ρ u2 A2 − ρ u1 A1 which can be rearranged to p1 − p2 2 2 = u2 − u1 (2.26) (2.27) ρ Since there is no heat transfer ( δ Q / dt = 0 ) and no shaft, body force, and viscous work, δ Wsh / dt = δ Wb / dt = δ Ws / dt = 0 , and problem is steady state, ∂ / ∂t = 0 , the energy equation (2.19) becomes A1 A2  i.e., A ρ (Vrel ⋅ n)  e +      2  Vrel p + + gz  dA = 0  2 ρ  (2.28) ρ u2 A2  e2 +  2  u2 p2  u2 p  +  = ρ u1 A1  e1 + 1 + 1   2 ρ 2 ρ    (2.29) Substituting eq. (2.23) into eq. (2.29), one obtains ρ Considering eq. (2.27), eq. (2.30) becomes e2 − e1 = 2 u1  A1  1 −  2 A2  2 e2 − e1 = p1 − p2 + 2 2 u1 − u2 2 (2.30) (2.31) 98 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2 2 4 The Sec d aw f Ther dy a ics The second law of thermodynamics requires that the entropy generation in a closed system (fixed-mass) must be greater than or equal to zero. The entropy change for a system with fixed-mass and only one phase can be obtained by setting Φ = S , φ = s in eq. (2.3), i.e., dS dt = system d dt  V ρ sdV +  A ρ (Vrel ⋅ n)sdA (2.32) The change of entropy in a closed system results from heat transfer and/or entropy generation: dS dt = system  A −q′′ ⋅ n dA + T  V q′′′ dV + T  V  ′′′ sgen dV (2.33) where, on the right-hand side of eq. (2.33), the first term represents the change of entropy due to heat transfer across the boundary of the control volume, and the second term represents the change of entropy due to internal heat generation in the control volume. The last term represents entropy generation, which should always be greater than or equal to zero, i.e.,  ′′′ sgen dV ≥ 0 (2.34)  V Combining eqs. (2.32) and (2.33) and applying eq. (2.34), one obtains the integral form of the second law of thermodynamics for single phase systems: d dt +  V ρ sdV +  A ρ (Vrel ⋅ n)sdA q′′′ dV = T  q′′ ⋅ n dA − AT  V  (2.35) V  ′′′ sgen dV ≥ 0 If the control volume includes Π phases, the second law of thermodynamics must be obtained by integrating over the two phases separately (Faghri and Zhang, 2006),   dt   k =1 Π d Vk (t ) ρk sk dV + Π  Ak (t ) ρk (Vk ,rel ⋅ n k )sk dA +  Ak (t ) q′′ ⋅ n k k dA − Tk  Vk (t ) Tk ′′′  qk dV   =  k =1 Vk (t )  ′′′ sgen ,k dV +  AI (t )  ′′ sgen , I dA ≥ 0 (2.36) The entropy generation for a control volume including Π phases consists of entropy generation in each phase, plus that in the interfaces. The second law of thermodynamics requires that each of these entropy generations be greater than or equal to zero. 2 2 5 Species The continuity equation states that the total mass for a closed system is constant. For a system containing one phase but more than one component, the total mass Chapter 2 Ge era i ed G ver i g Equati s 99 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press of the system is composed of different species. If the concentrations of each of these species are not uniform, mass transfer occurs in a way that makes the concentrations more uniform. Therefore, it is necessary to track the individual components by applying the principle of conservation of species mass. For a system of multiple components, each component can have its own mass density and velocity. The conservation of mass for the ith species in a single-phase system for a fixed-mass system is obtained by applying eq. (2.3) for the ith species with Φ = mi and φ = ρi / ρ (where ρi is the concentration of the ith component), i.e., dmi dt = system ∂ ∂t  V ρi dV +  A ρi (Vi ,rel ⋅ n)dA (2.37) If there is no chemical reaction, the total mass of the ith species for a closed system remains constant. Chemical reactions, on the other hand, will result in the production or consumption of the ith species, which can be modeled as a mass source or sink for the ith species, i.e., dmi dt = system  V  dV mi′′′ (2.38)  where mi'" is the mass production rate for the ith species, which can be either positive (mass production) or negative (mass consumption). For the cases  without chemical reaction, mi'" = 0 . The conservation of species mass for a control volume is therefore ∂ ∂t  V ρi dV +  A ρi (Vi ,rel ⋅ n)dA =  V  dV mi′′′ (2.39) which is valid for each component in the control volume. If the total number of species is N, summation of the conservation of mass of all species results in ∂ ∂t  N   ρi  dV + V  i =1    N   ( ρi Vi ,rel ) ⋅ n  dA = A   i =1   N    mi′′′ dV V   i =1  (2.40) The summation of the densities of the individual species is equal to the bulk density of the multicomponent substance: ρ= ρ i =1 N N i (2.41) The bulk velocity of the multicomponent substance is the mass-averaged velocity of the velocity of each individual species: ρ Vrel =  (ρ V i i =1 i ,rel ) (2.42) The right-hand side of eq. (2.40) must be zero, because the total mass of species produced is equal to the total mass of species consumed, i.e.,   m′′′ = 0 i i =1 N (2.43) Substituting eqs. (2.41) – (2.43) into eq. (2.40), the continuity equation (2.5) is obtained. This means that we can write conservation of species mass for 100 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press N species, but only N-1 of these equations are independent. In practice, one can use these N-1 equations for conservation of species mass in combination with the continuity equation, eq. (2.5), to describe multicomponent systems. In each equation of the conservation of species mass, there are two new unknowns: mass density of the ith species ρi, and the species mass velocity Vi ,rel . This yields more unknowns than the number of equations. So in order to have a properly-posed equation set, it is necessary to reduce the number of unknowns in each equation from two to one. The second term on the left-hand side of eq. (2.39) represents the species i mass flow across the surface of the control volume, which results from convection by bulk flow and diffusion relative to the bulk convection. ρi (Vi ,rel ⋅ n)dA = ρi (Vrel ⋅ n)dA + J i ⋅ ndA (2.44)  A  A  A where J i is the diffusive mass flux vector of species i, which includes mass fluxes due to ordinary diffusion driven by the concentration gradient, and thermal (Soret) diffusion (see Table 1.7). Substituting eq. (2.44) into eq. (2.39), one obtains an expression for the conservation of species mass that contains only one new additional variable, ρi: ∂ ∂t  V ρi dV +  A ρi (Vrel ⋅ n)dA = −  A J i ⋅ ndA +  V  dV mi′′′ (2.45) The above analysis is based on the assumption that the control volume contains only one phase. For a control volume containing Π phases, the conservation of species mass is (Faghri and Zhang, 2006).   ∂t     = − k =1  k =1 Π Π ∂ Vk (t ) ρk ,i dV +  Ak (t ) ρk ,i (Vk ,rel ⋅ nk )dA       ′′′ (J k ,i ⋅ n k )dA + mk ,i dV  Ak (t ) Vk (t )  (2.46)  23 icr sc pic (Differe tia ) ca I sta ce F r u ati The microscopic (differential) formulations to be presented here include conservation equations and jump conditions. The former apply within a particular phase, and the latter are valid at the interface that separates two phases. The phase equations for a particular phase should be the same as those for a singlephase system. Most textbooks (e.g., White, 1991; Incropera and DeWitt, 2001; Bejan, 2004; Kays et al., 2004) obtain the governing equations for a single-phase system by performing mass, momentum, and energy balances for a microscopic control volume. We will obtain the conservation equations by using the integral equations for a finite control volume that includes only one phase. Jump conditions at the interface will be obtained by applying the conservation laws at the interfaces. Chapter 2 Ge era i ed G ver i g Equati s 101 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press In order to obtain the differential formulations, it is necessary to apply the divergence theorem from vector calculus. For a general vector quantity, Ω , which is continuously differentiable, and for a control volume V, enclosed by a piecewise smooth control surface A, the divergence theorem states that Ω ⋅ ndA = ∇ ⋅ ΩdV (2.47)  A  V The desired differential equations may be obtained by applying this relationship to the integral form of the basic laws. Furthermore, since the control volume shape and size are fixed in time, Leibniz’s rule for the specific general quantity, φ , is also valid: d ∂ ( ρφ ) dV (2.48) ρφdV = dt  V  V ∂t 231C servati f ass The surface integral in eq. (2.5) may be converted to a volume integral by applying eq. (2.47) as follows: ρ (Vrel ⋅ n)dA = ∇ ⋅ ρ Vrel dV (2.49)  A  V Substituting eq. (2.49) into eq. (2.5) and considering eq. (2.48), the entire left-hand side of eq. (2.5) is included in a single volume integral, i.e.,  ∂ρ  + ∇ ⋅ ρ Vrel  dV = 0 (2.50)   V  ∂t  The only condition that ensures that eq. (2.50) is true for any arbitrary shape and size of the control volume is that the integrand must be equal to zero, i.e., ∂ρ + ∇ ⋅ ρ Vrel = 0 (2.51) ∂t which is the differential form of the law of the conservation of mass. Equation (2.51) can also be rewritten as Dρ + ρ∇ ⋅ Vrel = 0 Dt (2.52) where D/Dt is the substantial derivative (or material derivative) defined by D ∂ ≡ + Vrel ⋅∇ (2.53) Dt ∂t For a stationary reference frame, Vrel = V , and eq. (2.52) becomes Dρ + ρ∇ ⋅ V = 0 Dt (2.54) For incompressible flow, in which the density of the fluid is constant ( ρ ≡ const ), eq. (2.54) simplifies as ∇⋅V = 0 (2.55) For steady-state compressible flow, the continuity equation can be obtained by simplifying eq. (2.51), i.e., 102 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∇ ⋅ ρ Vrel = 0 (2.56) Example 2.2 For a stationary reference frame, present the continuity equation in the Cartesian coordinate system for both compressible and incompressible fluid flow. Solution: The ∇ operator in the Cartesian coordinate system is [see eq. (G.6) in Appendix G] ∇=i ∂ ∂ ∂ +j +k ∂x ∂y ∂z (2.57) The velocity vector in a Cartesian system is (2.58) where u, v , and w are the velocity components in the x-, y-, and zdirections, respectively. The continuity equation for compressible fluid flow in a stationary reference frame can be obtained by substituting eqs. (2.57) and (2.58) into eq. (2.54), i.e.,  ∂u ∂v ∂w  Dρ (2.59) +ρ + + =0 V = iu + jv + kw Dt  ∂x ∂y ∂z  where D ρ ∂ρ ∂ρ ∂ρ ∂ρ = +u +v +w ∂t ∂x ∂y ∂z Dt (2.60) For incompressible flow, eq. (2.59) simplifies as ∂u ∂v ∂w + + =0 ∂x ∂y ∂z (2.61) 232 e tu Equati The integral form of Newton’s second law for a control volume that includes only one phase is expressed by eq. (2.10). The surface integral terms in eq. (2.10) can be rewritten using the divergence theorem: ′ ′ τ rel ⋅ ndA = ∇ ⋅ τ rel dV (2.62)  A   A ρ (Vrel ⋅ n)Vrel  dA =  V V ∇ ⋅ ρ Vrel Vrel dV (2.63) Substituting eq. (2.62) and (2.63) into eq. (2.10), and considering eq. (2.48), the entire equation can be rewritten as a volume integral:  ′ ∇ ⋅ τ rel + V  ρ X i i =1 N i −  ∂ ( ρk Vrel ) − ∇ ⋅ ρ Vrel Vrel  dV = 0 ∂t   (2.64) As was the case for the continuity equation, the integrand must equal zero to assure the general validity of eq. (2.64); so, one obtains the desired differential form of the momentum equation: Chapter 2 Ge era i ed G ver i g Equati s 103 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂ ′ ( ρ Vrel ) + ∇ ⋅ ρ Vrel Vrel = ∇ ⋅ τ rel + ∂t ρ X i i =1 N i (2.65) The derivatives on the left-hand side of eq. (2.65) may be expanded to yield N  ∂V   ∂ρ  ′ Vrel  + ∇ ⋅ ρ Vrel  + ρ  rel + Vrel ⋅∇Vrel  = ∇ ⋅ τ rel + ρi Xi (2.66)  ∂t   ∂t   i =1 The first bracketed term on the left vanishes, as required by the continuity eq. (2.51). The second term may be written more simply in substantial derivative form, and the entire equation becomes ρ DVrel ′ = ∇ ⋅ τ rel + Dt ρ X i i =1 N i (2.67) ′ The stress tensor, τ rel , is the sum of an isotropic thermodynamic stress, − pI , and the viscous stress tensor, τ rel [defined in eq. (1.56)], i.e., τ ′ = − pI + τ rel (2.68) rel Substituting eq. (2.68) into eq. (2.67), the momentum equation becomes ρ DVrel = −∇p + ∇ ⋅ τ rel + Dt ρ X i i =1 N i (2.69) The viscous stress tensor measured in the reference frame, τrel, can be determined by using Newton’s law of viscosity [see eq. (1.59)]: 2 τ rel = 2μ Drel − μ (∇ ⋅ Vrel )I 3 (2.70) where Drel is the rate of strain tensor, i.e., Drel = 1 T ∇Vrel + ( ∇Vrel )   2 (2.71) and I in eq. (2.70) is the unit tensor that satisfies a·I = I·a = a for any tensor a. The diagonal components of I are equal to one and all other components are zero: 1 i = j I ij =  0 i ≠ j (i , j = 1, 2,3) (2.72) If the fluid is incompressible ( ρ =const ), the second term on the right-hand side of eq. (2.70) will be zero according to eq. (2.55). The momentum equation (2.67) then becomes ρ DVrel = Dt ρ X i i =1 N i − ∇p + ∇ ⋅ ( μ∇Vrel ) (2.73) where the left-hand side is the inertial term (mass per unit volume ρ times acceleration, DVrel/Dt). The three terms on the right-hand side represent body force per unit volume, pressure force per unit volume, and viscous force per unit volume, respectively. For DVrel / Dt = 0 , we have Stokes’ flow or creep flow, and eq. (2.73) becomes elliptic and is similar to the steady-state conduction equation. 104 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press In a Cartesian coordinate system, the vector form of the momentum equation, eq. (2.73), for incompressible and Newtonian fluid with constant viscosity can be written as three equations in the x-, y-, and z-directions: ρ ρ ρ Durel = Dt Dvrel = Dt ρ X i i =1 N i =1 N N i −  ∂ 2 urel ∂ 2 urel ∂ 2 urel ∂p +μ + +  ∂x 2 ∂x ∂y 2 ∂z 2  ∂p  ∂ 2 vrel ∂x 2     (2.74) (2.75)      ρ Y − ∂y + μ    ii + ∂ 2 vrel ∂y 2 + ∂ 2 vrel   ∂z 2   Dwrel = Dt  i =1 ρi Z i −  ∂ 2 wrel ∂ 2 wrel ∂ 2 wrel ∂p +μ + +  ∂x 2 ∂z ∂y 2 ∂z 2  (2.76) where X i , Yi , and Z i are the components of body force per unit volume acting on the ith species in the x-, y-, and z- directions, respectively. For the case that the only body force is gravity, Xi = g , eq. (2.73) becomes ρ DVrel = ρ g − ∇p + ∇ ⋅ ( μ∇Vrel ) Dt (2.77) For natural convection problem, it is often assumed that the fluid is incompressible except in the first term on the right-hand side of eq. (2.77); this is referred to as the Boussinesq assumption. The density of a mixture is a function of temperature and mass fractions of species. It can be expanded using a Taylor’s series near the vicinity of a reference point ( T , ω1 , ω2 , ω N ): ρ=ρ+ ∂ρ (T − T ) + ∂T i =1 N  ∂ω N ∂ρ i (ωi −ωi ) +  ≈ ρ − ρβ (T − T ) − ρ β i =1 (2.78) −ωi ) m (ωi where ρ is density at the reference point, β = −(∂ρ / ∂T ) / ρ is the coefficient of thermal expansion, and β m = −(∂ρ / ∂ωi ) / ρ is the composition coefficient of volume expansion. Substituting eq. (2.78) into eq. (2.77), the momentum equation for natural convection is obtained ρ N DVrel = ( −∇p + ρ g ) − ρ gβ (T − T ) − ρ g β m (ωi −ωi ) + ∇ ⋅ ( μ∇Vrel ) Dt i =1  (2.79) where the second and third terms on the right-hand side of eq. (2.79) describe the effect of buoyancy force due to temperature and composition variation within the system, respectively. 2 3 3 E ergy Equati The surface integrals on the right-hand side of eq. (2.20) can be rewritten as volume integrals using eq. (2.47), i.e., Chapter 2 Ge era i ed G ver i g Equati s 105 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  A ρ (Vrel ⋅ n)  e +    2 Vrel 2   dA =      V2 ∇ ⋅  ρ Vrel  e + rel  V 2    V    dV   (2.80) (2.81) (2.82)    A A −q′′ ⋅ ndA = ′ n ⋅ τ rel ⋅ Vrel  −∇ ⋅ q′′dV dA =  ∇ ⋅ ( τ ′ ⋅V V rel rel )dV Substituting eqs. (2.80) – (2.82) into eq. (2.20) and considering eq. (2.48) yields 2 ∂   Vrel   ρ  e + V ∂t   2       V2   + ∇ ⋅  ρ Vrel  e + rel   2          +∇ ⋅ q′′ − q′′′ − ∇ ⋅ ( τ ′ rel ⋅Vrel ) −        ρi Xi  ⋅ Vrel  dV = 0   i =1   N (2.83) For eq. (2.83) to be true for any arbitrary control volume, the integrand must be zero. The result is the general differential form of the energy equation: V2 ∂   ρ  e + rel ∂t   2     V2   + ∇ ⋅  ρ Vrel  e + rel   2          = −∇ ⋅ q′′ + q′′′ + ∇ ⋅ ( τ ′ rel ⋅Vrel ) +     V2  e + rel  2    ρi Xi  ⋅ Vrel  i =1  N (2.84) The left-hand side of eq. (2.84) can be rewritten as ∂    ∂ρ V2    + ∇ ⋅ ( ρ Vrel )  + ρ   e + rel   ∂t 2   ∂t    N i   V2  + Vrel ⋅∇  e + rel   2    i rel      = −∇ ⋅ q′′ + q′′′ + ∇ ⋅ ( τ ′ rel ⋅Vrel ) +     ρ X ⋅V   i =1 N (2.85) According to the continuity eq. (2.51), the first bracketed term on the lefthand side of eq. (2.85) is zero. The second term on the left-hand side may be written more simply in substantial derivative form, i.e., ρ V2 D  e + rel Dt  2     = −∇ ⋅ q′′ + q′′′ + ∇ ⋅ ( τ ′ rel ⋅Vrel ) +      ρ X ⋅V   i i i =1  rel (2.86) which is the total energy (including thermal and mechanical energies) balance equation. It would be more convenient to remove the mechanical energy terms from eq. (2.86). The mechanical energy balance equation can be obtained by forming a dot (scalar) product of the momentum equation (2.67) with the velocity vector Vrel, i.e., ρ  DVrel ′ ⋅ Vrel = ( ∇ ⋅ τ rel ) ⋅ Vrel +   Dt  ρ X ⋅V   i i i =1 N  rel (2.87) 106 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press which can be rearranged to obtain ρ D  Vrel 2  Dt  2     = ∇ ⋅ ( τ ′ ⋅ Vrel ) − ∇V : τ rel +  rel      ρ X ⋅V   i i i =1 N  rel (2.88) Subtracting eq. (2.88) from eq. (2.86) yields the following thermal energy equation: ρ De = −∇ ⋅ q′′ + q′′′ − p∇ ⋅ Vrel + ∇Vrel : τ rel Dt (2.89) where the left-hand side represents the rate of gain of energy per unit volume. The terms in the right-hand side are the rate of energy input by heat transfer per unit volume, the internal heat generation per unit volume, the reversible rate of energy increase per unit volume by compression, and the irreversible rate of energy increase per unit volume by viscous dissipation, respectively; the viscous dissipation is the contraction [see eq. (G.35)] of two tensors of rank two: ∇Vrel and τ rel , i.e., ∇Vrel : τ rel = ∂urel ∂u ∂u ∂v ∂v τ xx + rel τ xy + rel τ xz + rel τ yx + rel τ yy ∂x ∂y ∂z ∂x ∂y (2.90) ∂v ∂wrel ∂w ∂w τ zx + rel τ zy + rel τ zz + rel τ yz + ∂z ∂x ∂y ∂z Equation (2.89) is the energy equation expressed in terms of internal energy. In order to obtain an equation that contains enthalpy, the definition of enthalpy is employed: (2.91) ρ Substituting eq. (2.91) into eq. (2.89), and considering continuity equation (2.51), the energy equation in terms of enthalpy and temperature is obtained: h =e+ p ρ Dh Dp = −∇ ⋅ q′′ + + q′′′ + ∇Vrel : τ rel Dt Dt h= (2.92) For a multicomponent system, the enthalpy can be expressed as ω h i =1 N ii (2.93) where ωi and hi are the mass fraction and specific enthalpy of the ith component. Substituting eq. (2.93) into eq. (2.92), the energy equation for a multicomponent system becomes   ρ Dt (ω h ) = −∇ ⋅ q′′ + Dt + q′′′ + ∇V     D Dp ii i =1 N rel : τ rel (2.94) The heat flux vector q′′ can be obtained from eq. (1.70) q′′ = −k∇T +  i =1 N hi J i + cRu T  i =1 N xi x j DiT ρi Dij  Ji J j  −   ρi ρ j    (2.95) Chapter 2 Ge era i ed G ver i g Equati s 107 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Substituting eq. (2.95) into eq. (2.94), the energy equation becomes  D   ρ Dt (ωi hi )  = ∇ ⋅ (k∇T ) −   i =1 N ∇ ⋅(J h ) ii i =1 N  −∇ ⋅  cRu T    i =1 N xi x j DiT ρi Dij  J i J j   Dp + q′′′ + ∇Vrel : τ rel −  +  ρi ρ j   Dt   (2.96) As noted in Chapter 1, the second and third terms on the right hand side, which represent the contributions of interdiffusional convection and Dufour effect, are both negligible for most applications. For a pure substance, eq. (2.96) is reduced to: ρ Dh Dp = ∇ ⋅ (k ∇T ) + + q′′′ + ∇Vrel : τ rel Dt Dt (2.97) The enthalpy of a pure substance (single component) can be expressed as a function of temperature and pressure, h = f (T , p), i.e., Dh  ∂h  DT  ∂h  Dp = +   Dt  ∂T  p Dt  ∂p T Dt (2.98) Thermodynamic relations give us  ∂h  1− βT  ∂h    = cp ,   = ρ  ∂T  p  ∂p T (2.99) where β= 1  ∂ρ  ρ  ∂T  p   (2.100) is the coefficient of thermal expansion. Therefore, eq. (2.98) becomes Dh DT 1 − β T Dp (2.101) = cp + Dt Dt ρ Dt Substituting eq. (2.101) into eq. (2.97), the energy equation becomes ρcp DT Dp = ∇ ⋅ ( k ∇T ) + T β + q′′′ + ∇Vrel : τ rel Dt Dt (2.102) which can be simplified for different substances as demonstrated below. For ideal gases, substituting the equation of state for ideal gas [ ρ = p / ( RgT ) ] into eq. (2.100) yields β = 1 / T , and eq. (2.102) becomes ρcp DT Dp = ∇ ⋅ (k∇T ) + + q′′′ + ∇Vrel : τ rel Dt Dt (2.103) For incompressible fluids, the constant density yields β = 0 and eq. (2.102) becomes ρcp DT = ∇ ⋅ (k∇T ) + q′′′ + ∇Vrel : τ rel Dt (2.104) For solids, the density is constant and the velocity is zero, and the energy equation becomes 108 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ρ s c ps ∂Ts ′′′ = ∇ ⋅ (ks ∇Ts ) + qs ∂t (2.105) which is a heat conduction equation. Example 2.3: Determine the viscous dissipation for incompressible flow in the Cartesian coordinate system. The reference velocity can be assumed to be zero. Solution: Since the reference velocity is zero, the viscous dissipation is ∇V : τ . According to eq. (G.33),  ∂u  ∂x   ∂v ∇V =   ∂x  ∂w   ∂x ∂u ∂y ∂v ∂y ∂w ∂y ∂u  ∂z   ∂v   ∂z  ∂w   ∂z  The shear stress tensor for incompressible flow is   ∂u ∂v  ∂u  ∂u ∂w   μ +  μ +  2μ  ∂x  ∂z ∂x    ∂y ∂x     ∂v ∂u   ∂v ∂w   ∂v τ = 2μ D =  μ  + 2μ μ +   ∂y   ∂x ∂y   ∂z ∂y     ∂w   μ  ∂w + ∂u  μ  ∂w + ∂v  2μ     ∂x ∂z   ∂z   ∂y ∂z    (2.106) Therefore, the viscous dissipation is  ∂u  ∂x   ∂v ∇V : τ =   ∂x  ∂w   ∂x ∂u ∂y ∂v ∂y ∂w ∂y  ∂u ∂v   ∂u ∂w   ∂u   2μ ∂u μ +  μ +   ∂x  ∂z ∂x    ∂y ∂x  ∂z     ∂v ∂w   ∂v    ∂v ∂u  ∂v :μ + 2μ μ +    ∂z    ∂x ∂y  ∂y  ∂z ∂y    ∂w    ∂w ∂u   ∂w ∂v  ∂w   μ  μ 2μ + +  ∂z    ∂x ∂z  ∂z   ∂y ∂z    The contraction of two tensors of rank two yields 2  ∂u  2 ∂u ∂v   ∂u 2 ∂u ∂w   ∂u   + μ   + ∇V : τ = 2 μ   + μ    +  ∂z ∂x   ∂y  ∂y ∂x   ∂x   ∂z      2  ∂v 2 ∂v ∂u   ∂v 2 ∂v ∂w   ∂v  + μ   +  + 2 μ   + μ   +   ∂x  ∂x ∂y   ∂y   ∂z  ∂z ∂y      2  ∂w  2 ∂w ∂v   ∂w 2 ∂w ∂u   ∂w   + 2μ  + μ  +  + μ  +   ∂x ∂z  ∂y ∂z   ∂y   ∂z   ∂x      (2.107) Chapter 2 Ge era i ed G ver i g Equati s 109 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The final form of the viscous dissipation for incompressible flow is obtained by rearranging eq. (2.107), i.e., 2 2   ∂u  2  ∂v   ∂w  ∇V : τ = μ  2   + 2   + 2    ∂x  ∂z   ∂y     ∂v ∂u   ∂w ∂v   ∂u ∂w  + + +  + +  +   ∂x ∂y   ∂y ∂z   ∂z ∂x  2 2 2 (2.108)    The viscous dissipation in cylindrical and spherical coordinate systems can be found by similar means using information in Appendix G. 2 3 4 The Sec d aw f Ther dy a ics To obtain the differential form of the second law of thermodynamics, the surface integrals in eq. (2.35) can be rewritten as volume integrals: ρ (Vrel ⋅ n)sdA = ∇ ⋅ ( ρ Vrel s )dV (2.109)  A  V  q′′ ⋅ ndA = AT   q′′  ∇ ⋅   dV V T  (2.110) Substituting eqs. (2.109) and (2.110) into eq. (2.35) and considering eq. (2.48) yields  ∂  q′′  q′′′   ∂t ( ρ s ) + ∇ ⋅ ( ρ Vrel s ) + ∇ ⋅  T  − T  dV = V    V ′′′ sgen dV ≥ 0 (2.111) In order for eq. (2.111) to be true for any arbitrary control volume, the integrand in eq. (2.111) should always be positive, i.e., ∂ q′′′ q′′  ′′′ ( ρ s ) + ∇ ⋅ ( ρ Vrel s) + ∇ ⋅   − = sgen ≥ 0  ∂t T  T (2.112) Equation (2.112) can be rewritten as  ∂ρ   ∂s   q′′  q′′′  ′′′ = sgen ≥ 0 s  + ∇ ⋅ ( ρ Vrel )  + ρ  + Vrel ⋅∇s  + ∇ ⋅   −  ∂t   ∂t  T  T (2.113) Considering the continuity equation, eq. (2.51), and definition of the substantial derivative, eq. (2.53), eq. (2.113) can be reduced to ρ Ds  q′′  q′′′  ′′′ + ∇ ⋅  − = sgen ≥ 0 Dt T  T (2.114) where the three terms on the left-hand side represent rate of change of entropy per unit volume, rate of change of entropy per unit volume by heat transfer and internal heat generation, respectively. Equation (2.114) means that the entropy generation per unit volume must not be negative at any time or location. For a multicomponent system without internal heat generation ( q′′′ = 0 ), Curtiss and Bird (1999; 2001) obtained the entropy flux vector and the entropy generation as 110 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press s′′ =   ′′′ Tsgen = −  q′′ −   1  q′′ − T  N h J    ii i =1 N  N (2.115)   hi J i  ⋅∇ ln T −  i =1  N   cRu T  d  − τ : ∇V −  Ji ⋅ ρi i  i =1  M i =1 gi i  mi′′′ (2.116) where di is the diffusional driving force [see eq. (1.115) – (1.117)], gi is partial molar Gibbs free energy, q′′ is total heat flux, obtained by eq. (2.95), including conduction, interdiffusional convection, and Dufour effect. 2 3 5 Species The surface integrals in eq. (2.45) may be converted to volume integrals by applying eq. (2.47) as follows: ρi (Vrel ⋅ n)dA = ∇ ⋅ ρi Vrel dV (2.117)  A   J ⋅ ndA =  ∇ ⋅ J dV V A i V i (2.118) Substituting eqs. (2.117) and (2.118) into eq. (2.45) and considering eq. (2.48), the entire left-hand side of eq. (2.45) is included in a single volume integral, i.e.,  ∂ρi  + ∇ ⋅ ρi Vrel + ∇ ⋅ J i − mi′′′ dV = 0 (2.119)   V  ∂t  The only condition that makes eq. (2.119) true regardless of the shape and size of the control volume is that the integrand must equal zero, i.e., ∂ρi  (2.120) + ∇ ⋅ ρi Vrel = −∇ ⋅ J i + mi′′′, i = 1, 2, , N ∂t The first term on the left-hand side is the rate of increase of mass of the species i per unit volume, and the second term is net rate of additions of mass of the ith species per unit volume by convection. The terms on the right-hand side are the net rate of mass of ith species per unit volume by diffusion, and rate of production of species i by chemical reaction. Equation (2.120) is the equation of conservation of mass for species. If the total number of species is N, a total of N-1 independent equations for conservation of species mass can be obtained. After defining the mass fraction of species i as ρ ωi = i (2.121) ρ eq. (2.120) can be rewritten as  ∂ωi   ∂ρ   + ∇ ⋅ ωi Vrel  = −∇ ⋅ J i + mi′′′ (2.122)  + ∇ ⋅ ρ Vrel  ωi + ρ   ∂t   ∂t  According to eq. (2.51), the first bracket on the left-hand side of eq. (2.122) is zero. The second bracket on the left-hand side is the substantial derivative of Chapter 2 Ge era i ed G ver i g Equati s 111 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the mass fraction. Therefore, the conservation of species mass in terms of mass fraction becomes Dωi  ρ = −∇ ⋅ J i + mi′′′ (2.123) Dt Assuming binary system of A and B, one can use Fick’s law in eq. (2.123) to yield Dω A A ρ = − ρ∇ ⋅ ( D AB ∇ω A ) + m ′′′ (2.124) Dt which is useful in determining the diffusion in dilute liquid solution at constant A temperature and pressure. Equation (2.124), with m ′′′ = 0 , is similar to energy equation (2.104) with no internal heat source and viscous dissipation and therefore it is used for analogy between heat and mass transfer analysis. In the preceding discussion to develop eq. (2.120), the mass fraction and mass flux were used. The species equation can also be developed in terms of molar concentration (or molar fraction) and molar flux. By following a similar procedure, the species equation is ∂ci i  = −∇ ⋅ n′′ + ni′′′ ∂t (2.125) where the molar flux relative to the stationary coordinate axes can be obtained from eq. (1.104), i.e.,  i n′′ = ci V* + J* (2.126) i Substituting eq. (2.126) into eq. (2.125), we have ∂ci   + ∇ ⋅ (ci V* ) = −∇ ⋅ J* + ni′′′ i ∂t (2.127) where the first term on the left-hand side is rate of increase of mole of the species i per unit volume, and the second term is net rate of additions of mole of the ith species per unit volume by convection. The terms on the right-hand side are net rate of mole of ith species per unit volume by diffusion, and the molar rate of production of species i by chemical reaction. For a binary system of components A and B with constant pressure, eq. (2.127) reduces to ∂c A  A + ∇ ⋅ (c A V* ) = −c∇ ⋅ ( D AB ∇x A ) + n ′′′ (2.128) ∂t  where c is the mixture molar concentration and V* is molar-averaged velocity defined in Chapter 1. Equation (2.128) can be applied to low density gases with constant temperature and pressure.   The production rate of the ith species, mi′′′ (or ni′′′ ), is still unknown at this point, but it can be obtained by analyzing the chemical reaction. If the number of chemical reactions taking place in the system is Nc, the mass production rate is (Kleijn, 1991; Mahajan, 1996)  mi′′′ = a j =1 Nc ij M i ℜ j (2.129) 112 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where aij is the stoichiometric coefficient that describes the proportions of the mole numbers of reactants disappearing and mole numbers of products appearing as a result of the reaction process (see Section 2.3.3). The net reaction rate of the jth chemical reaction ℜ j is the difference between the forward and backward reactions, i.e., ℜ j = ℜ j+ − ℜ j− (2.130) If chemical reactions take place in gas mixture, the forward and backward reaction rates are ℜ j + = rj + ( p, T ) ℜ j − = r j − ( p, T ) ∏ i =1 Np i =1 Nr p  xi    Ru T  p  xi    RuT  aij (2.131) aij ∏ (2.132) where Nr and Np are number of reactants and products, respectively. The two reaction rate constants, rj + and rj − , depend on the specific chemical reaction under consideration. These constants are related by rj − ( p, T ) rj + ( p, T ) = 1   K j (T )  p0    a  Ru T  ij i =1 N (2.133) where p0 is the standard pressure and Kj(T) is the thermodynamic equilibrium constant of the jth chemical reaction:  −ΔG 0 (T )  j K j (T ) = exp   Ru T     (2.134) where ΔG 0 (T ) is the standard Gibbs energy for the jth chemical reaction; its j value depends on the specific chemical reaction considered. For the case in which only one chemical reaction ( N c = 1 ) leads to production of the ith component, the molar production rate can be obtained by a simplified from  ′′′ ni′′′ = kn cin (2.135) ′ where the index n represents the order of the reaction and kn′′ is a rate constant (1/sn) depends on the temperature. As before, the mass flux, J i , in eqs. (2.120) and (2.123) includes mass fluxes due to ordinary diffusion driven by the concentration gradient, pressure diffusion, body force diffusion, and thermal (Soret) diffusion [see eq. (1.117)]. For a binary mixture, the mass flux can be calculated by eq. (1.132). Equation (2.123) is valid for the case that chemical reaction occurs in the entire volume – referred to as homogeneous reacting system. For the case where the chemical reaction takes place on a surface – referred to as heterogeneous reaction – the source term in eq. (2.123) will not appear, and the rate of production will be accounted as a boundary condition. More discussion about Chapter 2 Ge era i ed G ver i g Equati s 113 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press modeling of the heterogeneous reacting system can be found in Faghri and Zhang (2006). Example 2.4 Couette flow is found in many engineering applications, such as lubricant flow in a journal bearing. Consider the laminar viscous Couette flow where two plates are of infinite extent, steady state conditions apply (see Fig. 2.4). The fluid is Newtonian and incompressible. Starting from the basic differential equations and using mathematical, physical arguments to make suitable simplifying assumptions, derive the differential equations governing the fluid velocity and temperature, as well as boundary conditions. Obtain the velocity and temperature distributions. The fluid properties can be assumed to be constant, but viscous dissipation should be accounted for. Solution: To solve the problem, the following assumptions are made: 1. The flow and heat transfer are steady state, ∂ / ∂t = 0 2. Thermophysical properties are constants 3. The flow is laminar 4. 5. 6. 7. No pressure gradient in the x-direction, ∂p / ∂x = 0 No edge effect in the z-direction, ∂ / ∂z = 0 No end effect in the x-direction, ∂ / ∂x = 0 Gravity is the only body force, X = Z = 0, Y = − g If the reference frame is stationary, the continuity, momentum, and energy equations can be simplified based on the above assumptions ∂u ∂v ∂w + + =0 ∂x ∂y ∂z (2.136) 0 Assumptions (6) 0 (5) Moving plate T2 Plate 2 y, v z, w Plate 1 x, u Figure 2 4 Couette Flow. U L T1 Stationary plate 114 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ρ  ∂ 2u ∂ 2u ∂ 2u   ∂u ∂u ∂u ∂u  ∂p +u +v +w  = ρX − +μ 2 + 2 + 2   ∂x ∂x ∂y ∂z  ∂x ∂y ∂z   ∂t   (2.137) 0 Assumptions (1) ρ 0 (6) 0 (5) 0 (7) 0 (4) 0 (6) 0 (5) (2.138)  ∂2v ∂2v ∂2v   ∂v ∂v ∂v ∂v  ∂p + u + v + w  = ρY − +μ 2 + 2 + 2   ∂x ∂x ∂y ∂z  ∂y ∂y ∂z   ∂t   0 0 Assumptions (1) (6) ρ 0 (5) −ρ g 0 (6) (7) 0 (5)  ∂2w ∂2w ∂2w   ∂w ∂w ∂w ∂w  ∂p +u +v +w +μ 2 + 2 + 2   = ρZ −   ∂x ∂y ∂z  ∂z ∂y ∂z   ∂t  ∂x 0 Assumptions (1) ρcp  0 (6) 0 (5) 0 (7) 0 (5) 0 (6) 0 (5) (2.139)  ∂ 2T ∂ 2T ∂ 2T  ∂T ∂T ∂T ∂T  +u +v +w  = k 2 + 2 + 2  ∂x ∂x ∂y ∂z  ∂y ∂z  ∂t    ∂u 2   + μ 2      ∂x    Assumptions 0 (1) 2 0 (6) 2 0 (5) 2 0 (6) 2 0 (5) 0 (6) 2  ∂v   ∂w   ∂v ∂u   ∂w ∂v   ∂u ∂w  + 2  + 2 +  + +  +  + +   ∂z   ∂x ∂y   ∂y ∂z   ∂z ∂x   ∂y   (2.140)   0 Assumptions (5) 0 (6) ∂v =0 ∂y 0 (5) 00 (5) (6) Therefore, eqs. (2.136) - (2.140) can become (2.141) (2.142) (2.143) (2.144) (2.145) ρv ρv ∂u ∂2u =μ 2 ∂y ∂y ∂v ∂p ∂2v = −ρ g − +μ 2 ∂y ∂y ∂y ρv ρcpv ∂w ∂2w =μ 2 ∂y ∂y   ∂v  2  ∂u  2  ∂w 2  ∂T ∂ 2T = k 2 + μ 2   +   +   ∂y ∂y   ∂y   ∂y   ∂y     The boundary conditions for eqs. (2.141) - (2.145) are u = w = 0, y = 0 (no-slip) v = 0, y = 0 (impermeable wall) T = T1 , y = 0 (no temperature jump) u = U , w = 0, y = L (no-slip) (2.146) (2.147) (2.148) (2.149) Chapter 2 Ge era i ed G ver i g Equati s 115 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (2.150) (2.151) The solution of continuity equation (2.141) with boundary conditions specified by eq. (2.147) and (2.150) is v=0 (2.152) Substituting eq. (2.152) into eq. (2.144), one obtains v = 0, y = L (impermeable wall) T = T2 , y = L (no temperature jump) ∂2w ∂y 2 =0 (2.153) The solution of eq. (2.153) with boundary conditions specified by eqs. (2.146) and (2.149) is w=0 (2.154) Substituting eqs. (2.152) and (2.154) into eqs. (2.142), (2.143) and (2.145), one obtains ∂2u ∂y 2 =0 (2.155) (2.156) 2 ∂p = −ρ g ∂y k ∂ 2T  ∂u  = −μ   2 ∂y  ∂y  (2.157) The solutions of eqs. (2.155) with boundary conditions specified by eqs. (2.146) and (2.149) is u= y U L (2.158) The solution of eq. (2.156) is p = p0 − ρ gy (2.159) where p0 is pressure at y = 0. Substituting eq. (2.158) into eq. (2.157) and considering eqs. (2.148) and (2.151), the solution for temperature distribution is obtained T − T1 y Ec Pr y  y  =+ 1 −  2 L L T2 − T1 L (2.160) where Ec = U 2 / [c p (T2 − T1 )] is Eckert number. Example 2.5 Simplify the energy equation for a multicomponent system, eq. (2.96), for the case in which the Dufour effect is negligible but the interdiffusion term needs to be considered. It is assumed that the reference frame is stationary (Vref = 0) and the viscous dissipation is negligible. Solution: With the Dufour effect and viscous dissipation neglected, eq. (2.96) becomes 116 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  or D   ρ Dt (ωi hi )  = ∇ ⋅ (k∇T ) −   i =1 N  ∇ ⋅ ( J h ) + Dt ii i =1 N i =1 N Dp (2.161)  ∂   ∂t (ωi ρ hi ) + ∇ ⋅ ( Vωi ρ hi )  = ∇ ⋅ (k∇T ) −   i =1 N  ∇ ⋅ ( J h ) + Dt ii Dp (2.162) The diffusive mass flux can be rewritten in terms of the velocity of one component, Vi , compared to the mass-averaged velocity, V [see eq. (1.99)] J i = ωi ρ ( Vi − V ) (2.163) The advective and diffusive term can be combined as ωi ρ V + J i = ωi ρ Vi (2.164) The mass averaged velocity, V, is a function of the component velocities V= ω V i i =1 N i (2.165) Using the component velocity, Vi, defined in eq. (2.164), in the energy equation, eq. (2.162), yields:   ∂t (ω ρ h ) + ∇ ⋅ (ω ρ V h ) = ∇ ⋅ (k∇T ) + Dt     i i i ii i =1 N ∂ Dp (2.166) For a stationary reference frame, the species conservation equation, (2.120), becomes ∂  (ωi ρ ) + ∇ ⋅ (ωi ρ V ) = −∇ ⋅ J i + mi′′′ ∂t (2.167) Substituting eq. (2.164) into eq. (2.167), the conservation of species equation can be written in terms of the velocity of the ith component ∂  (ωi ρ ) + ∇ ⋅ (ωi ρ Vi ) = mi′′′ ∂t (2.168) For convenience, let’s define a substantial derivative in terms of the velocity of the ith component Di ( Dt )= ∂ ( ∂t ) + Vi ⋅∇ ( ) (2.169) Using the species equation, eq. (2.168), the energy equation (2.166) can be rewritten as  The specific enthalpy of the ith component is a function of temperature and pressure, hi = hi ( T , p ) . The substantial derivative of enthalpy is [see eq. (2.101)] Di hi D T 1 − β i T Di p = c p,i i + (2.171) ρi Dt Dt Dt Di hi  Dp  ωi ρ Dt  = ∇ ⋅ (k∇T ) + Dt −  i =1  N   [m′′′ h ] i i i =1 N (2.170) Chapter 2 Ge era i ed G ver i g Equati s 117 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where βi = 1  ∂ρi  ρi  ∂T  p   (2.172) is the coefficient of thermal expansion for the ith component in the kth phase. Substituting eq. (2.171) into eq. (2.170) yields:  Di T Di p  Dp  ωi ρ c p,i Dt + (1 − β i T ) Dt  = ∇ ⋅ (k∇T ) + Dt −  i =1  N   [m′′′ h ] i i i =1 N (2.173) Since the mass-averaged velocity is a function of the component velocity, a simplification can be made on one of the pressure terms. Dp − Dt N  i =1 i N Di p = V ⋅ ∇p − Dt  i =1 N Vi ⋅∇p = −  (1 − ω ) V ⋅ ∇p i i i =1 N (2.174) Substituting this relation into the energy equation, eq. (2.173), yields  ω ρ c  i =1 N  p,i Di T  = ∇ ⋅ (k ∇T ) Dt   Di p   − (1 − ωi ) Vi ⋅ ∇p +  β i T Dt  −  i =1 i =1    N   [m′′′ h ] i i i =1 N (2.175) where the second and third terms on the right hand side of the energy equation above are due to the difference in pressure work from the mass average velocity, and the pressure work from each component. The last term is due to species generation or consumption. For a single component system, the energy equation reduces to ρcp DT Dp = ∇ ⋅ ( k ∇T ) + β T Dt Dt (2.176) Example 2.6 Determine the rate at which volatile liquid A evaporates into a stagnant gas B in a tube configuration as shown in Fig. 2.5, where the total pressure and temperature remain constant. Solution: This is the classic Stefan’s diffusion problem. The following assumptions are made to solve this problem: 1. Steady state 2. 1D problem 3. Constant total pressure and temperature 4. Liquid A is permeable to gas B 5. Ideal gas for gas mixture 6. Flat equilibrium condition at the liquid/vapor interface Based on the above six assumptions, the continuity and species equations are reduced to the following form 118 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Gas stream B Gases A&B H y Liquid A Figure 2 5 Evaporation of liquid A into a stagnant gas B. d ( ρ v) =0 dy (2.177) dω ρ v A = D AB dy  dω A  dρ   dy  dy (2.178) where ωA is the mass fraction of species A in the gas mixture and v is the velocity in the y-direction. It should be noted that energy and momentum equations are unnecessary since temperature and pressure are constant. However, mixture density is not constant. It is customary to solve this problem in terms of mole fraction, for simplicity, rather than mass fraction. The appropriate boundary conditions are At vapo Rc = 1 k  s2 s1 ds /liquid interface A( s ) vB = 0 y = 0 , ω A = ω A,δ (2.179) At the top of the tube (2.180) The velocity of the gas mixture, v, is related to the velocity of species, i.e., v = ω A v A + ω B vB (2.181) At the interface, vB = 0 , since liquid is impermeable to gas B. The mixture velocity at the interface is v = ω A vA (2.182) The rate of evaporation of species A at the interface is y = H , ω A = ω A, H Chapter 2 Ge era i ed G ver i g Equati s 119 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press A m ′′ = ρ A v A = ρ A v − ρ D AB dω A dy (2.183) Substituting eq. (2.182) into eq. (2.183) yields the fully related equation for evaporation at the interface dω A − ρ D AB dy A (2.184) m ′′ = ρ A v A = 1− ωA For a one-dimensional, steady-state problem A dm′′ =0 dy (2.185) Substituting eq. (2.184) into eq. (2.185) yields d  − ρ D AB  dω A =0   dy  1 − ω A  dy − ρ D AB dω A A = m′′ = constant 1 − ω A dy (2.186) (2.187) Since the gas mixture is assumed to be an ideal gas and isothermal, the total pressure is the sum of the partial pressures of components A and B p = p A + pB (2.188) where p = ρ RuT / M , p A = ρ A RuT / M A , and pB = ρ B RuT / M B Further simplification reduces eq. (2.188) to ρ1 ρ1 1 =A +B (2.189) M ρ M A ρ MB Using eq. (2.189) and definitions, with some simplification the mass fraction is ωA = p − p A (1 − M A / M B ) pA ( M A / M B ) (2.190) Similarly, the mixture density is given by ρ= M pM =B Ru T Ru T   M A   p − p A 1 −    M B    (2.191) Substituting eqs. (2.190) and (2.191) into eq. (2.187) yields  pM A  D AB  dp A A − = m′′ = constant    Ru T   p − p A  dy (2.192) Integrating eq. (2.192) and applying the boundary conditions yields  p − p A  y  p − p A, H  ln  (2.193) = ln   p− p  H  p− p    A ,δ  A ,δ    where p A,δ and p A, H are the partial pressures of species A at the liquid/vapor interface and at the top of the tube, respectively. Substituting eq. (2.193) into eq. (2.192) yields the rate of evaporation of species A 120 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press A m ′′ = pM A D AB  p − p A, H ln   p− p RuTH A ,δ      (2.194) where temperature T and pressure p are assumed to be uniform within the gas mixture. 2 3 6 C assificati f PDEs a d B u dary C diti s A general transport equation, whether it is mass, momentum, energy, or species, can be written as: ∂ (2.195) ( ρΦ ) + ∇ ⋅ ( ρ VΦ ) = ∇ ⋅ ( Γ∇Φ ) + F ( x, t , Φ,) ∂t where the dependent variable, Φ , is one for mass, any component of velocity ( u, v, w ) , enthalpy, h, or the species mass fraction, ωi . The gradient operator is the partial derivative of the dependent variable with respect to all the spatial directions x for a given coordinate system. A partial differential equation (PDE) is an equation of a function and its partial derivatives. In general, a PDE is classified by its linearity or nonlinearity, and by its order. Its order is considered by its highest derivative. A general second order PDE for two independent variables, η and ζ , is A ∂2Φ ∂η 2 + 2B ∂ 2Φ ∂2Φ ∂Φ ∂Φ +C +D +E + FΦ + G = 0 ∂η∂ζ ∂η ∂ζ ∂ζ 2 (2.196) This equation is linear if all the coefficients (A, B, C, D, E, F and G) are a function of η and ζ , a constant, or zero. If they are a function of Φ or any of its derivatives, then the PDE is nonlinear. A partial differential equation is called quasi-linear if it is linear in the highest derivatives. In eq. (2.196), this means that A, B and C are a function of η and ζ , a constant, or zero. If a differential equation is quasilinear, it can be classified as an elliptic ( AC − B 2 > 0 ), parabolic ( AC − B 2 = 0 ) or hyperbolic ( AC − B 2 < 0 ) equation. The classification of a PDE describes how disturbances or changes propagate through a domain. If you are interested in a single point on the domain, you need to know two things: what region in the domain affects that point, and, if you do something at that point, what region in the domain it will affect. What region affects a point and what region that point affects are called the zone of dependence and the zone of influence, respectively. A representation of the zone of dependence and the zone of influence of elliptic, parabolic, and hyperbolic equations are displayed in Fig. 2.6. A multidimensional PDE that varies with time may be classified by looking at only the relation of two independent variables at a time. For instance look at how the PDE would be classified only in the x and y directions. Then examine the PDE in x and t. Finally examine the PDE in y and t. The zone of dependence is where the zone Chapter 2 Ge era i ed G ver i g Equati s 121 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Elliptic Parabolic Hyperbolic η η η ζ (a) ζ η ζ η η ζ η T2 T1 T2 (b) ζ η T2 ζ η T2 ζ (c) ζ Flow across a cylinder Isotherms in steady-state heat conduction Wake formed by a boat ζ Figure 2 6 Propagation of disturbances for different types of PDEs displayed as the (a) Zone of Dependence, (b) Zone of influence, and (c) a physical example of dependence for the three cases intersects. The zone of influence is where the zone of influence for the three cases intersects. An example of an elliptic equation is the Laplace equation. This equation governs physical problems, such as two-dimensional steady-state heat conduction, electrical potential, the free stream characteristics in boundary layer equations, and the pressure field in a porous medium with Darcy’s assumption. The independent variables, η and ζ , are spatial coordinates in these problems, and Φ is the potential field, representing temperature, voltage, the stream function, or pressure, respectively, for the given cases. 122 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Some examples of parabolic problems are unsteady heat conduction, boundary layer equations for momentum, energy and species as well as small vibrations in an elastic beam. For the unsteady heat equation, η is a spatial coordinate, and ζ is time. For the boundary layer equations, η is a spatial coordinate normal to the free stream velocity, and ζ is a spatial coordinate parallel to the free stream. A hyperbolic function is often called the wave equation. The most common examples of the wave equation are the propagation of shock waves and sound waves. The interface between phases is also hyperbolic. Another less intuitive example of the use of a hyperbolic function is heat conduction in micro- and nanoscale. It becomes relevant in these situations because heat actually conducts at a finite rate, and the limits of these rates are relevant at such small scales. For any problem to be well defined, there are boundary/initial conditions that must be applied. There are three basic boundary conditions for second order PDEs. These boundary conditions are the Dirichlet [ Φ = f (η , ζ ) ], the Neumann [ ∂Φ / ∂η = f (ζ ) , or ∂Φ / ∂ζ = f (η ) ], and the mixed [ a (∂Φ / ∂η ) + bΦ = f (ζ ) , or a(∂Φ / ∂ζ ) + bΦ = f (η ) ] type. The number of boundary conditions that must be applied for each independent variable is equal to the highest order of that variable. For example, if the PDE is second order in both η and ζ , then two boundary conditions are needed for each independent variable for a total of four boundary conditions. If ζ is time and first order, then an initial Dirichlet-type boundary condition is needed, and if it is second order, then both a Dirichlet and Neumann boundary condition are needed. There are implications to experimental measurements and numerical analysis based on the classification of the governing PDE. For example, experimental measurements in an incompressible flow field with moderate to low Reynolds numbers are very difficult, because the governing equations are highly elliptic. The elliptic nature means that disturbances downstream greatly affect the upstream flow field. Therefore, any measurements that require a device in the flow field may change the nature of that flow field, and are inherently inaccurate. Numerical simulations are very reliable for these cases when the flow is in the laminar regime. In the laminar flow regime, the full Navier-Stokes equations can be directly solved as an elliptic problem with no approximations. In the compressible flow regime with high Reynolds numbers, the characteristic of the flow is parabolic or hyperbolic, depending on the Mach number. In these cases, disturbances produced downstream do not affect the upstream flow field; therefore measurement devices in the flow will give an accurate depiction of the flow field without the device. In numerical simulations, it is computationally efficient to solve a parabolic flow field, because disturbances propagate in one direction. However, a flow field that is truly parabolic usually has a high Reynolds number and therefore is turbulent, which means it is three-dimensional and transient. Turbulence modeling either involves Chapter 2 Ge era i ed G ver i g Equati s 123 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press a very fine mesh, which is very computationally expensive, or an averaging technique, which requires a closure problem. The classification of a particular problem will help reduce the time it takes to get reliable results. 2 3 7 u p a d B u dary C diti s at the I terfaces The conservation equations introduced above can be applied within each phase and up to an interface. However, they are not valid across the interface, where sharp changes in various properties occur. Appropriate boundary conditions at the interface must be specified in order to solve the governing equations for heat, mass, and momentum transfer in the two adjoining phases. The interface conditions will serve as boundary conditions for the transport equations in the adjacent phases. Jump conditions at the interface can be obtained by applying the basic laws (conservation of mass, momentum, energy, and the second law of thermodynamics) at the interface. It is the objective of this subsection to specify mass, momentum, and energy balance at a non-flat liquid-vapor interface (see Fig. 2.7), as well as species balance in solid-liquid-vapor interfaces. For solidliquid or solid-vapor interfaces, these jump conditions can be significantly simplified. Mass Balance At a liquid-vapor interface, the mass balance is  ′′ mδ = ρ (V − VI ) ⋅ n = ρ v (Vv − VI ) ⋅ n (2.197)  δ is mass flux at the interface due to phase change and VI is the ′′ where m velocity of the interface. For a three-dimensional interface, there are three . g y Figure 2 7 Shape of the liquid-vapor interface near a vertical wall. 124 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press components of velocity: the normal direction, and two tangential directions, denoted by n, t1, and t2, respectively. Therefore, the velocity components should be defined according to these directions, as follows: V ⋅ n = Vn (2.198) V ⋅ t1 = Vt1 (2.199) V ⋅ t 2 = Vt 2 (2.200) (2.201) The interfacial mass balance can be rewritten in these terms by:  ′′ mδ = ρ (V ,n − VI ,n ) = ρ v (Vv,n − VI ,n ) Momentum Balance For applications involving thin film evaporation and condensation, the effects of surface tension and disjoining pressure will create additional forces on the interface. The momentum balance at the interface is ( τ ′ − τ ′v ) ⋅ n + σ (T )  1 1 +  n − pd n  RI RII   dσ   ′′ −  (∇Tδ ⋅ t )t = mδ (V − Vv )  dT  (2.202) On the left-hand side, the first term is the stress tensor, the second term is the capillary pressure, the third term is the disjoining pressure, and the fourth term is the Marangoni stress (Faghri and Zhang 2006). The right-hand side is the momentum transfer due to inertia. In this equation, the tangential direction, t, can either be t1 or t2. The stress tensor is: τ ′ = − pI + 2 μ D − μ ( ∇ ⋅ V ) I 2 3 (2.203) The deformation tensor can be written for a reference frame that is adjusted to the interface:  ∂Vn  ∂xn   ∂Vt1  1  ∂V 1 T D =  ∇ V + ( ∇V )  =   n +   2  2  ∂xt1 ∂xn   1  ∂Vn ∂Vt 2  2  ∂x + ∂x  n   t2 1  ∂Vn ∂Vt1 +  2  ∂xt1 ∂xn          ∂Vt1 ∂xt1 1  ∂Vt1 ∂Vt 2 +  2  ∂xt 2 ∂xt1                ∂Vt1 ∂Vt 2   1 +     2  ∂xt 2 ∂xt1    ∂Vt 2   ∂xt 2  1  ∂Vn ∂Vt 2 +  2  ∂xt 2 ∂xn  (2.204) The normal direction of the interface is [1 0 0], the first tangential direction is [0 1 0] and the second tangential direction is [0 0 1]. Therefore, Chapter 2 Ge era i ed G ver i g Equati s 125 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press τ ′ ⋅ n = [ − p 0 0] +  ∂V 2μ  n  ∂xn  1  ∂Vn ∂Vt1 +  2  ∂xt1 ∂xn      1  ∂Vn ∂Vt 2 +  2  ∂xt 2 ∂xn   2   − μ [∇ ⋅ V 0 0]  3      (2.205) This can be reduced to the three components to obtain: τ ′ ⋅ n ⋅ n = − p + 2μ ∂Vn 2 4 ∂V 2  ∂Vt1 ∂Vt 2 − μ∇ ⋅ V = − p + μ n − μ  + 3 ∂xn 3  ∂xt1 ∂xt 2 ∂xn 3   ∂Vn ∂Vt1 +  ∂xt  1 ∂xn         (2.206) τ ′ ⋅ n ⋅ t1 = μ  τ ′ ⋅ n ⋅ t2 = μ  (2.207) (2.208)  ∂Vn ∂Vt 2 +  ∂xt  2 ∂xn The momentum equation balance at the interface is then broken into its three components, as follows: Normal Direction 4  ∂V,n − p + pv +  μ 3  ∂x ,n  2   ∂V,t1 ∂V ,t 2 −  μ  +  3   ∂xt1 ∂xt 2  − μv ∂Vv,n   ∂x v,n         ∂Vv,t1 ∂Vv,t 2 +  − μv    ∂xt ∂xt 2 1   (2.209) 1 1  ′′ +σ  +  − pd = mδ (V,n − Vv,n )  RI RII  Tangential 1 μ   ∂V ,n ∂V ,t1 +  ∂xt ∂xn 1    ∂Vv,n ∂Vv,t1 +  − μv    ∂xt ∂xn 1     dσ −   dT    ∂Tδ    ∂xt1    ′′  = mδ V ,t1 − Vv,t1   ( ) ) (2.210) Tangential 2 μ   ∂V ,n ∂V ,t 2 +  ∂xt ∂xn 2    ∂Vv,n ∂Vv,t 2 +  − μv    ∂xt ∂xn 2     dσ −   dT    ∂Tδ    ∂xt 2    ′′  = mδ V,t 2 − Vv,t 2   ( (2.211) The non-slip condition at the liquid-vapor interface requires that V,t1 = Vv,t1 and V,t 2 = Vv,t 2 . The momentum balance at the tangential directions becomes μ  μ   ∂V,n ∂V,t1 +  ∂xn  ∂xt1   ∂Vv,n ∂Vv,t1 +  = μv    ∂xn   ∂xt1   ∂Vv,n ∂Vv,t 2 +  = μv    ∂xt ∂xn 2     dσ +    dT   dσ +   dT    ∂Tδ     ∂xt1   ∂Tδ    ∂xt 2          (2.212) (2.213)  ∂V,n ∂V,t 2 +  ∂xt ∂xn 2  126 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  ′′ For most applications, the evaporation or condensation rate, mδ , is not very  ′′ high; therefore, it can be assumed that mδ (V − Vv ) ≈ 0 . If the liquid and vapor phases are further assumed to be inviscid ( τ  = τ v ≈ 0 ), the momentum equation at the interface can be reduced to 1 1 pv − p = σ (T )  +  − pd RI RII   (2.214) Energy Balance Neglecting the work done by surface tension and disjoining pressure, the energy balance at the interface is ( q′′ − q′′ ) ⋅ n − (n ⋅ τ ′ ) ⋅ V,rel + (n ⋅ τ′v ) ⋅ Vv,rel v  2  Vv,rel  ′′ = mδ  ev +  2    V2,rel  −  e +  2      (2.215) If the velocity of the reference frame is taken as the interfacial velocity VI, eq. (2.215) can be rewritten as ( kv ∇Tv − k ∇T ) ⋅ n − (n ⋅ τ ′ ) ⋅ (V − VI ) + (n ⋅ τ ′v ) ⋅ (Vv − VI )  (V − VI ) 2   (V − VI ) 2     ′′  = mδ  ev + v  − e +  2 2       (2.216) where the heat fluxes in the liquid and vapor phases have been determined by Fourier’s law of conduction. Equation (2.216) can be rewritten in terms of enthalpy ( h = e + p / ρ ), i.e., ( kv ∇Tv − k ∇T ) ⋅ n − (n ⋅ τ ′ ) ⋅ (V − VI ) + (n ⋅ τ ′v ) ⋅ (Vv − VI ) (2.217)   p p 1 2 1  ′′ = mδ hv −  v −   + Vv − Vv ⋅ VI − V2 + V ⋅ VI  2    ρ v ρ  2   where hv is the difference between the enthalpy of vapor and liquid at saturation, i.e., the latent heat of vaporization. The stress tensor in eq. (2.217) can be expressed as τ ′ = − pI + 2μ D − μ ( ∇ ⋅ V ) I = − pI + τ 2 3 (2.218)  ′′ Since the relative velocities at the interface satisfy (V − VI ) ⋅ n = mδ / ρ  ′′ and (Vv − VI ) ⋅ n = mδ / ρv , we have  p p   ′′ p ( V − VI ) ⋅ n − pv ( Vv − VI ) ⋅ n = mδ  −  v −    ρ v ρ      (2.219) Substituting eqs. (2.218) and (2.219) into eq. (2.217), the energy balance can be written as Chapter 2 Ge era i ed G ver i g Equati s 127 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ( kv∇Tv − k ∇T ) ⋅ n − ( n ⋅τ  ) ⋅ ( V − VI ) + ( n ⋅τ v ) ⋅ ( Vv − VI ) 12 1   ′′  = mδ  hv + Vv − Vv ⋅ VI − V2 + V ⋅ VI  2 2   (2.220) To simplify the energy equation, the kinetic energy terms are considered negligible and no-slip conditions are assumed at the interface, (2.221) V ,t = Vv,t = VI ,t Therefore, the energy equation can be rewritten as  ∂Tv ∂T  − k    kv ∂xn   ∂xn  4 ∂V,n 2  ∂V,t1 ∂V,t 2  ′′ = mδ  hv + ν  − ν  + ∂xn ∂xt 2 3 3  ∂xt1    4 ∂Vv,n 2  ∂Vv,t1 ∂Vv,t 2 − νv + νv  + ∂xn ∂xt 2 3 3  ∂xt1          (2.222) where ν  and ν v are kinematic viscosities of liquid and vapor phases, respectively. The energy balance at the interface can be simplified by assuming that the change in the kinetic energy across the interface is negligible, i.e.,  ′′ (2.223) ( kv∇Tv − k ∇T ) ⋅ n = mδ hv Equations (2.197), (2.214), and (2.223) are widely used in the analysis of evaporation and condensation. Species Balance For a general interface between phases k and j in a multi-component system, a local balance in mass flux of species i must be upheld. The total species mass  flux, mi′′ , at an interface is:  mi′′ = ρk ,i ( Vk ,i − VI ) ⋅ n = ρ j ,i V j ,i − VI ⋅ n ( ) (2.224) The velocity of species i in phase k and phase j is Vk ,i and V j ,i , respectively. These velocities are defined as: ρk ,i Vk ,i = J k ,i + ωk ,i ρk Vk ρ j ,i V j , i = J j ,i + ω j ,i ρ j V j (2.225) (2.226) Using the interfacial species mass balance, and substituting the definition of the species velocity, the interfacial species mass flux is:  mi′′ = J k ,i ⋅ n + ωk ,i ρk ( Vk − VI ) = J j ,i ⋅ n + ω j ,i ρ j ( V j ,k − VI ) (2.227) Remembering the overall mass conservation at the interface,  m ′′ = ρk ( Vk − VI ) ⋅ n = ρ j ( V j − VI ) ⋅ n The interfacial species mass flux is: (2.228) 128 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press    mi′′ = J k ,i ⋅ n + ωk ,i m ′′ = J j ,i ⋅ n + ω j ,i m ′′ (2.229) In some problems, the species mass flux will be specified, and the total mass flux is simply a sum of all the species mass fluxes.  m ′′ =   m′′ i i (2.230) If the species mass flux is not specified, the total mass flux at an interface can be calculated from the interfacial species mass flux equation, eq. (2.229):  m ′′ = ( J j ,i − J k ,i ) ⋅ n ωk ,i − ω j ,i (2.231) Since the species equation is second order in space, and there are two phases, there needs to be a secondary condition relating the species mass fraction in phases k and j. This can be done by expressing the mass fraction of species i in phase k as a function of the species mass fraction i in phase j, or vice versa. ωk ,i = ωk ,i (ω j ,i ) (2.232) ω j ,i = ω j ,i (ωk ,i ) (2.233) Note that, in general, the species mass fraction is not a continuous function, and there is almost always a jump condition at the interface when two or more species are present. To make this point clear, the simplest case of a binary mixture is considered. For a binary mixture, the species diffusion flux, J, can be calculated by Fick’s law. J k ,1 ⋅ n = − ρk Dk ,12∇ωk ,1 ⋅ n (2.234) (2.235) J j ,1 ⋅ n = − ρ j D j ,12∇ω j ,1 ⋅ n There are several examples of different phenomena of a binary mixture in the presence of a liquid/vapor, solid/vapor, or liquid/solid interface. In other examples, the species flux in one of the phases is unimportant; therefore, only one phase must be considered. A classical example of species flux is the condition of zero mass flux due to an impermeable surface where species A does not diffuse to the stationary media ( ∇x A ⋅ n = 0 ) or ( ∂x A / ∂y = 0 ) . The case of constant surface concentration is more typical; however, all types of mass transfer boundary conditions are demonstrated in three categories below. a. The mass transfer from solid or liquid to a gas stream Evaporation and sublimation are typical examples of mass transfer from liquid or solid to a gas mixture, as shown in Fig. 2.8, where x A, , x A,s , x A, g are ′ molar fractions in liquid, solid, and gas, respectively, and m ′A is mass flux for species A. There are two cases shown in Fig. 2.8. In case I, a pure liquid/solid is evaporating/sublimating into a gaseous mixture and the evaporation /sublimation rate is controlled purely by the species gradient in the gas phase. In case II, the liquid/solid is not pure, therefore concentration gradients also exist Chapter 2 Ge era i ed G ver i g Equati s 129 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1 I. x A, I. x A, s II. x A , xA II. x A , s A m′′ Gas I. x A, g Liquid Solid A m′′ Gas 0 y (a) Liquid to Gas II. x A, g y I. x A, g II. x A, g (b) Solid to Gas Figure 2 8 Species concentration and mass transfer from solid and or liquid to a gas mixture. in the liquid/solid phases. The mass transfer rate is limited by both phases. Also note that the concentration gradient in the solid [Fig. 2.8(b)] is steeper and does not penetrate as deeply when compared to the liquid concentration gradient [Fig. 2.8(a)]. Also, the liquid concentration gradient is steeper and does not penetrate as far compared to the gas. These trends are due to the ratio of the mass diffusivities in each of the phases. A simplified relationship for the mass concentration condition at the interface can be obtained assuming the gas mixture is approximated by the ideal gas and the solid or liquid has a high concentration of A. The condition within the gas stream is of interest, because the main resistance to mass transfer is within that region. With the preceding assumptions, the partial vapor pressure of A in the gas mixture at the interface can be approximated from Raoult’s Law: p A y = 0− = x A ,  (2.236) ( p A,sat ) Liquid to gas y = 0+ pA y = 0− = x A,s y = 0+ ( p A,sat ) Solid to gas (2.237) y = 0+ where pA is the partial vapor pressure of A in the gas mixture, x A, mole fraction of species A in liquid, x A, s y = 0+ is the is the mole fraction of species in solid, and pA,sat is the normal saturation pressure of species A at the surface interface. Clearly, if the solid or liquid is made of pure species A, then x A , = x A, s = 1 and eqs. (2.236) and (2.237) reduce to y =0 y =0 pA y = 0− = p A,sat (2.238) This means that the partial vapor pressure of species A in the gas mixture at the interface is equal to the normal saturation pressure for species A, which is a function of interfacial temperature and can be obtained from thermodynamic tables in Appendix E. At the interface, species A and B in the gas phase must be in equilibrium with species A and B in the liquid or solid phase, except in 130 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press extreme circumstances. Knowing p A y = 0− and total pressure p, one can easily calculate the mole fraction x A, g and mass fraction ω A, g of species A at the interface on the gas side by the following relation x A, g = pA y = 0− ω A, g = p x A, g M A (2.239) (2.240) x A, g M A + x B , g M B where MA and MB are molecular weights of species A and B. A conventional example for the case in Fig. 2.8(a) is the evaporation of water to a water-air mixture; obviously, one can imagine that the water absorbed a small amount of air. Sublimation of iodine in air is an application for the case of Fig. 2.8(b). Mass transfer from a solid to a gaseous state sometimes requires the specification of diffusion molar flux, rather than concentration, at the solid surface. J * = −cD AB A ∂x A ∂y (2.241) where J * can be a function of concentration and not necessarily constant. A  One application is related to catalytic surface reaction. For this case, m ′′ = 0 because the mass flux of reactants equals the mass flux of products. Therefore,  from eq. (2.229), mi′′ = J k ,i ⋅ n . Catalytic surfaces are used to promote heterogeneous reactions, which occur at the surface; the appropriate boundary condition is ′ J * = −k1 c A A where k1′ is the reaction rate constant and n is the order of reaction. b. Mass transfer from gases to liquids or solids There are two major forms of mass transfer from gas to liquid, shown in Fig. 2.9. Case I is for condensable species in the liquid phase, or for species that can deposit vapor into the solid phase. In this case, a species of low concentration in the gas phase can condense at a higher concentration in the liquid phase. This phenomenon is important in distillation processes. Case II is for species that are weakly soluble in a liquid or a solid. In the liquid ( x A, is small), Henry’s law will relate the mole fraction of species A in the liquid to the partial vapor pressure of A in the gas mixture at the interface by following relation [Fig. 2.9(a)] x A, y =0 + ( y =0 ) n (2.242) = pA y = 0− H (2.243) Henry’s constant, H, is dependent primarily on the temperature of the aqueous solution and the absorption of two species. The pressure dependence is usually small and in fact, negligible up to pressure equal to +5 bar. Table E.7 Chapter 2 Ge era i ed G ver i g Equati s 131 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1 Liquid I. x A, Gas II. x A, g I. x A, g I. x A, s Gas II. x A, g xA II. x A, I. x A, g Solid II. x A, s 0 a. Gas to Liquid y b. Gas to Solid Figure 2 9 Species concentration and mass transfer from gas mixture to liquid or solid provides H for selected aqueous solutions. For binary soluble gases, Henry’s law is not appropriate, and solubility data are usually presented in terms of gas-phase partial pressure to liquid-phase mole fraction. Such data are given in Appendix B for NH3-water and SO2-water systems. In chemical engineering applications, there are many cases in which gas is absorbed into a liquid; two examples are the absorption of hydrogen sulfide from an H2S-air mixture into liquid water and the absorption of O2 or chlorine into liquid water. The concentration of gas in a solid at the interface is usually obtained by the use of a property known as the solubility, S, defined below c A, s = S p A, g (2.244) − y = 0+ y =0 where c A,s is the mole concentration of species A in the solid at the interface and p A, g y = 0− is the vapor partial pressure of species A in the gas at the interface. The values of S ( kmol/Pa-m3 ) for vapor gas-solid combinations are given in chemistry handbooks; some are in Appendix E (Tables E.6, E.8 and E.9). Special care should be exercised when using eq. (2.244) for the variation in form and units – two different versions of S are presented in Appendix E. An example of this case is the diffusion of helium in a glass. The dissolution of gas in metals is much more complex and depends on the type of metal used. The dissolution of gas in some metals can be reversed, such as the case of hydrogen into titanium. In contrast to hydrogen, oxygen dissolution in titanium is irreversible but is complicated because it forms a layer of scale TiO2 on the surface. This topic is beyond the scope of this book, and interested readers should refer to a metallurgy or materials handbook to consult appropriate phase diagrams. 132 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 1 I. x A, s Liquid Solution Solid I. x A, II. x A, Liquid Solution xA II. x A, s Solid II. x A, s II. x A, I. x A, I. x A, s (b) Liquid to Solid 0 (a) Solid to Liquid Figure 2 10 Species concentration and mass transfer from solid to liquid and liquid to solid c. Mass transit from solid to liquid or liquid to solid Two cases are also presented for the case of mass transfer between a solid and a liquid (Fig. 2.10). Case I represents a pure solid dissolving into a liquid, or a pure liquid diffusing into a solid. Case II represents a solid mixture melting into a liquid mixture, or a liquid mixture solidifying into a solid mixture. For the first case, the mass transfer is related to the solubility, which can be found in chemistry handbooks (Lide, 2009) and some are reproduced in Appendix E. A conventional example for this case is the dissolution of salt into water. For a liquid/solid interface during melting and solidification, the ratio of the species concentration of the solid in liquid phases is called the partition ratio, Kp. Kp = c A, s c A, (2.245) The partition coefficient is a function of temperature. Furthermore, when the liquid and solid lines are nearly straight, it is a constant. This information can be obtained from phase diagrams for different solid/liquid mixtures. Example 2.7 A container with an open top is exposed to superheated water vapor at temperature Tv, with the bottom and the side walls insulated. At time t=0, the bottom surface of the vessel comes into contact with the air at temperature T∞, which is below the vapor temperature Tv. The sides and top of the vessel remain insulated for t > 0 [see Fig. 2.11(a) and (b)]. Formulate the governing equations for the condensation problem. If the vapor is saturated (Tv = Tsat), find the instantaneous thickness of the condensate. Solution: This is a one-dimensional condensation problem because the sides of the vessel are insulated. The thickness of the bottom wall is very small and therefore the temperature drop across the bottom wall can be neglected. In addition, it is assumed that the thermal properties of both Chapter 2 Ge era i ed G ver i g Equati s 133 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press liquid and vapor are constants. The interfacial velocity is related to the instantaneous liquid thickness by dδ wI = (2.246) dt The continuity equation at the interface can be obtained by using eq. (2.197) in the z-direction, i.e., ρ (w − wI ) = ρ v (wv − wI ) The condensate is stationary, so w = 0 . Therefore, the vapor phase velocity is ρ  dδ wv =   − 1 (2.247)  ρv  dt For the liquid phase, the governing energy equation and the corresponding initial and boundary conditions are ∂T ∂ 2T = α ∂t ∂z 2 T ( z , 0) = Tsat (2.248) (2.249) (a) Insulated (t < 0) . (b) Bottom cooled ( t > 0 ) Figure 2.11 Thin-walled vessel containing two-phase mixture. 134 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press k ∂T = h[T (0, t ) − T∞ ] ∂z (2.250) For the vapor phase, the governing energy equation and the corresponding initial and boundary conditions are  dδ ∂Tv ∂Tv  ρ ∂ 2Tv + − 1 = αv (2.251) ∂t  ρ v ∂z 2  dt ∂z Tv ( z , 0) = Tv (2.252) Tv (∞, t ) = Tv (2.253) The liquid and vapor temperatures at the interface must be the same, i.e., Tv (δ , t ) = T (δ , t ) = Tsat (2.254) The energy balance at the interface can be obtained by using eq. (2.220). Since the velocity in the liquid is zero and the velocity in the vapor phase is not very high, the kinetic energy and the viscous dissipation in eq. (2.220) can be neglected. The resulting energy balance equation at the interface is ∂T (δ , t ) ∂T (δ , t ) dδ −kv v + k  = ρ hv (2.255) ∂z ∂z dt The analytical solution of this problem will be very difficult. If the vapor is saturated, the governing equations for the vapor phase will not be necessary because the vapor temperature will be uniformly equal to the saturation temperature. This is a one-domain problem because only the liquid domain needs to be investigated. The problem can be further simplified by introducing the quasi-steady state assumption. This step eliminates the transient term in the energy equation for the liquid phase because condensation is a very slow process and the change of internal energy in the liquid can be neglected. The energy equation in the liquid phase is thus reduced to 0 = α ∂ 2T ∂z 2 (2.256) which is subject to the boundary conditions specified by eqs. (2.250) and (2.254). The energy balance at the interface is reduced to ∂T (δ , t ) dδ = ρ hv (2.257) k  The temperature in the liquid film can be obtained by integrating eq. (2.256) twice: T ( z ) = c1 z + c2 (2.258) where c1 and c2 are integral constants. They can be determined using eqs. (2.249) and (2.250), i.e., k c1 = h ( c2 − T∞ ) (2.259) Tsat = c1δ + c2 (2.260) which can be solved to yield ∂z dt Chapter 2 Ge era i ed G ver i g Equati s 135 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Figure 2.12 Temperature distribution for a two-phase mixture in a thin-walled vessel. c1 = h (Tsat − T∞ ) k + hδ (t ) c2 = Tsat − h(Tsat − T∞ ) δ (t ) k + hδ (t ) (2.261) Substituting eq. (2.261) into eq. (2.258) yields the temperature distribution in the liquid film T ( z , t ) = Tsat + h (Tsat − T∞ ) [ z − δ (t )] k + hδ (t ) (2.262) The temperature distribution in the vessel can be illustrated by Fig. 2.12. At this point, eq. (2.262) can be substituted into eq. (2.257). An ordinary differential equation for the liquid thickness is obtained as a result: dδ k h (Tsat − T∞ ) = ρ hv (2.263) dt k + hδ (t ) which is subject to the following initial condition: δ (0) = 0 (2.264) Rearranging eq. (2.263) so that δ appears only on the left-hand side, and t appears only on the right-hand side, gives h  h (Tsat − T∞ ) dt  δ + 1 dδ = ρ hv  k  (2.265) Integrating eq. (2.265) with the boundary condition, eq. (2.264) yields h (Tsat − T∞ ) h2 t =0 δ +δ − 2k ρ hv (2.266) which is a quadratic equation in the form of ax 2 + bx + c = 0 . This allows for one positive and one negative solution of δ , and the latter does not physically make sense. Therefore, the liquid film thickness is  2h 2 (Tsat − T∞ )t  2 hδ (t ) = −1 + 1 +  k k ρ hv     1 (2.267) 136 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2 3 8 Rarefied Vap r Se f Diffusi de The discussions so far are limited to the case where the density of the fluid is sufficiently high to permit the continuum assumption. There are cases in which the assumption of continuum is not valid ( Kn ≥ 0.001 , see Section 1.4.1). For example, in the early stage of heat pipe startup from the frozen state, the vapor pressure and density are very small in the heat pipe core. Because of the low density, the vapor in the rarefied state is somewhat different from the conventional continuum state. Also, the vapor density gradient is very large along the axial direction of the heat pipe. The vapor flow along the axial direction is mainly caused by the density gradient via vapor molecular diffusion. The low-density vapor state that has partly lost its continuum characteristics is referred to as rarefied vapor. Neglecting the presence of noncondensable gases, the rarefied vapor flow can be simulated by a self-diffusion model. The term self-diffusion here means the interdiffusion of particles of the same mass due to a gradient in density. The governing equations for startup of a heat pipe from the frozen state were derived by applying the principles of the conservation of mass and energy in a differential cylindrical control volume in conjunction with the definition of mass flux (Cao and Faghri, 1993). The mass self-diffusion equation is ∂ρ ∂  ∂ρ  1 ∂  ∂ρ  −  Dv (2.268) −  rDv =0 ∂t ∂z  ∂z  r ∂z  ∂r  and the energy equation is ∂ ( ρ cvT ) 1 ∂ ∂t + r ∂r  ′′ (mr c p rT ) + 1∂ ∂ ∂T  ∂  ∂T   ′′ (m z c p T ) =  rkv  +  Dv  r ∂z  ∂z ∂r  ∂z  ∂z  (2.269) ′ ′ where the mass fluxes mr′ and mz′ are ∂ρ ∂r ∂ρ  ′′ mz = ρ w = − Dv ∂z  ′′ mr = ρ v = − Dv (2.270) (2.271) where Dv is the self-diffusion coefficient, and kv is the vapor molecular conductivity. The evaluation of low-density properties such as Dv is carried out using the kinetic theory of gases. The coefficient of self-diffusion is obtained from the relation based on the Chapman-Enskog kinetic theory (Hirschfelder et al., 1966): Dv = 2.628 × 10−7 T 3 / MW pσ 2 Ω D (kbT / ε ) (2.272) where p is the pressure in atmospheres, σ is the collision diameter in Å, ε is the maximum energy of attraction between a pair of molecules, kb is the Boltzmann constant, and Ω D is the collision integral for mass diffusion. For Chapter 2 Ge era i ed G ver i g Equati s 137 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press sodium, both σ and ε / kb can be found from the table for constants of the Lennard-Jones potential model (Edwards et al., 1979), which gives σ = 3.567 Å, and ε / kb = 1375 K. The value of the collision integral Ω D can also be found from the same reference, which is listed as a function of ε / kb . 2 3 9 A Exte si C busti Combustion is an exothermic chemical reaction process between fuel and oxidant. If combustion involves a liquid fuel, the liquid fuel does not actually burn as a liquid; it is vaporized first and diffuses away from the liquid-vapor surface. Meanwhile, the gaseous oxidant diffuses toward the liquid-vapor interface. Under the right conditions, the mass fluxes of vapor fuel and gaseous oxidant meet and the chemical reaction occurs at a certain location known as the flame (Lock, 1994; Avedisian, 1997, 2000). The flame is usually a very thin region with a color dictated by the temperature of the combustion. The temperature and mass concentration distributions during a combustion process can be represented as shown in Fig. 2.13. The initial temperature of a liquid fuel is T0, and the temperature of the liquid–vapor interface, TI, is at the dew point temperature of the fuel. The temperature reaches a maximum at the location of the flame, and decreases with increasing x until it reaches T∞. The mass fraction of fuel, ωf, is maximal at the liquid-vapor interface and decreases as location of the flame is approached. The mass fraction of the oxidant, ωo, on T, ω Liquid-vapor interface T∞ ωo TI ωf ωp T0 0 sI sf x Figure 2.13 Combustion near a planar surface. 138 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the other hand, is maximal at infinity and decreases as location of the flame is approached. The mass fractions of both vapor fuel and oxidant reach their minimum values at the location of the flame. The mass fractions of the products of combustion, ωp, are at their maximum at the location of the flame and decrease as x either increases or decreases. The transport phenomena involved in combustion include heat and mass transfer in both the liquid fuel and the gaseous mixture. The governing equations for the gaseous phase will be discussed below. Firstly, the gaseous mixture must satisfy the continuity equation, i.e., Dρ + ρ∇ ⋅ V = 0 (2.273) Dt and the gas flow is governed by the momentum equation, ρ DV = ∇ ⋅ τ + ρg Dt (2.274) where the shear stress can be obtained from eq. (2.70) and gravity is the only body force considered. The energy equation for combustion is ρ Dh = ∇ ⋅ k ∇T Dt (2.275) which is a simplified version of eq. (2.92), with internal heat generation and viscous dissipation neglected. The effect of the substantial derivative of pressure on the energy is also neglected. The specific enthalpy for the mixture is related to the specific enthalpy of each component in the mixture by h= ω h i =1 N ii (2.276) It is a common practice in thermodynamic and heat transfer analyses to consider the change in enthalpy during a chemical reaction process rather than the absolute values of enthalpy. For a process that does not involve chemical reaction, we can choose any reference state for an individual substance and define the enthalpy at that reference state as zero. For example, the reference state for water is often chosen at the triple point, while the reference state for an ideal gas is often chosen as zero K. However, when chemical reaction is involved in a process, as is the case with combustion, the composition of the system changes during the chemical reaction; therefore, the reference states for all reactants and products must be the same. One convenient option in such a situation is to segregate the enthalpy of any component into two parts: (1) the enthalpy due to its chemical composition at the standard reference state (at 25 ˚C and 1 atm), and (2) the sensible enthalpy due to any temperature deviation from the standard reference state. Therefore, the enthalpy for the ith component in the mixture can be expressed as hi = hi + hiT (2.277) Chapter 2 Ge era i ed G ver i g Equati s 139 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Tab e 2 1 Enthalpy of formation ho for selected substances at 25 ˚C and 1 atm Substance Acetylene Benzene Carbon Carbon monoxide Carbon dioxide Ethyl alcohol Ethyl alcohol Ethylene Ethane Hydrogen Methane Nitrogen Nitrogen n-Dodecane n-Octane n-Octane Oxygen Propane Water Water vapor Formula (phase) C2H2 (g) C6H6 (g) C (s) CO (g) CO2 (g) C2H5OH (g) C2H5OH (  ) C2H4 (g) C2H6 (g) H2 (g) CH4 (g) N2 (g) N (g) C12H26 (  ) C8H18 (g) C8H18 (  ) O2 (g) C3H8 (g) H2O (  ) H2O (g) Enthalpy of formation (kJ/kmol) 226,730 82,930 0 -110,530 -393,520 -235,310 -277,690 52,280 -84,680 0 -74,850 0 472,650 -291,010 -208,450 -249,950 0 -103,850 -241,820 -285,830 where hi is the enthalpy of formation for the ith component, i.e., the enthalpy due to its chemical composition at the standard reference state. The enthalpy of formation for selected substances is shown in Table 2.1. The sensible enthalpy, hiT , is related to temperature by hiT =  T T c pi dT (2.278) Substituting eqs. (2.277) and (2.278) into eq. (2.276), the enthalpy of the mixture becomes h= ω h +   ii N T i =1 T c p dT (2.279) where cp is the average specific heat of the mixture, defined as cp = ω c i =1 N i pi (2.280) Substituting eq. (2.279) into eq. (2.275), the energy equation for combustion becomes N D (c p T ) Dωi = ∇ ⋅ k∇T − ρ hi (2.281) ρ Dt  i =1 Dt  If the fuel is consumed at a rate of m ′f′′ per unit volume, hc – the heat of combustion – is defined by 140 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press hc = − h  m′′′  f i =1 ρ N i Dωi Dt (2.282) Substituting eq. (2.277) into eq. (2.282) gives N N ρ ρ  Dωi T Dωi hc = − f m ′′′ h i =1 i Dt − f m ′′′ h i =1 i Dt (2.283) The contribution of the second term on the right-hand side of eq. (2.283) is negligible, since hc  hiT . Dropping the second term from the right-hand side of eq. (2.283) and substituting the result into eq. (2.276), the energy equation becomes ρcp DT f = ∇ ⋅ k∇T + m′′′hc Dt (2.284) The mass fraction of each component (fuel, oxidant, and product) is dominated by Dωi  ρ (2.285) = −∇ ⋅ J i + mi′′′ Dt where the subscript i can be f (fuel), o (oxidant), or p (product). The ratio of the rates of oxygen and fuel consumption is defined as the oxygen/fuel ratio: γ=  ′′′ mo f m ′′′ (2.286) The above analysis applies to combustion occurring on a planar surface. For many applications, combustion of the liquid fuel is usually preceded by breaking up a fuel jet into liquid droplets so that combustion occurs around a spherical liquid droplet. Combustion of a falling liquid droplet, as shown in Fig. 2.14, is analyzed here. The liquid fuel droplet vaporizes at the dew point, Td, which is the saturation temperature corresponding to the partial pressure of the fuel vapor in T, ω Tf Td T∞ ωf ωp ω0 rI rf r Figure 2 14 Combustion near a liquid fuel droplet. Chapter 2 Ge era i ed G ver i g Equati s 141 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press the mixture near the liquid-vapor interface. To simplify the analysis, it is assumed that the temperature in the liquid fuel droplet is uniformly equal to the saturation temperature of the fuel at the total system pressure, i.e., the mass fraction of the fuel at the liquid-vapor interface equals one. It is further assumed that the shapes of both the liquid fuel droplet and the flame are spherical, which allows for application of a one-dimensional symmetric model (Lock, 1994), which is presented here. If the combustion process is assumed to be in a quasisteady state (neglecting the transient term in the governing equation), the energy equation for combustion becomes 1d 2 1 d  2 dT (r ρ c p uT ) = 2 r k 2 dr dr r r dr    ′′′  + m f hc  (2.287) where u is the velocity of the gaseous mixture in the radial direction. The mass fraction of the oxidant can be obtained by simplifying eq. (2.285), i.e., dωo  1d 2 1 d2  ′′′ (r ρ uωo ) = 2 (2.288) r ρ Do  − mo 2 r dr r dr  dr  where Do is the mass diffusivity of the oxidant in the gaseous mixture. Multiplying eq. (2.288) by hc / γ and combining the resulting equation with eq. (2.287) yields  h 1 d2 r ρ c p u  T + ωo c 2  γ cp r dr      1 d  2 dT d  hc   rk + r 2 ρ Do  =  ωo   (2.289)   r 2 dr  dr dr  γ    where eq. (2.286) was used to eliminate the rates of fuel and oxygen consumption. Assuming that the Lewis number [ Le = k / ( ρ c p Do ) ] equals one and the specific heat is constant, eq. (2.289) can be simplified as follows: 1d 2 1 d  2 dT *  r ρ c p uT * = 2 r k  dr  r 2 dr r dr    ( ) (2.290) where T * = T + ωo hc γ cp (2.291) is a modified temperature in the combustion process. Equation (2.290) can be rearranged to the following form: d  2 dT *   * r  ρ c p uT − k  = 0 dr   dr     (2.292) The continuity equation requires that  4π r 2 ρ u = 4π rI2 ( ρ u) I = m (2.293) Since it has been assumed that the mass fraction of the fuel in the mixture at  the liquid-vapor interface equals one, the mass flow rate m reflects the mass 142 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press flow rate of the fuel. The mass flux is often used in combustion analysis, and it is defined as  m′′ =  m Equations (2.292) and (2.293) can be rewritten in terms of mass flux, i.e., d  2 dT *   *  r  m ′′c pT − k  = 0 dr   dr       ′′  4π r 2 m′′ = 4π rI2 mI = m 4π r 2 (2.294) (2.295) (2.296) Integrating eq. (2.295) from the liquid fuel droplet surface ( r = rI ) to an arbitrary radius ( r > rI ) yields  dT *  dT *  2 *   r 2  m ′′c pT * − k  = rI  m′′c pT − k    dr  dr  I     (2.297) Substituting eq. (2.296) into eq. (2.297), one obtains r 2k  dT *  dT *  ′′ = rI2 mI c p (T * − TI* ) + rI2  k  dr   dr  r = r (2.298) I ′ The mass flow rate at the surface of the liquid fuel droplet, mI′ , is the same  as the fuel burning rate, m′f′ , because the mass fraction at the surface of the droplet is one. Introducing the excess modified temperature, θ = T * − TI* eq. (2.298) becomes r 2k  dθ 1 f = rI2 m′′ c p θ + f dr m′′ c p    dθ  k   dr r =rI     (2.299) (2.300) Introducing a new dependent variable, 1  dθ  ϕ =θ + k  eq. (2.300) becomes 1 f m′′ c p  dr r =rI f rI2 m′′ c p 1 dr k r2 2 (2.301) ϕ Integrating eq. (2.302) over the interval of (r , ∞) , one obtains ln dϕ = (2.302) f ϕ∞ rI m′′ c p 1 = (2.303) ϕ k r The fuel burning rate, Gf, in eq. (2.303) is related to the heat transfer at the liquid droplet surface as follows  dθ  f (2.304) m ′′ hv = qI ′′ = k    dr r =rI Substituting eq. (2.304) into eq. (2.301), the new dependent variable ϕ becomes Chapter 2 Ge era i ed G ver i g Equati s 143 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ϕ =θ + hv cp (2.305) Substituting eq. (2.305) into eq. (2.303) and letting r = rI , an equation for the fuel burning rate is obtained:  c pθ ∞  k f m ′′ = ln  1 + (2.306)  rI c p  hv  Equation (2.306) can also be used to determine the transient liquid fuel droplet size, because the fuel burning rate is related to the size of the droplet by d ( ρ f 4π rI3 / 3) / dt dr f m′′ = − = −ρ f I (2.307) 2 dt 4π rI where ρ f is the density of the liquid fuel. Substituting eq. (2.307) into eq. (2.306), one obtains  c pθ ∞  dr k rI I = − ln  1 + (2.308)  dt hv  ρ f cp  Integrating eq. (2.307) and considering the initial condition rI = ri at t = 0, the liquid droplet radius becomes  c pθ ∞  2kt ln 1 + rI2 = ri2 − (2.309)  hv  ρ f cp  Equation (2.309) can be used to estimate the time needed to completely burn the liquid fuel droplet: ρ f c p ri2 (2.310) tf =  c pθ ∞  2k ln  1 +  hv   Example 2.8 Liquid hydrocarbon fuel is suspended in air in the form of small droplet, and the oxygen/fuel ratio is γ = 3.5 . Before the liquid fuel is sprayed into the air, the mass fraction of oxygen is ωo,∞ = 0.21 , and the temperature of the air is T∞ = 25 o C . The latent heat of evaporation for the fuel is hv = 360 kJ/kg , and the enthalpy of combustion is hc = 4.4 × 104 kJ/kg . The specific heat and thermal conductivity of the gaseous mixture are c p = 1.15 kJ/kg-K and k = 0.07 W/m-K, respectively. The initial temperature of the liquid droplet can be assumed to be equal to the saturation temperature of the fuel, 187 °C. The density of the liquid fuel is ρ f = 700 kg/m3 . Estimate the time needed to completely burn a liquid fuel droplet with a diameter of 40 μm . 144 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Solution: The excess modified temperature at the locations far away from the fuel droplet can be obtained by substituting eq. (2.291) into eq. (2.299): * θ∞ = T∞ − TI* = (T∞ − TI ) + (ωo,∞ − ωo, I ) hc γ cp = (25 − 187) + (0.21 − 0) × The time required to completely burn the liquid fuel droplet is then obtained by using eq. (2.310), i.e., ρ f c p ri2 700 × 1.15 × 103 × (20 × 10−6 ) 2 = = 0.00112s tf =  c pθ ∞   1.15 × 2133.65  2k ln 1 +  hv   2 × 0.07 × ln 1 +  360 4.4 × 104 = 2133.65 o C 3.5 × 1.15   2 4 V u e Averaged de s 241 verview f Averagi g Appr aches The objectives of the various averaging methods are twofold: (1) to define the average properties for the multiphase system and correlate the experimental data, and (2) to obtain solvable governing equations that can be used to predict the macroscopic properties of the multiphase system. This chapter will address the application of averaging methods to the governing equations. Based on the physical concepts used to formulate multiphase transport phenomena, the averaging methods can be classified into three major groups: (1) Eulerian averaging, (2) Lagrangian averaging, and (3) Molecular statistical averaging. These averaging techniques are briefly reviewed below. 2 4 1 1 Eu eria Averagi g Eulerian averaging is the most important and widely-used method of averaging, because it is consistent with the control volume analysis that we used to develop the governing equations in the preceding section. It is also applicable to the most common techniques of experimental observations. Eulerian averaging is based on time-space description of physical phenomena. In the Eulerian description, changes in the various dependent variables, such as velocity, temperature, and pressure, are expressed as functions of time and space coordinates, which are considered to be independent variables. One can average these independent variables over both space and time. The integral operations associated with these Chapter 2 Ge era i ed G ver i g Equati s 145 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press averages smooth out the local spatial or instant variations of the properties within the domain of integration. For a generalized function Φ = Φ ( x , y, z , t ) , the most widely-used Eulerian averaging includes time averaging and volumetric averaging. The Eulerian time average is obtained by averaging the flow properties over a certain period of time, Δt, at a fixed point in the reference frame, i.e., Φ= 1 Δt  Δt Φ ( x , y, z , t )dt (2.311) for this equation, the time period Δt is chosen so that it is larger than the largest time scale of the local properties’ fluctuation, yet small enough in comparison to the process macroscopic time scale. During this time period, different phases can flow through the fixed point. Eulerian time averaging is particularly useful for a turbulent multiphase flow as well as for the dispersed phase systems (Faghri and Zhang, 2006). Eulerian volumetric averaging is usually performed over a volume element, ΔV, around a point ( x , y, z ) in the flow. For a multiphase system that includes Π different phases, the total volume equals the summation of the individual phase volumes, i.e., ΔV =  ΔV k =1 Π k (2.312) The volume fraction of the kth phase, ε k , is defined as the ratio of the elemental volume of the kth phase to the total elemental volume for all phases, i.e., εk = ΔVk ΔV =1 (2.313) The volume fraction of all phases must sum to unity: ε k =1 Π k (2.314) Eulerian volume averaging is expressed as Φ= 1 ΔV  k =1 Π ΔVk Φ k ( x , y, z , t )dV (2.315) where the volume element ΔV must be much smaller than the total volume of the multiphase system so that the average can provide a local value of Φ in the flow field. The volume element ΔV must also be large enough to yield a stationary average. Since the volume element includes different phases, information about the spatial variation of Φ for each individual phase is lost and Φ represents the average for all phases. For any variable or property that is associated with a particular phase, Φ k , the phase-average value of any variable or property for that phase is obtained with the following equations 146 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Intrinsic phase average: Φk k = 1 ΔVk 1 ΔV  ΔVk Φ k dV (2.316) Extrinsic phase average: Φk =  ΔVk Φ k dV (2.317) Intrinsic means that it forms to the inherent part of a phase and is independent of other phases in the volume element. In contrast, extrinsic means it is a property that depends on the phase’s relationship with other phases in the volume element. While the intrinsic phase average is taken over only the volume of the kth phase in eq. (2.316), the extrinsic phase average for a particular phase is taken over an entire elemental volume in eq. (2.317). These two phase-averages are related by k Φk = ε k Φk (2.318) The intrinsic and extrinsic phase averages defined in eqs. (2.316) and (2.317) are related to the volume average defined in eq. (2.315) by Φ=  k =1 Π Φk = ε k =1 Π k Φk k (2.319) The deviation from a respective intrinsic phase-average value is (2.320) When the products of two variables are phase-averaged, the following relations are needed: Φk Ψ k k ˆ Φk ≡ Φk − Φk k = Φk k k Ψk Ψk k ˆˆ + Φk Ψk k k (2.321) (2.322) Φk Ψ k = ε k Φk ˆˆ + Φk Ψ k ΔV nk dAk Figure 2 15 Control volume for volume averaging. Chapter 2 Ge era i ed G ver i g Equati s 147 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press In order to obtain the volume-averaged governing equations, the volume average of the partial derivative with respect to time and gradient must be obtained. For a control volume ΔV shown in Fig. 2.15, the volume averaging of the partial derivative with respect to time is obtained by the following general transport theorem: ∂ Ωk ∂Ωk 1 = − ∂t ∂t ΔV  Ak Ωk VI ⋅ n k dAk (2.323) where Ak is the is the interfacial area surrounding the kth phase within control volume ΔV , ΔVk is the volume occupied by the kth phase in the control volume and ΔV , VI is the interfacial velocity, and nk is the unit normal vector at the interface directed outward from phase k (see Fig. 2.15). The volume average of the gradient is ∇Ωk = ∇ Ωk + 1 ΔV 1 ΔV  Ak Ωk n k dAk (2.324) and the volume average of a divergence is ∇ ⋅ Ω k = ∇ ⋅ Ωk +  Ak Ωk ⋅ n k dAk (2.325) The general quantity Ωk in eqs. (2.323) and (2.324) can be a scalar, vector, or tensor of the second order. It can be a vector or tensor of the second order in eq. (2.325). The formulation of macroscopic equations for multiphase systems can be classified into two groups: (1) the multi-fluid model (Section 2.4.2), and (2) the homogeneous model (Section 2.4.3), also known as the mixture or diffuse model. If the averaging is performed for each individual phase within a multiphase control volume, as shown in eqs. (2.316) and (2.317), one obtains the multifluid model, in which Π sets of averaged conservation equations – each set includes continuity, momentum and energy equations – describe the flow of a Π − phase system. The equations will also include source terms that account for the transfer of momentum, energy, and mass between phases. If only two phases are present, the multifluid model is referred to as the two-fluid model. However, if spatial averaging is performed over both phases simultaneously within a multiphase control volume, as indicated in eq. (2.315), the homogeneous model is obtained; in this case the mixture of a two-phase fluid would be considered a whole. The governing equations for the homogeneous model comprise a single set of equations including continuity, momentum, and energy equations, with one additional diffusion equation to account for the concentration change due to interphase mass transfer by phase change. Continuity, momentum, and energy equations for the mixture model can be obtained by adding together the governing equations for the multifluid models; a diffusion model must be developed to account for mass transfer between phases. In this section, it is assumed for the sake of simplicity that the reference frame is stationary. 148 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2 4 1 2 agra gia Averagi g Lagrangian averaging is directly related to the Lagrangian description of a system, which requires tracking the motion of each individual fluid particle. Therefore, Lagrangian averaging is a very useful tool when the dynamics of individual particles are of interest. To obtain Lagrangian time averaging, it is necessary to follow a specific particle and observe its behavior for a certain time interval. Then, the behavior of this particle is averaged over the time interval. For a generalized function Φ = Φ ( X , Y , Z , t ) , X, Y, and Z are material coordinates moving with the particle, and X, Y, Z are functions of the spatial coordinates x, y, z, and time t, i.e., X = X ( x , y, z , t ), Y = Y ( x , y, z , t ), Z = Z ( x , y, z , t ) The most widely used Lagrangian averaging is time averaging, where the time average of the function Φ in time interval of Δt is Φ= 1 Δt  Δt Φ ( X , Y , Z , t )dt (2.326) Lagrangian time averaging is performed for a distinct particle moving in the field; therefore, X, Y, and Z in the time interval Δt are not fixed in space. This focus on specific particles moving in space and time distinguishes Lagrangian averaging from Eulerian time averaging, which treats a fixed point in space relative to the reference frame. An example from daily experience will serve to illustrate this difference. In order to monitor traffic on the highway, the speed of all cars passing a point can be measured and averaged over a certain time interval – a case of Eulerian averaging. To catch an individual speeder, the police must follow the vehicle of interest to measure its speed as it moves in space over a certain time interval – a case of Lagrangian averaging. 2413 ecu ar Statistica Averagi g When the collective mechanics of a large number of particles is of interest, molecular statistical averaging may be employed. This relies on the concept of particle number density, which is the number of particles per unit volume. For a system with a large number of particles, the behavior of each individual particle is random because random collisions occur. To describe the behavior of each particle, it is necessary to track the motion resulting from each collision – an impractical and often unnecessary task. Although the behavior of each particle is random, the collection of particles may demonstrate some statistical behaviors that are different from those of the individual particles. When the number of molecules involved in the averaging process is large enough, the statistical average value becomes independent of the number of molecules involved. The statistical average value of the microscopic properties for a large number of molecules is related to the macroscopic properties of the system. For example, temperature is a statistical measure of the kinetic energy of the individual Chapter 2 Ge era i ed G ver i g Equati s 149 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press molecules, and the pressure of a gas in a container is the result of many molecules’ collisions with the wall. For some engineering problems, the macroscopic properties of the fluid as well as the microscopic properties are required for design or analysis. Most numerical codes are based on the Navier-Stokes equations, which treats a fluid as a continuous field. It is well known that a fluid is made of a discrete number of particles or molecules. Since the number of molecules is extremely large (Avogadro’s number = 6.022×1023 atoms/mole) for almost all practically sized systems, it may never be computationally viable to track each particle and its interactions with other particles. The number of molecules in a given region and the molecular interaction are described through the fluid’s density and transport coefficients (i.e., viscosity) in the continuous model. Modeling the individual molecules for a small system over a small period of time has been achieved by molecular dynamic simulations (MDS). The computational requirements needed in these simulations can be greatly reduced if the degrees of freedom of the system are reduced. Also, instead of considering individual molecules, groups of molecules can be considered. The degrees of freedom can be reduced by restricting the movement of the molecules to a lattice. A lattice is simply a predefined direction in which a molecule can move. From this standpoint the independent variables are space, velocity and time, while the dependent variable is a molecular distribution function for species i, fi ( x, c, t ) . The Boltzmann equation relates the distribution function at ( x, c, t ) to the distribution function at (x + Δx, c + Δc, t + Δt ). The location in space is x, and the particle velocity is c. It is important to note that the particle velocity is directly related to the mass average velocity, V, that is used throughout this book. This distribution function can be related to the Navier-Stokes equations as well as other transport equations; these relationships give insight to the origin of transport coefficients such as viscosity. A detailed presentation of Boltzmann statistical averaging including the discussion of Lattice Boltzmann model for both single and multiphase systems can be found in Faghri and Zhang (2006). 2 4 2 V u e Averaged u ti F uid de s If spatial averaging is performed for each individual phase within a multiphase control volume, the multi-fluid model is obtained. Additional source terms are needed in these equations to account for the interaction between phases. C ti uity Equati The volume average of the continuity equation for the kth phase is obtained by taking extrinsic phase averaging from eq. (2.51) 150 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂ρk + ∇ ⋅ ρk Vk = 0 ∂t (2.327) where the two terms on the left-hand side can be obtained by using eqs. (2.323) and (2.325), i.e., ∂ ρk ∂ρk 1 = − ρk VI ⋅ nk dAk ∂t ∂t ΔV  Ak ∇ ⋅ ρk Vk = ∇ ⋅ ρk Vk + 1 ΔV  Ak ρk Vk ⋅ n k dAk Substituting the above expressions into eq. (2.327), the volume-averaged continuity equation becomes ∂ ρk 1 ρk (Vk − VI ) ⋅ n k dAk (2.328) + ∇ ⋅ ρk Vk = − ∂t ΔV  Ak The right-hand side of eq. (2.328) represents mass transfer per unit volume from all other phases to the kth phase due to phase change; it can be rewritten as − 1 ΔV  Ak ρk (Vk − VI ) ⋅ n k dAk = j =1( j ≠ k )  Π  jk m′′′ (2.329)  ′′ where m′jk represents mass transfer per unit volume from the jth to the kth phase  ′′ due to phase change. The value of m ′jk depends on the phase change process that takes place in the multiphase system, and the conservation of mass requires  ′′  ′′′ that m ′jk = −mkj . The extrinsic phase-averaged density, averaged density, ρk k ρk , is related to the intrinsic phase- , by ρk = ε k ρk k (2.330) Furthermore, the intrinsic phase-averaged density is equal to the density ρk . Substituting eqs. (2.329) and (2.330) into eq. (2.328), and considering eq. (2.322), the continuity equation for the kth phase becomes ∂ ε k ρk ∂t ( k ) + ∇ ⋅ ε k ρk ( k Vk k ˆˆ + ρk Vk ) = Π  jk m′′′ (2.331) j =1( j ≠ k ) The dispersive term in eq. (2.331), k k ˆˆ ρk Vk , is generally small compared with ε k ρk Vk ; it is assumed that it can be neglected. The continuity equation for the kth phase becomes ∂ ε k ρk ∂t ( k ) + ∇ ⋅ (ε k ρk k Vk k )=  Π  jk m ′′′ (2.332) j =1( j ≠ k ) Chapter 2 Ge era i ed G ver i g Equati s 151 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press e tu Equati The extrinsic phase-averaged momentum equation for the kth phase can be obtained by performing extrinsic phase-averaging on the momentum equation (2.65): ∂ ( ρk Vk ) + ∇ ⋅ ( ρk Vk Vk ) = ∇ ⋅ τ ′ + ρk Xk (2.333) k ∂t where the body force per unit mass is assumed to be the same for different species for sake of simplicity. After evaluating each term in eq. (2.333), the multi-fluid volume-averaged momentum equation becomes (Faghri and Zhang, 2006) ∂ ε k ρk ∂t ( k Vk k k ) + ∇ ⋅ (ε k k ρk k Vk Vk Π k ) = ∇ ⋅ ε k τ′ k ( ) + ε k ρk Xk + j =1( j ≠ k )  (  jk F jk + m ′′′ Vk , I k ) (2.334) where Vk , I k is intrinsic phase-averaged velocity of the kth phase at the interface. The difference between two adjacent phases results solely from the density difference between the two phases. F jk is an interactive force between the jth and the kth phase, and depends on the friction, pressure, and cohesion between different phases. Newton’s third law requires that the interactive forces satisfy (2.335) F jk = − Fkj The interactive force can be determined by F jk = K jk where Kjk is the momentum exchange coefficient between phases j and k. Determining the momentum exchange coefficient is a very challenging task because interphase momentum exchange depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k, as is the case with the dispersed phase system summarized in Table 1.8, one can assume that the secondary phase is spherical in shape and an appropriate empirical correlation can be used to obtain the momentum exchange coefficient. Since liquid-vapor flow is widely used in various applications, we will use liquid-vapor flow as an example to explain the determination of the momentum exchange coefficient. If liquid is the primary phase and vapor is the secondary phase, the vapor phase is dispersed in the liquid as vapor bubbles. If vapor is the primary phase and liquid is the secondary phase, the liquid phase is dispersed in the vapor as liquid droplets. Boysan (1990) suggested that the momentum exchange coefficient could be estimated by ε j ρk j 3 k K jk = C D V j − Vk (2.337) 4 dj (V j j − Vk k ) (2.336) 152 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where phase k is the primary phase, and phase j is the secondary phase, and dj is the diameter of vapor bubbles or liquid droplets of the secondary phase j. CD is the drag coefficient based on the relative Reynolds number, which obtained by the following empirical correlations:  24 0.687 )  (1 + 0.15 Re C D =  Re 0.44  Re ≤ 1000 Re > 1000 k (2.338) where Re = ρk Vj j − Vk dj μk (2.339) Example 2.9 A mixture of water and its vapor at 1 atm flows in a 0.1 m ID tube with a mass flow rate of 0.1 kg/s. The liquid water is dispersed in the vapor phase in the form of 0.1- mm -diameter droplets and the quality of the mixture is x=0.8. While the vapor is saturated, the liquid droplets, are subcooled, both at 95 °C. The volume fraction of the liquid phase is ε  = 0.01 . Find the interactive force between the liquid and vapor phases. Solution: Since the liquid droplets are dispersed in the vapor phase, vapor ( v ) is the primary phase and liquid (  ) is the secondary phase. At 1 atm, the properties of the saturated vapor are μ v = 16.698 × 10−6 N-s/m 2 ρ v = 0.596 kg/ m3 The density of the subcooled water at 95 °C is ρ = 965.35kg/m3 The mass flow rates of the liquid and vapor phases are, respectively,   mv = xm = 0.8 × 0.1 = 0.08kg/s   m = (1 − x )m = (1 − 0.8) × 0.1 = 0.02kg/s The total cross-sectional area of the tube is 1 1 A = π D 2 = × π × 0.12 = 7.85 × 10−3 m 2 4 4 The mass flow rates of the liquid and vapor phases are, respectively,   m = ρ w A = ρ w ε  A  mv = ρ v wv Av = ρ v wv v (1 − ε  ) A The average velocities of the liquid and vapor phases can be obtained by w wv  = v 0.02 = 0.264m/s 965.35 × 0.01× 7.85 × 10−3  mv 0.08 = = = 17.27m/s ρ v (1 − ε  ) A 0.596 × (1 − 0.01) × 7.85 × 10−3  m ρ ε  A = The relative Reynolds number of the liquid droplets is Chapter 2 Ge era i ed G ver i g Equati s 153 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Re = = ρ v w  − wv v d μv 16.698 × 10−6 = 60.70 0.596 × 0.264 − 17.27 × 0.1×10−3 The drag coefficient is then CD = 24 24 (1 + 0.15 Re0.687 ) = × (1 + 0.15 × 60.700.687 ) = 1.39 Re 60.70 The momentum exchange coefficient can be obtained from eq. (2.337), i.e., v ε  ρv 3 v  K v = = 4 CD d w − wv 3 0.01× 0.596 × 1.39 × × 0.264 − 17.27 = 1056.63kg/(m3 -s) −3 4 0.1× 10 The interactive force between liquid droplets and vapor is Fv = K v (w   − wv v ) = 1056.63 × ( 0.264 − 17.27 ) = −1.797 × 104 N/m3 which is negative because the liquid droplet moves slower than the vapor phase. E ergy Equati The extrinsic phase-average of the energy equation, (2.92), is ∂ ( ρk hk ) + ∇ ⋅ ρk Vk hk ∂t ′′′ = − ∇ ⋅ q′′ + qk + k ∂pk + Vk ⋅∇pk + ∇Vk : τ k ∂t (2.340) which can be used to obtain the volume-averaged energy equation and the result is (Faghri and Zhang, 2006) ∂ ε k ρk ∂t +ε k D pk Dt ′′ q′jk ( k k hk k ) + ∇ ⋅ (ε k ρk k Π Vk hk k ) = −∇ ⋅ q′′ k ′′′ + qk k + ∇ Vk : τ k + j =1( j ≠ k )   q′′′ + m ′′′ h  jk k , I  jk  (2.341)   where is the intensity of heat exchange between phase j and k. It can be hc ΔA j obtained by using Newton’s law of cooling: q′′′ = jk (T j j − Tk k ΔV j ) (2.342) 154 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where hc is the convective heat transfer coefficient, ΔAj is the area of the interface between phases j and k, and ΔVj is the volume of the secondary phase in the volume element. Like the momentum exchange coefficient, the interphase heat transfer also depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k – as in the dispersed phase summarized in Table 1.8– the following empirical correlation recommended is widely used: Nu = 2 + (0.4 Re 1/2 + 0.06 Re 2/3 0.4 ) Prk  μk  μ  k ,s     1/4 (2.343) where the Reynolds number, Re, is obtained by eq. (2.339), the Nusselt number is defined as Nu = hc d j kk (2.344) Tk except and all thermal properties of the primary phase are evaluated at μk , s , which is evaluated at T j . Equation (2.343) is valid for 3.5 < Re < 7.6 × 104 and 0.71 < Prk < 380 , which covers a wide variety of problems. If the secondary phase is liquid and the primary phase is vapor (gas), eq. (2.343) can be simplified to 1/3 Nu = 2 + 0.6 Re1/2 Prk (2.345) Example 2.10 Estimate the interphase heat transfer rate between the liquid and vapor phases in Example 2.9. Solution: The thermal conductivity and the Prandtl number of the vapor phase are kv = 0.0248W/m-K and Prv = 0.984 , respectively. Since the primary phase is vapor ( v ) and the secondary phase is liquid (  ), eq. (2.345) can be used to estimate the heat transfer coefficient, i.e., 1/3 Nu = 2 + 0.6 Re1/2 Prv = 2 + 0.6 × 60.701/2 × 0.9841/3 = 6.65 The heat transfer coefficient is then hc = Nukv 6.65 × 0.0248 = = 1649.2W/m 2 -K −3 d 0.1× 10 Therefore, the interphase heat transfer between the liquid and vapor phases can be obtained by eq. (2.342): ′′′ qv = hc A T (  − Tv v v = 6hc T ( V  ) = h πd ( T c 2    − Tv v − Tv d ) = 6 ×1649.2 × (95 − 100) = −4.95 ×10 W/m 8 3 π d ) /6 3 0.1× 10−3 Chapter 2 Ge era i ed G ver i g Equati s 155 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Species If the fluid undergoing phase change involves multiple components, it is also necessary to write the equation for conservation of the species mass in the kth phase. The extrinsic phase-average of conservation of species mass can be obtained by ∂ρk ,i  ′′′ + ∇ ⋅ ρk ,i Vk = − ∇ ⋅ J k ,i + mk ,i (2.346) ∂t where each term can be obtained by ∂ ρk ,i ∂ρk ,i 1 ∂t = ∂t − ΔV ∇ ⋅ ρk ,i Vk = ∇ ⋅ ρk ,i Vk + ∇ ⋅ J k ,i = ∇ ⋅ J k ,i  Ak ρk ,i VI ⋅ n k dAk 1 ΔV 1 + ΔV  Ak ρk ,i Vk ⋅ n k dAk  Ak J k ,i ⋅ n k dAk Substituting the above expression into eq. (2.346), one obtains ∂ ρk ,i  ′′′ + ∇ ⋅ ρk ,i Vk = −∇ ⋅ J k ,i + mk ,i ∂t 1 − ΔV  1 ρk ,i (Vk − VI ) ⋅ n k dAk + Ak ΔV  (2.347) Ak J k ,i ⋅ n k dAk where the third and fourth terms in the right-hand side of eq. (2.347) represent mass source (or sink) of the ith component in the kth phase due to phase change from other phases to the kth phase, as well as mass transfer at the interface due to diffusion, respectively, i.e., − 1 ΔV  Ak ρk ,i (Vk − VI ) ⋅ n k dAk + 1 ΔV  Ak J k ,i ⋅ nk dAk =  jk where m ′′′ ,i represents the mass source (or sink) of the i component in phase k j =1( j ≠ k ) th  Π  jk m ′′′ ,i (2.348) due to phase change from phase j to phase k, as well as diffusive mass transfer at the interface between phases j and k. Substituting eq. (2.348) into eq. (2.347) and considering eq. (2.322), the volume-averaged species equation becomes ∂ ε k ρk ,i ∂t ( k ) +∇⋅ ε ( k ρk ,i k Vk k )  ′′′ = −∇ ⋅ J k ,i + mk ,i + j =1( j ≠ k )  Π  jk m ′′′ ,i (2.349) where the three terms on the right-hand side represent the effects of mass diffusion, species source/sink due to chemical reaction, and phase change at the interfaces. Example 2.11 A cooling tower is a multiphase system. Warm water is sprayed at the top, in the form of water droplets. Gravity drives these droplets 156 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press y, vk k x, z, Water in uk k Air/water vapor wk k g Domain to model Insulated wall Air in Water Figure 2 16 Schematic of the cooling of droplets. Air in downward. The water droplets are cooled by evaporated effects as they flow through the air. Air enters the cooling towers from the bottom on the sides. The system being investigated has a planar geometry. The schematic for this system is presented in Fig. 2.16. The temperature of the water collected at the bottom of the tank is desired, as well as the amount of water lost due to evaporation. Develop the governing equations by the volume-average method, as well as the boundary conditions. Solution: This problem can be solved using a volume-averaged multifluid formulation. The basic assumptions are that the flow field is steady. The planar geometry is long enough that the changes in the zdirection can be neglected; therefore a two-dimensional geometry is utilized. Also, since this is the case, a plane of symmetry exists half the distance between the two adiabatic walls. The pressure effects on energy are considered negligible. The final assumption is that intrinsic phase average of the products of deviation is zero, ˆˆ Vk Vk k ˆˆ = Vk hk  k ˆˆ = Vk ωk  k = 0.   The volume-averaged continuity equation, eq. (2.332), for two phases is: ∂ ε  ρ ∂x ( u ) + ∂∂y (ε  ρ v  ) = m′′′ v (2.350) Chapter 2 Ge era i ed G ver i g Equati s 157 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂ ε v ρv ∂x ( v uv v ) + ∂∂y (ε  v ρv v vv v   ) = m′′′ = −m′′′ v v (2.351) The volume-averaged momentum equations for both phases in the x- and y- directions are obtained from eq. (2.334): Liquid x-momentum ε  ρ −  u   ∂ u ∂x  + ε  ρ  v  ∂ u ∂y  =  ∂ε  p ∂x + ∂ u ∂  ε  μ ∂x  ∂x   ∂ ∂ u  +  ε  μ  ∂y  ∂y     + Fv, x + mv  ′′′   u , I  (2.352) Liquid y-momentum ε  ρ =−  u  ∂ v ∂x  + ε  ρ  v  ∂ v ∂y  ∂ε  p ∂y  + ∂ v ∂  ε  μ ∂x  ∂x  v, I v    ∂ ∂ v  +  ε  μ  ∂y  ∂y         (2.353)  ′′′ + Fv, y + mv − ε  ρ g Vapor x-momentum ε v ρv + v uv v ∂ uv ∂x v + ε v ρv v vv ∂ uv ∂y v v =− ∂ε v pv ∂x v ∂ uv ∂  ε v μv ∂x  ∂x   ∂ ∂ uv  +  ε v μv  ∂y  ∂y     + Fv, x + mv  ′′′   uv, I v (2.354) Vapor y-momentum ε v ρv =− v uv v ∂ vv ∂x v + ε v ρv v vv v ∂ vv ∂y v ∂ε v pv ∂y v ∂ vv ∂ +  ε v μv ∂x  ∂x  vv, I v v  ∂ ∂ vv  +  ε v μv  ∂y  ∂y   v v     (2.355)  ′′′ + Fv, y + mv − ε v ρv g The volume-averaged energy equations for the liquid and vapor phases are obtained from eq. (2.341): Liquid energy 158 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ε  ρ  u ∂ h ∂x   + ε  ρ  v  ∂ h ∂y  = ∂ T ∂  ε  k ∂x  ∂x   ∂ ∂ T  +  ε  k  ∂y  ∂y   v v    + qv + mv hv  ′′′ ′′′   v (2.356) Vapor energy ε v ρ v uv v ∂ hv ∂x v + ε v ρv vv ∂ hv ∂y v = ∂ Tv ∂  ε v kv ∂x  ∂x   ∂ ∂ Tv  +  ε v kv  ∂y  ∂y    v   + qv + mv hv  ′′′ ′′′   (2.357) where the enthalpy is defined as: h hv v  = c N = c p dT pv dT (2.358) (2.359) There are multiple components in the vapor phase; therefore, the specific heat is a function of the mass fraction of those components. c pv = ω c i =1 i pi (2.360) Air is considered to be one component, and water vapor is considered to be the other component. Therefore, only one species equation is needed. The species equation for water vapor, eq. (2.349), is: v v v v ∂ ωv v v ∂ ωv ε v ρ v uv vv + ε v ρv = ∂x ∂y v ∂ ωv ∂  ε v Dv,air ∂x  ∂x  v  ∂ ∂ ωv  +  ε v Dv,air  ∂y  ∂y     + mv  ′′′   (2.361) It is also known that ε + εv = 1  ′′′ mv = (2.362) (2.363) − hm Av (ωsat − ωv,i ) ΔV (1 − ωsat ) v ωsat = p0 Rg Tv v ρv h exp  − v  Rg  1  T v v 1 T0     (2.364) (2.365)  ′′′ qv = hc Av ( Tv − T ) ΔV Fv = − Fv = K v (V v v − V  ) (2.366) Chapter 2 Ge era i ed G ver i g Equati s 159 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ρv v = P Rg Tv v (2.367) The surface to volume ratio can be calculated by assuming that the liquid is in the form of spherical droplets. Therefore, the surface to volume ratio is defined in terms of the liquid droplet diameter, d , and the volume fraction of the liquid. Av 6 = ε d ΔV (2.368) There are 15 equations and 16 unknowns (9 equations and 6 constitutive relations). To solve this problem, it is assumed that the liquid is in droplets, and there is no pressure gradient within each droplet. Therefore the liquid pressure gradient is assumed to be equal to the vapor pressure gradient. ∂ p ∂x  = ∂ pv ∂x v , ∂ p ∂y  = ∂ pv ∂y v (2.369) Now the liquid pressure is eliminated, therefore we have 13 equations and 13 unknowns, which is well posed. The boundary conditions are: At walls uv v = vv  v = u v  = v v  =0 =0 (2.370) (2.371) (2.372) (2.373) (2.374) (2.375) (2.376) (2.377) ∂ T ∂x = ∂ Tv ∂x v = ∂ ωv ∂x v At air inlet pv = ρv v gh Tv = Tin ε = 0 At top T  = Ttop v ε  = ε ,in pv =0 k At symmetry ∂φ = 0 , φ = vk ∂x uk k k , Tk , ε  , pk k , ωv v (2.378) (2.379) (2.380) = 0 , k = , v At bottom surface v ωv = ωv,sat 160 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press − ε v ρv v Dv,i ,air ∂ ωv v v v 1 − ωv ∂y = ε v ρv v v vv v (2.381) (2.382) (2.383) (2.384) −ε v kv ∂ Tv ∂y uv = ε v ρv v vv v hv = u v  ≈0  Tv = T 2 4 3 V u e Averaged H ge e us de The multi-fluid model presented above is obtained by performing phase averaging as defined in eqs. (2.316) and (2.317). If spatial averaging is performed for all phases within a multiphase control volume, the homogeneous (or mixture) model can be obtained. The relationship between volume averaging and phase averaging is given in eq. (2.319), which indicates that the homogeneous model can be obtained by summing the individual phase equations of the multi-fluid model. C ti uity Equati The continuity equation for phase k in the multifluid model is expressed by eq. (2.332). Summing the continuity equations for all Π phases together, one obtains ∂  ε k ρk ∂t  k =1   Π k Π   + ∇ ⋅ ε k ρk  k =1   k Vk k =  Π Π  jk m ′′′ (2.385) k =1 j =1( j ≠ k ) The right-hand side of equation (2.385) must be zero because the total mass of all phases produced by phase change must equal the total mass of all phases consumed by phase change. Considering this fact and eq. (2.330), the continuity equation becomes Π ∂ρ k k Vk = 0 (2.386) + ∇ ⋅ ε k ρk ∂t  k =1 Π The bulk velocity of the multiphase mixture is the mass-averaged velocity of all the individual phases: 1  V= ρ ε k =1 k ρk k Vk k (2.387) Substituting eq. (2.387) into eq. (2.386), the final form of the continuity equation for a multiphase mixture is Chapter 2 Ge era i ed G ver i g Equati s 161 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂ρ ∂t  +∇⋅ ρ V = 0 ( ) (2.388) It can be seen that eq. (2.388) has the same form as the local continuity equation (2.51), except that the volume-averaged density and velocity are used in eq. (2.388), where ρ= ε k =1 Π k ρk k . e tu Equati The momentum equation for phase k in the multi-fluid model is expressed in eq. (2.334). By adding together the momentum equations for all Π phases, one obtains ∂  ε k ρk ∂t  k =1   Π Π k Vk + Π k  Π  + ∇ ⋅  ε k ρk     k =1  k k Vk Vk Π k     jk = ∇⋅ ε k =1 k τ′ k k  (ε k =1 k ρk )X +   ( F Π k =1 j =1( j ≠ k )  jk + m ′′′ Vk , I k ) (2.389) The stress tensor of the multiphase mixture is τ′ = ε k =1 Π k τk k 2    = − p I + μ ∇V + ∇V T  − μ (∇ ⋅ V )I  3 (2.390) The summation of all interphase forces must be zero since F jk = − Fkj , i.e.,  ∂ ∂t  ( ρ V ) + ∇ ⋅  εk   k =1 Π Π Π F jk = 0 (2.391) k =1 j =1( j ≠ k ) Considering eqs. (2.387), (2.390) and (2.391), the momentum equation becomes  Π ρk k Vk Vk k    = ∇ ⋅ τ ′ + ρ X + M ′′′ I   Vk , I k (2.392) where  M ′′′ = I  Π  jk m′′′ (2.393) k =1 j =1( j ≠ k ) Equation (2.393) represents the momentum production rate due to interaction between different phases along their separating interfaces. It must be specified according to the combination of phases in the multiphase system that is under consideration. 162 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press E ergy Equati By summing the energy equations for all Π phases in the multifluid model, eq. (2.341), one obtains ∂  ε k ρk  ∂t  k =1  εk Π k hk + Π k  Π  + ∇ ⋅  ε k ρk     k =1  k Vk hk Π k  Π q′′  = −∇ ⋅  k     k =1  Π Π ′′′ qk +   k =1  +  k =1 Π D pk Dt k  k =1 ∇ Vk k : ε k τk k +  k Π q′′′ + jk k =1 j =1( j ≠ k )  Π  jk m ′′′ hk , I k k =1 j =1( j ≠ k ) (2.394) The mass average enthalpy of the multiphase mixture is 1  h= ρ ε k =1 Π k ρk k hk (2.395) The fifth term on the right-hand side of eq. (2.394) is for summation of all interphase heat transfer and it must be zero. The last term on the right-hand side of eq. (2.394) accounts for contribution of interphase phase change energy flux due to phase change; it can be defined as  Π Π  jk m ′′′ hk , I k ′′′ = qI  m ′′′ = 0 . (2.396) k =1 j =1( j ≠ k ) It is usually not zero although  k Π Π k =1 j =1( j ≠ k ) Considering eqs. (2.395) and (2.396), the energy equation (2.394) becomes ∂ ∂t  ( ρ h) + ∇ ⋅ ε   k =1  Π k ρk k Vk hk  Dp  ′′′ + q′′′ + ∇V : τ + qI  = −∇ ⋅ q′′ +  Dt  (2.397) Species Summing the equations for conservation of species mass, eq. (2.349), for all phases yields ∂  ε k ρk ,i ∂t  k =1   Π k  Π  + ∇ ⋅  ε k ρk ,i     k =1  k Vk Π k     Π J k ,i = −∇ ⋅    k =1  Π Π  ′′′ mk ,i + +  k =1  k =1 (2.398)    jk m ′′′ ,i j =1( j ≠ k ) By applying eq. (2.330) to the mass density of the ith component, one obtains Chapter 2 Ge era i ed G ver i g Equati s 163 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ρi = th ε k =1 Π k ρk ,i k (2.399) In accordance with the conservation of mass, the mass source (or sink) of the i component due to phase change in all phases must add up to zero, i.e.,  Π Π  jk m ′′′ ,i = 0 (2.400) k =1 j =1( j ≠ k ) Substituting eqs. (2.399) and (2.400) into eq. (2.398), and using the massaveraged velocity defined in eq. (2.387), the conservation of species mass becomes ∂ ρi   (2.401) + ∇ ⋅ ρi V = −∇ ⋅ J i + mi′′′ ∂t Example 2.12 Solve problem in Example 2.11 using the homogeneous model using the same assumptions made in Example 2.11. The water droplets are small; therefore, make further reasonable assumptions to find a relationship between the liquid velocity and vapor velocity so that the homogeneous model can be applied. The flow variable to be solved is the   mass averaged velocities, u and v , as well as the average temperature, T. Solution: Use the same assumptions listed in Example 2.10. To make this problem a homogenous problem, more assumptions need to be made. The inertial terms in the liquid momentum equation can be assumed to be negligible because the droplets are very small. It can also be assumed that the evaporation rates and pressure gradients do not affect liquid momentum equation. Also, assume that the temperature of the liquid and  v vapor are in local equilibrium, T = Tv . Based on the assumptions above, the liquid momentum equation in the x-direction reduces to a statement that the velocities of the liquid and vapor are equivalent. (2.402) Similarly, the liquid momentum equation in the y-direction reduces to a statement that drag force between the liquid and vapor is equivalent to the gravitational force on the liquid droplets.  ερ g  v (2.403) v = vv −   u = uv K v  v  =u Therefore, the y component of the mass averaged velocity is 164 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (ε  v= v ρv v + ε  ρ  (2.404) Therefore, the continuity equations, eq. (2.388), for the liquid phase and the mixture are: Liquid 2  ε2 ρ g ε  ρ ∂ ∂   ) +  ε  ρ v +    ( ε  ρ u  ∂x ∂y  K v ρ    (ε )v v v v 2 − ε  ρ ( 2 v ρv + ε  ρ  ) ) −1 gK v = vv v − 2 ε  ρ ( 2 ρ K v ) g   ′′′ − 1  = mv   (2.405) Mixture ∂ ∂   ( ρ u ) + ∂y ( ρ v ) = 0 ∂x (2.406) The homogeneous momentum equations, eq. (2.392), are:  2  2 ε  ρ g  ∂u   v  εv v+  ρ K v  ∂y    ∂p   ∂ ∂u  ∂  ∂u  =− +  ε v μv   +  ε v μv ∂x ∂x  ∂x  ∂y  ∂y    2  2  2 2 ε  ρ g ε  ρ g ∂ v v    ρv u  v +  + ε v ρv  v +  ρ K v  K v ∂x         v ∂u  ρv u + ε v ρv ∂x ( ) (2.407) εv ( ) ( )    2   2  2 ε  ρ ε l2 ρ g g  ∂p ∂ ∂ ∂   ⋅ v+ +  ε v μv =− v+  ρ K v  ρ K v   ∂y  ∂y ∂x  ∂x            2  ε l2 ρ g    ∂ ∂  (2.408) +  ε v μv v+   − ε v ρv g ∂y  ∂y  ρ K v        ( ) ( ) ( ) The homogenous energy equation, eq. (2.397), is:   u  ε v ρv   +ε  ρ  v ∂ hv ∂x v + ε  ρ  ∂ h ∂x    +ε ρ vv   2  2 ε  ρ g∂ h  v v  v+  ρ K v  ∂y    ( ) v   ε ρ g  ε  ρ v+       K v ρ     ∂h ∂T  ∂ − 1   =  ( ε v kv + ε  k )    ∂y ∂x  ∂x   + ∂T  ∂  ′′′  ( ε v kv + ε  k )  + mv hv ∂y  ∂y  (2.409) Chapter 2 Ge era i ed G ver i g Equati s 165 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The enthalpies of the liquid and the vapor are: h hv  = v c = c N p dT pv dT (2.410) (2.411) There are multiple components in the vapor phase; therefore the specific heat is a function of the mass fraction of those components. c pv = ω c i =1 i pi (2.412) The homogenous species equation, eq. (2.401), is:  2  2 ε  ρ g∂ ω v  v v v   + ε v ρv  v + ε v ρv u  ∂x ρ K v  ∂y    v v ∂ ωv ∂ ωv ∂ ∂  ε v Dv,air  +  ε v Dv,air  + mv  ′′′ = ∂x  ∂x  ∂y  ∂y      ∂ ωv v ( ) (2.413) Also ε + εv = 1  ′′′ mv = hm Av ωsat − ωv ΔV (1 − ωsat ) v ( v ) v (2.414) (2.415) − ωsat = p0 Rg Tv ρv v v ρv = h exp  − v  Rg  P v 1  T v 1   T0    (2.416) (2.417) Rg Tv The surface to volume ratio can be calculated by assuming that the liquid is in the form of spherical droplets. Therefore, the surface to volume ratio is defined in terms of the liquid droplet diameter, d , and the volume fraction of the liquid. Av 6 = ε ΔV d (2.418) The problem was reduced to 6 equations with 3 constitutive relations from 9 equations. In total, the present formulation has 9 equations and 9 unknowns and is well posed. At walls  u=v=0 (2.419) ∂T ∂x = ∂ ωv ∂x v =0 (2.420) At top 166 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press T = Ttop p =0 (2.421) (2.422) (2.423) ε  = ε ,in At air inlet p = ρv v gh (2.424) (2.425) (2.426) v T = Tin ε = 0 At symmetry ∂φ  = 0, φ = v, T , ε  , ωv ∂x  u=0 (2.427) (2.428) (2.429) ∂ ωv v At bottom surface ωv  ρ v = − ε v ρv v = ωv,sat  ( v + ε  ρ )1− ω Dv,air v v ∂y − ε  ρ  ε  ρ K v  g (2.430) (ε v kv + ε  k ) ∂T ∂y =εv  2  2 ε  ρ g  v  ρv  v + h ρ K v  v    ( ) (2.431) From the above problem formulation, the homogenous model can still be used for counter-flowing phases. However, more assumptions were made to solve the problem as a homogeneous model. So the trade off for reducing the number of equations is the accuracy of the solution. 2 4 4 A Exte si P r us edia Transport in porous media is applicable to a wide range of fields, including mechanical, chemical, environmental, and petroleum engineering, as well as geology. Porous media can be found naturally in rocks and sand beds, and also can be fabricated as wicks and catalytic pellets. They are essential component in high-technology devices such as fuel cells and heat pipes. A fundamental formulation of the governing equations within a porous media will be presented here. The variation of some of the assumptions made has been quantified in recent studies. For these analyses of the variants within a particular porous media model, see Alazmi and Vafai (2000, 2001), and Nield and Bejan (1999). A porous medium is a solid matrix with a several voids, or pores, which are continuously connected. The voids are filled with one or more fluids that can pass through the medium because the voids are interconnected. The void fraction in the solid matrix is most frequently referred to as the porosity, ε . Chapter 2 Ge era i ed G ver i g Equati s 167 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Conversely, 1 − ε is the volume fraction of the solid matrix, ε sm , that makes the porous medium. If g( x , y, z , t ) is a geometric function that is identically equal to 1 in a void, and 0 in the solid matrix, the porosity is defined as: ε= 1 ΔV  ΔV gdV . (2.432) ΔV is an elemental volume of the porous zone. When there is more than one phase present in a porous medium, there is a porosity of a single phase, ε k , which is simply the volume fraction of a particular phase in an elemental volume. εk = 1 ΔV  ΔVk gdV (2.433) The subscript, k, pertains only to phase k. The volume fraction of the solid matrix in an elemental volume is represented by ε sm . Therefore, 1 = ε + ε sm = ε k =1 Π−1 k + ε sm (2.434) To model a porous zone, it is imperative that one be familiar with the different length scales of the zone. The first dimension is the particle or void length scale, d, and the second dimension, L, is the system or porous zone length scale. If d is on the order of L, such as in a very thin porous layer on a heat transfer surface, where the bulk flow is normal to the surface, the porous zone can be directly modeled with minimal assumptions. However, this situation is usually not the case, and more often d  L . When d  L , direct simulation of transport characteristics in an individual pore is not practical, therefore the local mean velocity is used. Because of this difference in scale, the volume averages defined in eqs. (2.317) and (2.318) are often used to describe a porous system. (a) (b) Figure 2 17 Comparison of velocity distributions in an elemental volume in a porous medium: (a) actual velocity, and (b) volume-averaged velocity. 168 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Figure 2.17 (a) shows how the flow in an individual pore can be very difficult to model when a larger scale is needed, and therefore a volume-averaged velocity, seen in Fig. 2.17 (b), is often more useful. For analyses in a porous system, one of two velocities are often used: the intrinsic average velocity, k Vk , which is the volume-averaged velocity over the volume filled with a particular phase, or the extrinsic average velocity, Vk , which is the extrinsic averaged velocity over an entire volume including a solid. These two velocities are related by the Dupuit-Forchheimer relationship: k Vk = ε k Vk (2.435) The governing equations for the fluid phase in a porous medium can now be developed. It is assumed that the porous medium is saturated with the fluid phase and that the relative velocity, Vf,rel, only refers to the relative motion of a fluid phase f with respect to the solid matrix. In this section, the solid matrix is assumed to be stationary, therefore V f = V f ,rel . The conservation laws will be derived for a porous zone by using a volume-averaging technique for a single fluid occupying the pores. The momentum equations will be derived in a manner in which the viscous interactions between a fluid and the solid matrix are modeled with a property of the porous wick called the permeability, K. The energy equation will be derived in such a way that local thermal equilibrium between the solid matrix and the fluid can be assumed, or local thermal nonequilibrium conditions can be assumed. C servati f ass The fluid-saturated porous media can be considered as a two-phase system that contains fluid (f) and solid matrix (sm). According to eq. (2.332), the continuity equation for the fluid phase (f) in the porous media is ∂ ε ρf ∂t ( f ) + ∇ ⋅ (ε ρ f f Vf f )=0 (2.436) where Vf f is the intrinsic phase-averaged velocity. The right-hand side is zero, because there is no phase change between the porous matrix and the fluid. If the fluid phase in the porous media can be assumed as incompressible, the volume-averaged density is equivalent to the density of a single phase, ρf f = ρ f and eq. (2.436) can be simplified as ∂ ερ f + ∇ ⋅ ρ f V f ∂t ( ) ( )=0 (2.437) where Vf = ε Vf f is the extrinsic phase-averaged velocity. For an incompressible flow in a porous media of uniform porosity, the continuity equation is Chapter 2 Ge era i ed G ver i g Equati s 169 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∇ ⋅ Vf = ∇ ⋅ ε Vf ( f ) = 0. (2.438) C servati f e tu Since a macroscopic approach is the most feasible way to model transport, there needs to be a way to model the bulk resistance to the flow caused by the porous zone. In 1856, Henry Darcy experimentally measured the resistance to a steady, one-dimensional, gravitationally-driven flow through an unconsolidated, uniform, rigid, and isotropic solid matrix. He came up with a relationship for the pressure gradients/resistance as a function of the dynamic viscosity, μf, the extrinsic phase-averaged velocity, and the permeability, K, now known as Darcy’s law. μf 1 ρ g − ∇(ε p) = Vf (2.439) K ε where the permeability, K, is the resistance to the flow, which is analogous to the thermal conductivity in the Fourier’s law of heat conduction. It has units of meters squared. Darcy’s law can be expanded into vectorial form for an anisotropic solid matrix. μf 1 ρ g − ∇(ε p) = Vf (2.440) K ε where the permeability, K, becomes a second-order symmetric tensor. When the porous medium is isotropic (having equal resistance in all directions), the permeability reduces to a scalar (single value), K. Darcy’s law is valid for Re<1, or the creeping flow regime, where the viscous forces dominate. The Reynolds number of the porous media is defined by the mean volumetric velocity, and the average characteristic length scale of the voids, d . ρ f Vf d Re = (2.441) μf In this flow regime, the wall effects are confined to within one or two particle diameters of the wall. Also, in the region where the flow enters the porous zone, the entrance region is confined to approximately three particle diameters. Therefore flow becomes very closely uniform in one main direction, and the walls bounding the porous zone have minimal effect. As the Reynolds number increases from unity to about 10, the drag smoothly transitions from linear to nonlinear. This linear to nonlinear transition of drag in a porous zone represents the flow moving from the creeping flow to the laminar flow regime. To account for this nonlinearity, when Re>10, Ward (1964) came up with a quadratic drag term dependent of ρ/K1/2. μf Cf 1 ρ g − ∇(ε p) = Vf + 1 2 ρ f Vf Vf (2.442) K ε K 170 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where Cf is a dimensionless drag constant, initially thought by Ward (1964) to be a universal constant, 0.55. Since the drag smoothly transitions from linear to nonlinear when 0<Re<10, the significance of the quadratic drag term can be formed as a linear function of Re. μf Cf ( Re− 1) 1 ρ g − ∇(ε p) = Vf + 1 2 ρ f Vf Vf (2.443) 9 K ε K It was later found to vary based on the nature of the porous medium. From the differential momentum equation, eq. (2.65), the advective term is ρ f ( V f ⋅∇ ) V f . From a scaling analysis, the velocity of the fluid is on the order of the volume weighted average velocity over the porosity, for a single phase, V f ~ V f / ε . Therefore, the advective term is on the order of ρ f V f 2 / (ε 2 L ) . A scaling analysis between the advective term and the quadratic drag term gives a ratio proportional to K1 2 / (C f ε 2 L) , where L is the macroscopic characteristic length scale. The ratio of these two terms is usually very small, therefore the advective term often can be assumed to have a negligible effect on the momentum equations. An extension of the Stokes drag force on a sphere in an infinite domain is to include the effects of the neighboring spheres, as with the case of packed beds made of spheres. The Stokes drag force is the drag exerted by a single sphere in an infinite domain on a flow with a Re  1 (negligible inertia). This was accomplished by superimposing the Stokes flow on the Darcy flow, i.e., μf 1 ρ g − ∇(ε p) = V f − μ ′f ∇ 2 V f (2.444) K ε where the effective viscosity denoted by μ ′f = μ ′f ( μ f , ε , ℑ) is a function of the dynamic viscosity, the porosity, and the tortuosity ℑ of the porous media. The tortuosity is a measure of the connectivity of void space in a porous zone. Equation (2.444) is called the Brinkman equation, and it is usually noted that it is valid for high porosity, ε > 0.8 . When the porosity is low, the stresses felt on a fluid in one pore are communicated to the fluid in another pore mainly by pressure because the solid matrix prevents direct viscous interaction of the fluid in separate pores. Even though this equation is only valid for higher porosities, the Laplacian operator of V f is needed when the boundaries of the porous media affect the flow field. An empirical momentum equation can be heuristically obtained if Darcy’s law, eq. (2.439), is combined with the inertial drag component, eq. (2.442), and Brinkman’s equation, eq. (2.444). μf Cf 1 ρ g − ∇(ε p) = − μ ′f ∇ 2 V f + Vf + 1 2 ρ f Vf Vf (2.445) K ε K However, for a more complete understanding, this equation is compared to a volume-averaged momentum equation derived in Section 2.5.2. If the fluid Chapter 2 Ge era i ed G ver i g Equati s 171 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press phase is incompressible, i.e., ∂ ερ f V f ∂t = ∇ ⋅ ε τ ′f ρf Vf sm , f f f = ρ f , the volume-averaged momentum Vf f equation for the fluid phase, eq. (2.334), becomes ( f ( f ) + ερ g + F f ) + ∇ ⋅ (ερ f ) ˆˆ − ∇ ⋅  ερ f V f V f   Vf f    Vf f (2.446) ˆˆ + Vf Vf f where Vf Vf f in eq. (2.334) is replaced by f and Xf is replaced by g, i.e., gravity is the only body force. Since the flow in a porous media is usually laminar because of the characteristic diameter, the product of the velocity deviations, ˆˆ Vf Vf f , is also small and therefore neglected henceforth. The third term on the right hand side of eq. (2.446) represents interaction between the solid matrix and the fluid phase. For an incompressible fluid, the volume average of the stress tensor is τ ′f f = − pf f  I + μ f ∇ V f  f + ∇ Vf ( fT )  (2.447) Since the local static pressure is needed in the equation of state for an ideal gas, or for the temperature at a liquid vapor interface, the local pressure is more accurately described as the intrinsic phase average pressure, not the extrinsic phase averaged pressure. If the deviation of pressure within a fluid volume ˆ element is small, pk ≈ 0 , then the volume-averaged pressure is pf = ε pf f =εp. (2.448) Substituting eq. (2.447) into eq. (2.446) and considering eq. (2.448), the momentum equation becomes ∂ ερ f V f ∂t ( f μ f ∇ ε Vf 2 ( f ) + ∇ ⋅ (ερ V ) + ερ g + F f f f f Vf f ) = −∇ (ε p) + (2.449) sm , f Using the continuity equation, eq. (2.437), and assuming constant porosity, the momentum equation becomes ρf   1 ∂ Vf ε  ∂t + Vf ε2 ∇ ⋅ Vf ( )   1  (2.450) ∇ 2 V f + ρ f g + Fsm , f ε ε The empirical momentum relationships are heuristically related to the volumeaveraged momentum equation through reasonable observations. In the second term on the right hand side of eq. (2.450), the value of μ f / ε represents μ ′f in eq. (2.445). The fourth term on the right-hand side of eq. (2.449) is replaced by the Darcian relationship and the quadratic drag term. The Darcian drag term = −∇p + μf 172 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press comes from the deformation stress caused by the solid/fluid interaction, and the quadratic drag term comes from the pressure coefficient associated with the fluid flowing around the solid matrix. ρf  +  1 ∂ Vf ε  ∂t  ∇ ⋅ V f  = −∇p  ε  μf ρfCf + ρfg− Vf − 1 2 K K + Vf 2 ∇2 V f Vf Vf (2.451) ε This equation reduces to eq. (2.445) if the inertia terms are negligible. In eq. (2.451), the viscous interaction between the fluid and the solid matrix is replaced by Darcy loss and inertial loss terms. The shear losses between the fluid and the solid matrix are accounted for without having any information on the characteristics of the flow in each pore. Tab e 2 2 Properties of porous media materials (Nield and Bejan, 1999; Faghri 1995) Solid Matrix Black Slate Powder Brick Cigarette Cigarette filters Coal Concrete (ordinary mixes) Copper powder (hotcompacted) Cork board Fiberglass Hair (on mammals) Hair felt Leather Limestone (dolomite) Sand Screen, SST, 200 mesh Screen, Nickel, 50 mesh, sintered Felt, Sintered, SST Felt, Nickel, A30 Powder, Sintered, Nickel Beads, Monel, 70-80 Mesh Sandstone (“oil sand”) Silica grains Silica powder Soil Wire crimps Porosity 0.57 – 0.66 0.12 – 0.34 0.17 – 0.49 0.02 – 0.12 0.02 – 0.07 0.09 – 0.34 0.88 – 0.93 0.95 – 0.99 0.56 – 0.59 0.04 – 0.10 0.37 – 0.50 0.733 0.625 0.916 0.815 0.597 0.4 0.08 – 0.38 0.65 0.37 – 0.49 0.43 – 0.54 0.68 – 0.76 3.3×10-10 – 1.5×10-9 2.4×10-11 – 5.1×10-11 8.3×10-10 – 1.2×10-9 9.5×10-14 – 1.2×10-13 2×10-15 – 4.5×10-14 2.0×10-11 – 1.8×10-10 5.2×10-11 6.63×10 -10 μf Permeability K (m2) 4.9×10-14 – 1.2×10-13 4.8×10-15 – 2.2×10-13 1.1×10-9 Effective Pore Radius reff (m) ( m −1 ) ΔV 7×105 - 8.9×105 ΔAsm (5.6 – 7.7)×104 (1.2 – 1.6)×106 (1.5– 2.2)×104 58×10 305×10-6 94×10-6 120×10-6 58×10-6 96.9×10-6 (6.8 – 8.9) ×105 (2.9– 4.0)×103 -6 5.46×10-10 3.06×10-11 7.75×10-11 5×10-16 – 3×10-12 1.3×10-14 – 5.1×10-14 2.9×10-13 – 1.4×10-11 3.8×10-9– 1.0×10-8 Chapter 2 Ge era i ed G ver i g Equati s 173 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The permeability, K, can be obtained by empirical relations, which are well tabulated for various types of porous media. For packed beds of spherical particles of diameter d, the permeability is d 2ε 3 K= (2.452) 150(1 − ε ) 2 Further details on permeability and porosity of various porous materials, including porous wicks in heat pipes, can be found in Kaviany (1995) and Faghri (1995). For convenience, the porous zone properties for a variety of applications are presented in Table 2.2, including K, ε , ΔAsm / ΔV and reff. Example 2.13 A cylindrical tube has a porous zone inside it. The porous zone consists of N smaller cylindrical tubes within it as shown in Fig. 2.18. The large cylinder has a radius of R0, while the smaller internal cylinders all have a radius of r0. Derive the permeability from the momentum equation derived by the volume averaging technique. Assume that the flow is incompressible, a fully-developed flow within each tube, and that the volume-averaged velocity, Vk , is constant throughout the larger cylinder. Solution: There are two parts in solving this problem, (1) knowing the frictional losses in each pore, and (2) then substituting these losses into the simplified porous media momentum equation, eq. (2.439). Since the flow is steady, incompressible and one-dimensional, continuity is already specified by having a uniform flow through the cylinder. Also the volume averaged momentum equation, eq. (2.334), can be simplified by: ∂ ερ f V f ∂t ( f ) + ∇ ⋅ (ερ f Vf Vf f ) = ∇ ⋅ (ε τ ) − Δ1V  f ΔAsm ρ f V f ( V f − VI ) ⋅ ndA 0 0 + 1 ΔV 0  Ak τ ⋅ ndAk + ερ f X f 0 x, i R0 2r0 Figure 2.18 Cross-sectional and side view of a cylinder with cylindrical pores. 174 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press With constant properties, and the assumed fully developed profile, the above equation becomes ε ∂ pf ∂x f = 2μ f ΔV 2 0  ( D ⋅ n ) dA ΔAsm f (2.453) The fully developed velocity profile in a single tube is: Vf ⋅ i = (r − r 2  dp f  −  4 μ  dx  ) (2.454) where i is the unit vector in the x-direction. From the velocity profile, we can obtain the rate of deformation, D f ⋅ n , with respect to the smaller cylindrical tubes surface. D f ⋅n = 1 ∂V f ⋅ i r0 dp f = 2 ∂r 4μ dx (2.455) Integrating this value over all of the pore wall areas gives the total frictional loss. 2 2μ f N  r0  dp f (2.456) ( D f ⋅ n ) dA = N   ΔV  ΔAsm  R0  dx The extrinsic average velocity is equal to the area average velocity. Vf ⋅ i = 2N 2 R0  r0 0 V f ⋅ irdr = Nr04  dp f  −  2 8μ R0  dx  (2.457) V f , and Now the frictional loss can be converted to a function of substituted into the original momentum equation, eq. (2.453). ε d pf dx f = −8μ f V f ⋅ i r02 Nr02 2 R0 (2.458) The porosity is: ε= f (2.459) To maintain the form of Darcy’s law, − d pf dx = μ K Vf ⋅ i = Nr04 2 8 R0 2 8 R0 μ Nr04 Vf ⋅ i (2.460) Therefore, the permeability is K= (2.461) Example 2.14 When using Darcy’s law in a porous media, the effect of a no-slip boundary condition cannot be accounted for because of the order of the equation. The shear stress caused by a wall (Fig. 2.19) can be captured using Brinkman’s equation. Use the analysis of a semi-infinite, fully Chapter 2 Ge era i ed G ver i g Equati s 175 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Elemental Control Volume x, uf y, vf Figure 2.19 Porous region near a wall developed velocity in a porous zone in no gravitational field to determine the criteria when the wall effects should be incorporated in the governing equations. Develop a solution using a no-slip wall condition and use this solution to determine the significance of the wall for different porous properties. Also, find the volume average velocity profile of a flow through a porous media bounded by parallel plates under a given pressure gradient. Use no-slip boundary conditions, and compare the mean flow rate to the mean flow rate using Darcy’s law. Solution: Brinkman’s equation in a one-dimensional porous zone is: − 2  uf 1 ∂ uf dp = μf  − K ε ∂y 2 dx      (2.462) where uf The extrinsic phase-averaged velocity in the x-direction and the effective viscosity is assumed to be μ ′ = μ / ε . With boundary conditions at the wall and at the free stream of: y = 0, u f = 0 , u f = 0 (2.463) y = ∞, ∂ Vf ∂y =0 (2.464) Assume that the volume average velocity at the wall is zero (knowing it may not be correct). The solution to the semi-infinite Brinkman’s equation is: 176 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press − μf uf   ε  = 1 − exp  − y    dp K     K dx (2.465) which is plotted in Fig. 2.20. It can be seen in Fig. 2.20 that the wall effects are confined to a distance of approximately 6 K / ε . The average pore length in the ydirection is about K / ε ; therefore, the wall effects will only be noticeable in cells that are within approximately six pore lengths from the wall. For computational considerations, if the grid is not resolved within six average pore lengths, then there is no advantage to using Brinkman’s equation over Darcy’s law. To quantify the overall size of the solution to determine what scale is large enough such that wall effects can be neglected, we must solve Brinkman’s equation for a bounded flow between two parallel plates ( 0 ≤ y ≤ h ). − 2  uf 1 ∂ uf dp = μf  − K ε ∂y 2 dx      (2.466) The boundary conditions are: y =0, u =0 y=h, u =0 (2.467) (2.468) The general solution to this equation after applying boundary conditions in eq. (2.467) and (2.468) is: ε  K dp  ε  uf = − (2.469) 1 − A exp   K y  − B exp  − K y       μ f dx      − μf uf K dp dx y ε /K Figure 2 20 Solution to Brinkman’s equation in a semi-infinite domain. Chapter 2 Ge era i ed G ver i g Equati s 177 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press where the constants A and B are:    ε  ε    ε  A = exp  − h  − 1 exp  − h  − exp       K h   K   K           ε     ε  ε B = 1 − exp  h   exp  − h  − exp  K  K  K h              −1 (2.470) (2.471) uf −1 The average velocity found by integrating between the plates is uf = 1 h over the distance (2.472)  h 0 u f dy Therefore, the nondimensional average velocity is: 2 ( exp ( −η ) − 1) ( exp (η ) − 1) μu − = 1− dp exp ( −η ) − exp (η ) η K dx (2.473) where η is equal to the non-dimensional distance, h ε / K . Equation (2.473) is presented in Fig. 2.21. By analyzing Fig. 2.21, it can be noted that the wall effects are insignificant when the overall length of the system η is greater than 103, and therefore Darcy’s law is a good approximation. However, when η is less than 103, Brinkman’s equation should be applied. If η is on the order of 1, the volume averaging approach is not good, and the local velocities should be modeled. − μf uf K dp dx Figure 2.21 Non-dimensional average velocity vs. non-dimensional distance. 178 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press E ergy Equati The derivation of the energy equation in a porous medium is similar to that of continuity and momentum equations. The volume-averaged energy equation (2.341) is valid for both fluid phase and the solid matrix. Assuming the fluid phase is incompressible, and neglecting the viscous dissipation, the energy equation for the fluid phase becomes ∂ ερ f h f ∂t ( f ) + ∇ ⋅ (ερ f Vf f hf f ) = −∇ ⋅ q′′ f ′′′ + q′′′ + qsm , f f Vf f (2.474) where Vf h f f in eq. (2.341) is replaced by hf f in eq. (2.474), i.e., the product of the deviations were neglected. ′′′ qsm , f is heat transfer from solid matrix to the fluid phase. By applying the continuity equation (2.437), eq. (2.474) becomes ερ f ∂ hf ∂t f + ερ f V f ∂ hsm ∂t f ⋅∇ hf f ′′′ = −∇ ⋅ q′′ + q′′′ + qsm , f f f (2.475) The volume-averaged energy equation in the solid matrix is: ε sm ρ sm sm sm ′′′ ′′′ = −∇ ⋅ q′′ + qsm − qsm , f sm (2.476) where the convection term is dropped since the solid matrix is in the solid state and has no velocity. The thermal boundary conditions at the fluid-solid interface, ΔAsm, are: (2.477) T f = Tsm q′′ ⋅ n f = q′′ ⋅ n sm f sm (2.478) These boundary conditions between the fluid and the solid matrix are what make ′′′ up the heat transfer term between the fluid and solid, qsm , f . Equations (2.477) and (2.478) are very general, and no assumptions about the thermal equilibrium between the fluid and the solid matrix have been made. If the fluid and the solid matrix are considered to be in local thermal equilibrium, Tf f = Tsm sm = T , the energy equations for the solid and the fluid can be added to find the total porous media energy equation: ερ f ∂ hf f ∂t ∂t ′′′ = −∇ ⋅ q′′ − ∇ ⋅ q′′ + q′′′ + qsm f sm f + ε sm ρ sm sm ∂ hsm sm + ερ f V f f ⋅∇ h f f (2.479) The energy equation can be further simplified by assuming that the enthalpy is only a function of temperature and constant specific heat: ( ρ c p )eff ∂T ′′′ + ( ρ c p ) f V f ⋅ ∇T = ∇ ⋅ keff ∇T + qeff ∂t (2.480) Chapter 2 Ge era i ed G ver i g Equati s 179 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press The effective heat capacity, thermal conductivity and heat generation rates per unit volume, denoted by the subscript “eff,” are defined by (2.481) ( ρ c p ) = ε f ( ρ c p ) f + ε sm ( ρ c p )sm eff keff = ε keff , f + ε sm keff , sm ′′′ qeff = ε q′′′ f f (2.482) (2.483) ′′′ + ε sm qsm sm where keff , f and keff , sm are, respectively, the effective thermal conductivity of the fluid and solid-matrix in the porous media, and both of them depend on porosity and pore structure in the porous media. When the thermal conductivities of the fluid and the solid matrix are close to each other, k f  ksm , one can assume that keff , f ≈ k f and keff , sm ≈ ksm , and use eq. (2.482) to evaluate the effective thermal conductivity. In a case where thermal conductivities of the fluid and the solid matrix differ significantly, the following equation should be used (Hadley, 1986): keff ε f0 + (1 − ε f0 )ksm / k f = (1 − α 0 ) kf 1 − ε (1 − f 0 ) + ε (1 − f0 ) (2.484) 2(1 − ε )(ksm / k f ) 2 + (1 + 2ε )ksm / k f + α0 (2 + ε )ksm / k f + 1 − ε where f0 = 0.8 + 0.1ε (2.485) −4.898ε 0 ≤ ε ≤ 0.0827   (2.486) log α 0 = −0.405 − 3.154(ε − 0.0827) 0.0827 ≤ ε ≤ 0.298  −1.084 − 6.778(ε − 0.298) 0.298 ≤ ε ≤ 0.580  The above energy equation, eq. (2.480), is only valid when the fluid and the solid matrix are assumed to be in local thermal equilibrium. While local thermal equilibrium is an often used hypothesis for heat transfer in porous media, it is well recognized that local equilibrium between different phases in a system of solid and fluid cannot be achieved when thermal properties for different phases differ widely, or during rapid heating or cooling. The fluid and solid-matrix energy equations need to be kept separate, and this situation is referred to as local thermal non-equilibrium. In this situation, heat transfer from the solid matrix to ′′′ the fluid phase, qsm, f , in eqs. (2.475) and (2.476) can be modeled as a convective boundary condition, with a heat transfer coefficient of hc . When Fourier’s law of heat conduction is applied, the energy equations for the fluid and solid matrix are  ∂T f  + V f ⋅∇T f  (ρcp ) f  ε  ∂t  = ∇ ⋅ (ε keff , f ∇T f ) + ε q′′′ f f + h f , sm ΔAsm Tsm − T f ΔV ( ) (2.487) 180 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ε sm ( ρ c p ) sm ∂Tsm ΔAsm ′′′ Tsm − T f (2.488) = ∇ ⋅ (ε sm keff ,sm ∇Tsm ) + ε sm qsm − hc ∂t ΔV ( ) Equations (2.487) and (2.488) are more accurate representations of the thermal field in a porous zone. However, they require an empirical relationship for the heat transfer coefficient between the solid matrix and the fluid phase. Species The volume-averaged species eq. (2.349) can be readily applied to the fluid phase in the porous media: ∂ ε ρ f ,i ∂t ( f Applying the continuity equation (2.437) and the definition of the mass fraction, eq. (1.77), the averaged-species equation becomes ερ f ∂ ω f ,i ∂t f ) +∇⋅ ε ρ ( f f ,i Vf f ) = −∇ ⋅ J f ,i f + m′′′,i (2.489) + ρ f v f ⋅∇ ω f ,i f f = −∇ ⋅ J f ,i + m ′′′,i (2.490) 2 5 Fu da e ta s f Turbu e ce 2 5 1 Descripti f Turbu e ce Turbulent flow is the most commonly encountered type of viscous flow but its theoretical treatment is not as well developed as that for the laminar flow. In this section, turbulent flow will be examined with particular attention paid to mechanism of the momentum and energy transfer. Hinze (1975) was the first one who presented the structure of the turbulent flow in detail. The constituencies of a turbulent flow are eddies and vortices at different sizes. Although the effect of the molecular viscosity are important within these eddies and vortices, the interactions of these eddies and vortices are dominant on the transport phenomena in the turbulent flow. For turbulent air flow at a speed of 100 m/s, the typical size of eddies are on the order of 1 mm; this is much larger than the mean free path of the air molecules, which is on the order of 66 nm standard condition (1 atm pressure and 25 °C). In a cube with 1 mm on a side, there are 2.46 × 1016 air molecules (see Example 1.1 in page 15). Therefore, eddies in the turbulent flow are considered to be “large” from the microscale point of view. On the other hand, they are too small from the macroscopic point of view. When the velocity of the turbulent air flow is measured using a hot wire anemometer, the local instantaneous velocity varies about 10% of its average velocity due to the Chapter 2 Ge era i ed G ver i g Equati s 181 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u u(x,y,z,t) u u(x,y,z,t) u ( x, y , z ) time u ( x, y , z , t ) time (b) unsteady mean flow (a) steady mean flow Figure 2 22 Time dependence of velocity component in the x-direction in a turbulent flow random eddies and vortices. In other words, for an air flow with mean velocity of 100 m/s, a velocity fluctuation on the order of 10 m/s can be expected. Compared with the average random molecular velocity of 465 m/s, this fluctuation is relatively small. The frequency of the velocity fluctuation is about 1 to 104 Hz, whereas the frequency of molecular collision is about 7 × 109 Hz. It appears that the interactions between eddies in a turbulent flow are very different from the molecular collision which is the mechanism of the viscous effect for laminar flow. Since the number of eddies in a turbulent flow is enormous, it is very difficult to model or measure the flow structure within each eddy. Therefore, time averaging technique is usually used to obtain the governing equations that are manageable. In turbulent flow, the transport phenomena variables (i.e., u, v, w, T, p, etc.) always vary with time. For example, the magnitude and direction of the instantaneous velocity vector are different from those of the time-averaged velocity. While the instantaneous velocity in a turbulent flow is always timedependent, the time averaged velocity can be either time-independent or time dependent [for the velocity component in the x-direction, see Fig. 2.22 (a) and (b)]. For any given location (x, y, z) and time t, the local instantaneous velocity can be expressed as a summation of its mean value and its fluctuation, i.e., u = u ( x , y, z , t ) + u′( x , y, z , t ) (2.491) where u ( x , y, z , t ) = 1 Δt  t +Δt t u( x , y, z , t )dt (2.492) is the time-averaged velocity at point (x, y, z). The time interval for the timeaverage, Δt , must be very long compared with the duration of fluctuation. The mean value of the fluctuation must be zero, i.e., u′ = 1 Δt  t +Δt t u′( x , y, z , t )dt = 0 (2.493) Similarly, the velocity components in the y- and z- directions can be expressed as v = v ( x , y, z , t ) + v′( x , y, z , t ) (2.494) 182 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press (2.495) In addition, the local instantaneous pressure, temperature and mass fraction can also be expressed as sums of their mean values and fluctuations: p = p( x , y, z , t ) + p′( x , y, z , t ) (2.496) T = T ( x , y, z , t ) + T ′( x , y, z , t ) (2.497) ′( x , y, z , t ) ω = ω ( x , y, z , t ) + ω (2.498) Similar to eq. (2.493), the volume averages of the fluctuations for these variables are zero, i.e., v′ = w ′ = p ′ = T ′ = ω ′ = 0 (2.499) For any variables φ and ψ , we have the following properties about the averaging of their fluctuations (Oosthuizen and Naylor, 1999): ∂φ ∂φ ∂φ ′ ∂φ ′ ∂n = ∂n w = w( x , y, z , t ) + w′( x , y, z , t ) , ∂n = ∂n =0 φψ = φψ , φφ ′ = φ φ ′ = 0, φ + ψ = φ + ψ (2.500) where n can be any of the special variable (x, y, or z). Although the mean values of the fluctuations are zero, these mean values do contribute to the mean value of some physical quantities. The kinetic energy per unit volume is 1 ke = [(u + u′) 2 + (v + v′) 2 + (w + w′) 2 ] 2 1 = [(u 2 + 2uu′ + u′2 ) + (v 2 + 2vv′ + v′2 ) + (w 2 + 2ww′ + w′2 )] 2 Since uu′ = uu′ = 0, vv′ = vv′ = 0, and ww′ = ww′ = 0 , the above expression becomes 1 ke = [(u 2 + v 2 + w 2 + u′2 + v′2 + w′2 )] 2 (2.501) which indicates that the time average of the squares of the velocity fluctuations are not zero. Similarly, the products of the fluctuations of different variables are generally not zero: u′T ′ ≠ 0, u′v′T ′ ≠ 0 , etc. Equation (2.501) shows that both mean velocity and velocity fluctuation contribute to the total kinetic energy in a turbulent flow. The level or intensity of turbulence is defined as (Welty et al., 2000) ke = ( u ′ 2 + v′ 2 + w ′2 ) / 3 U∞ (2.502) where the numerator is the root mean square value of fluctuation and the denominator, U∞, is the mean velocity of the flow. The intensity of turbulence will have very important effect on transition and separation of boundary layer, and the rate of heat and mass transfer. The intensity of turbulence is another quantity in addition to Reynolds number that we will need to describe the turbulent flow. Chapter 2 Ge era i ed G ver i g Equati s 183 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2 5 2 Ti e Averaged G ver i g Equati s The transport phenomena equations including the Navier-Stokes and the energy equation are exact for solving the turbulent problem if their complete unsteady forms are solved. In general, this approach is not taken because of complexity and significant computer time requirement. In addition, the turbulent problems, in their gross aspect, are assumed steady state as noted in the previous subsection. Time averaging of transport phenomena equations should provide the net effect of the turbulent perturbation. For an incompressible binary system with constant properties, the continuity, Navier-Stokes, energy, and conservation of mass species equations in a Cartesian coordinate system are ∂u ∂v ∂w + + =0 ∂x ∂y ∂z  ∂2u ∂2u ∂2u  ∂u ∂u ∂u ∂u 1 ∂p +u +v +w =− +ν  2 + 2 + 2   ∂x ∂t ∂x ∂y ∂z ρ ∂x ∂y ∂z     ∂2v ∂2v ∂2v  ∂v ∂v ∂v ∂v 1 ∂p +u +v +w =− +ν  2 + 2 + 2   ∂x ∂t ∂x ∂y ∂z ρ ∂y ∂y ∂z     ∂2w ∂2w ∂2w  ∂w ∂w ∂w ∂w 1 ∂p +u +v +w =− +ν  2 + 2 + 2   ∂x ∂t ∂x ∂y ∂z ρ ∂z ∂y ∂z     ∂ 2T ∂ 2T ∂ 2T  ∂T ∂T ∂T ∂T +u +v +w =α 2 + 2 + 2   ∂x ∂t ∂x ∂y ∂z ∂y ∂z     ∂ 2ω ∂ 2ω ∂ 2ω  ∂ω ∂ω ∂ω ∂ω +u +v +w = D 2 + 2 + 2   ∂x ∂t ∂x ∂y ∂z ∂y ∂z    (2.503) (2.504) (2.505) (2.506) (2.507) (2.508) Since we made no assumption about the nature of the flow in the above equations, the local instantaneous parameters in a turbulent flow satisfy eqs. (2.503) – (2.508). However, applying the above equations to turbulent flow will require very fine grid size to be able to capture the detailed flow structure within each eddy – which is very challenging and, in most cases, unnecessary. For most practical engineers, the mean values of various parameters in a turbulent flow are of interest. Therefore, we will obtain the governing equations in term of timeaveraged parameters. For turbulent flow, the velocity components, pressure, temperature and mass fraction of species can be expressed as sums of their mean values and fluctuations as it was indicated in eqs. (2.491), and (2.494) – (2.498). Substituting eqs. (2.491), (2.494), and (2.495) into the continuity equation (2.503), one obtains: ∂ ∂ ∂ (u + u′) + (v + v′) + (w + w′) = 0 ∂x ∂y ∂z (2.509) By taking time averaging of eq. (2.509), we have 184 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press ∂ ∂ ∂ (u + u′) + (v + v′) + (w + w′) = 0 ∂x ∂y ∂z ∂u ∂v ∂w + + =0 ∂x ∂y ∂z (2.510) Considering eq. (2.500), the continuity equation (2.510) becomes (2.511) which is the continuity equation in terms of time averaged velocity. It can be seen that it has the exactly same form as that in terms of local instantaneous velocity, eq. (2.503). In other words, the contributions of the fluctuations of velocity components are averaged out during the time averaging. Substituting eqs. (2.491), and (2.494) – (2.496) into the momentum equation (2.504) in the x-direction and taking time average, one obtains ∂u ∂u′ ∂u ∂u′ ∂u ∂u′ ∂u ∂u′ + +u +u + u′ + u′ +v +v ∂t ∂t ∂x ∂x ∂x ∂x ∂y ∂y + v′ ∂u ∂u′ ∂u ∂u′ ∂u ∂u′ + v′ +w +w + w′ + w′ ∂y ∂y ∂z ∂z ∂z ∂z =−  ∂ 2 u ∂ 2 u′ ∂ 2 u ∂ 2 u′ ∂ 2 u ∂ 2 u′  1 ∂p 1 ∂p′ − +ν  2 + 2 + 2 + 2 + 2 + 2   ∂x ρ ∂x ρ ∂x ∂x ∂y ∂y ∂z ∂z    (2.512) Considering eq. (2.500), the momentum equation in the x-direction becomes: ∂u ∂u ∂u ∂u +u +v +w ∂t ∂x ∂y ∂z =−  ∂ 2 u ∂ 2 u ∂ 2 u   ∂u′ ∂u′ ∂u′  1 ∂p ′ ′ + ν  2 + 2 + 2  −  u′   ∂x   ∂x + v ∂y + w ∂z  ρ ∂x ∂y ∂z     (2.513) where fluctuation components only appear in the last parenthesis on the righthand side. In order to simplify the terms that contain fluctuation component, the continuity equation (2.503) is multiplied by u, i.e., u ∂u ∂v ∂w +u +u =0 ∂x ∂y ∂z Substituting eqs. (2.491), (2.494), and (2.495) into the above equation and performing time averaging to the resultant equation, we have: u +u′ ∂u ∂u′ ∂u ∂u′ ∂v ∂v′ +u + u′ + u′ +u +u ∂x ∂x ∂x ∂x ∂y ∂y ∂v ∂v′ ∂w ∂w′ ∂w ∂w′ + u′ +u +u + u′ + u′ =0 ∂y ∂y ∂z ∂z ∂z ∂z which can be simplified to u ∂u ∂v ∂w ∂u′ ∂v′ ∂w ′ +u +u + u′ + u′ + u′ =0 ∂x ∂y ∂z ∂x ∂y ∂z Considering the continuity equation in terms of mean velocity component, eq. (2.511), the above equation is further simplified to: Chapter 2 Ge era i ed G ver i g Equati s 185 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press u′ ∂u′ ∂v′ ∂w′ + u′ + u′ =0 ∂x ∂y ∂z (2.514) Substituting the eq. (2.514) to eq. (2.513), the momentum equation in the xdirection becomes ∂u ∂u ∂u ∂u +u +v +w ∂t ∂x ∂y ∂z =−  ∂ 2 u ∂ 2 u ∂ 2 u   ∂ u′2 ∂ v′u′ ∂ w′u′  1 ∂p  +ν  2 + 2 + 2  −  + +   ρ ∂x ∂y ∂z  ∂y ∂z   ∂x  ∂x   (2.515) By following the same procedure, the momentum equations in the y- and zdirections become: ∂v ∂v ∂v ∂v +u +v +w ∂t ∂x ∂y ∂z =−  ∂ 2 v ∂ 2 v ∂ 2 v   ∂ u′v′ ∂ v′2 ∂ w ′v′  1 ∂p  +ν  2 + 2 + 2  −  + ++  ∂x ∂x ∂z  ρ ∂y ∂y ∂z   ∂y    (2.516) ∂w ∂w ∂w ∂w +u +v +w ∂t ∂x ∂y ∂z =−  ∂ 2 w ∂ 2 w ∂ 2 w   ∂ u′w′ ∂ v′w ′ ∂ w′2  1 ∂p  +ν  2 + 2 + 2  −  + +   ∂y ∂z  ρ ∂x ∂y ∂z   ∂x  ∂x   (2.517) It can be seen from eqs. (2.515) – (2.517) that extra terms associated with the fluctuations of the velocity components appeared in the momentum equations. These extra terms account for the momentum transfer caused by the velocity fluctuations and are termed as turbulent or Reynolds stress terms. The difference between the viscous stress and Reynolds stress is that the former occurs at the molecular scale whereas the latter is caused by the interaction between eddies. To obtain the energy equation in terms of mean values of various variables, eqs. (2.491), and (2.494) – (2.497) are substituted into the energy equation (2.507) and time averaging is performed to the resultant equations to obtain ∂T ∂T ′ ∂T ∂T ′ ∂T ∂T ′ ∂T ∂T ′ + +u +u + u′ + u′ +v +v ∂t ∂t ∂x ∂x ∂x ∂x ∂y ∂y + v′ ∂T ∂T ′ ∂T ∂T ′ ∂T ∂T ′ + v′ +w +w + w′ + w′ ∂y ∂y ∂z ∂z ∂z ∂z  ∂ 2T ∂ 2T ′ ∂ 2T ∂ 2T ′ ∂ 2T ∂ 2T ′   =α 2 + + + + +  ∂x ∂x 2 ∂y 2 ∂y 2 ∂z 2 ∂z 2    (2.518) Considering eq. (2.500), the energy equation becomes: ∂T ∂T ∂T ∂T +u +v +w ∂t ∂x ∂y ∂z 186 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press  ∂ 2T ∂ 2T ∂ 2T =α 2 + 2 + 2  ∂x ∂y ∂z    ∂T ′ ∂T ′ ∂T ′  ′ ′  −  u′    ∂x + v ∂y + w ∂z    (2.519) By following the similar procedure that eq. (2.514) was derived, the following equation can be obtained (see Problem 2.60): T′ ∂u′ ∂v′ ∂w′ +T′ +T′ =0 ∂x ∂y ∂z (2.520) Substituting the above equation into eq. (2.519), the energy equation becomes ∂T ∂T ∂T ∂T +u +v +w ∂t ∂x ∂y ∂z  ∂ 2T ∂ 2T ∂ 2T =α 2 + 2 + 2  ∂y ∂z  ∂x   ∂ u′T ′ ∂ v′T ′ ∂ w′T ′  + + −    ∂y ∂z    ∂x (2.521) By following the same procedure, the conservation of species mass can be expressed as (see Problem 2.61) ∂ω ∂ω ∂ω ∂ω ∂t +u ∂x +v ∂y +w 2 ∂z  ∂ ω ∂ ω ∂ ω   ∂ u′ω ′ ∂ v′ω ′ ∂ w′ω ′  = D 2 + 2 + 2  − + +    ∂x  ∂y ∂z  ∂y ∂z   ∂x  2 2 (2.522) Therefore, transport phenomena in turbulent flow of a binary system with constant properties can be described by continuity equation (2.511), momentum equations (2.515) – (2.517), energy equation (2.521), and conservation of mass species equation (2.522). There are six time-averaged unknowns, u , v , w, p, T , and ω and 12 products of their fluctuations in these six equations. To use these equations to solve turbulent transport phenomena, an appropriate turbulent model must be employed to express the partial derivatives of the time-averaged product of the fluctuations in terms of mean values of various variables. While many turbulent models are almost entirely empirical, some models based on additional partial differential equations can be used to describe these terms. Although the empiricisms cannot be completely eliminated, these extra equations will allow us to introduce the empiricisms in a logical manner. We will show below the derivation of the turbulence kinetic energy equation, which is one of the most widely used additional equations. Substituting eqs. (2.491), and (2.494) – (2.496) into eq. (2.504) and multiplying u′ to the resultant equation yield: u′ ∂ ∂ (u + u′) (u + u′) + u′(u + u′) ∂t ∂x ∂ (u + u′) ∂ (u + u′) +u′(v + v′) + u′(w + w′) ∂y ∂z =−  ∂ 2 (u + u′) ∂ 2 (u + u′) ∂ 2 (u + u′)  u′ ∂ + + ( p + p′) + ν u′   2 ρ ∂x ∂y 2 ∂z 2  ∂x  Chapter 2 Ge era i ed G ver i g Equati s 187 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press which can be time averaged to get u′ ∂u′ ∂u′ ∂u ∂u′ ∂u′ ∂u ∂u′ ∂u′ + uu′ + u ′2 + u ′2 + vu′ + u′v′ + u′v′ + wu′ ∂t ∂x ∂x ∂x ∂y ∂y ∂y ∂z  ∂ 2 u′ ∂ 2 u ′ ∂ 2 u′  ∂u ∂u′ 1 ∂p′ + u′w′ = − u′ + ν u′  2 + 2 + 2  ρ ∂x ∂z ∂z ∂y ∂z   ∂x +u′w′ The above equation can be rearranged to ∂u ∂u ∂u 1 ∂ u ′2 1 ∂ u ′2 1 ∂ u ′2 1 ∂ u ′2 +u +v +w + u ′2 + u′v′ + u′w′ ∂x ∂y ∂z 2 ∂t 2 ∂x 2 ∂y 2 ∂z  ∂u′ ∂u′ ∂u′  1 ∂p′ +  u ′2 + u′v′ + u′w′ + ν u′∇ 2 u′  = − u′   ∂x ∂y ∂z  ρ ∂x  (2.523) In order to further simplify eq. (2.523), the last term on the left hand side can be rewritten as  2 ∂u′ ∂u′ ∂u′  + u′v′ + u′w′  u′   ∂x ∂y ∂z     ∂ u′3 ∂ u′2 v′ ∂ u′2 w′  1  ∂u′ ∂v′ ∂w′  1  − u ′2  = + + + +  2  ∂x ∂y ∂z  2  ∂x ∂y ∂z    It can be shown that the second term in the above equation is zero by multiplying u′2 to the continuity equation (2.503) and time averaging (See Problem 2.62). Thus, eq. (2.523) becomes 1 ∂ u ′2 1 ∂ u ′2 1 ∂ u ′2 1 ∂ u ′2 ∂u ∂u ∂u +u +v +w = −u ′2 − u′v′ − u′w′ 2 ∂t 2 ∂x 2 ∂y 2 ∂z ∂x ∂y ∂z ∂p′ 1  ∂ u′3 ∂ u′2 v′ ∂ u′2 w ′   + ν u′∇ 2 u′ − + + (2.524) ∂y ∂z  ρ ∂x 2  ∂x   By multiplying v′ to the momentum equation (2.505) and performing time − 1 u′ averaging, the following equation can be obtained: ∂v ∂v ∂v 1 ∂ v′ 2 1 ∂ v′2 1 ∂ v′2 1 ∂ v′ 2 +u +v +w = −u′v′ − v ′2 − v′w′ ∂x ∂y ∂z 2 ∂t 2 ∂x 2 ∂y 2 ∂z − 1 ρ v′ ∂p′ 1  ∂ v′2 u′ ∂ v′3 ∂ v′2 w′   + ν v′∇ 2 v′ − + + ∂y 2  ∂x ∂y ∂z    (2.525) Similarly, the following equation can be obtained by multiplying w′ to the momentum equation (2.506) and performing time averaging: 1 ∂ w ′2 1 ∂ w ′2 1 ∂ w ′2 1 ∂ w ′2 ∂w ∂w ∂w +u +v +w = −u′w′ − v′w′ − w ′2 2 ∂t 2 ∂x 2 ∂y 2 ∂z ∂x ∂y ∂z − 1 ρ w′ ∂p′ 1  ∂ w′2 u′ ∂ w′2 v′ ∂ w′3   + ν w′∇ 2 w′ − + + ∂z 2  ∂x ∂y ∂z    (2.526) 188 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Adding eqs. (2.524) to (2.526) and defining the kinetic energy associated with velocity fluctuation, K= 12 ( u ′ + v ′ 2 + w ′2 ) 2 (2.527) the following additional equation of K is obtained: ∂K ∂K ∂K ∂K +u +v +w ∂t ∂x ∂y ∂z  ∂u ∂u ∂u ∂v ∂v ∂v = −  u ′2 + u′v′ + u′w′ + u′v′ + v′2 + v′w′ ∂x ∂y ∂z ∂x ∂y ∂z  + u′w′ ∂w ∂w ∂w  1 + v′w′ + w ′2 − ∂x ∂y ∂z  ρ  ∂p′ ∂p′ ∂p′  + v′ + w′  u′  ∂x ∂y ∂z   1  ∂ u′3 ∂ u′2 v′ ∂ u′2 w′ ∂ v′2 u′ ∂ v′3 ∂ v′2 w′ − + + + + + ∂y ∂z ∂x ∂y ∂z 2  ∂x  ∂ w′2 u′ ∂ w′2 v′ ∂ w′3   + ν (u′∇ 2 u′ + v′∇ 2 v′ + w′∇ 2 w′) + + + ∂x ∂y ∂z   (2.528) In eq. (2.528), the first term on the right hand side represents the rate of generation of turbulent kinetic energy by the interaction of the Reynolds stress with the mean shear. In other words, this term represents the shear production due to turbulence. This lowers the kinetic energy of the mean field. The last three terms on the right-hand side represent the spatial transport of turbulent kinetic energy. The first two terms represent the transport by turbulence, while the third term represents viscous transport. Additional terms that need to be considered that are not included in this equation include buoyant production and viscous dissipation terms. Buoyant generation of turbulent kinetic energy lowers the potential energy of the mean field. On the contrary, buoyant destruction of turbulent energy increases the potential of the mean field. This is in contrast to the shear generation of turbulent kinetic energy, which lowers the kinetic energy of the mean field. In order to use eq. (2.528) to determine the turbulence terms in the momentum equations (2.515) – (2.517), additional approximate relations will be needed to relate some terms to the time averaged quantities. While the above governing equations are in terms of time averaged variables, numerical methods that do not employ any averaging and turbulent models – known as Direct Numerical Simulation (DNS) have been developed recently (Moin and Krishnan, 1998; Yu et al., 2005). The DNS does not depend on any turbulent model (often empirical) and are based on the local instantaneous governing equations (2.503) – (2.508). The DNS is the closest to the real transport phenomena, but it is computationally most intensive. Since DNS is not practical for many engineering problems, alternative approaches like Reynolds Averaged Numerical Solution (RANS) or Large Eddy Simulation (LES) can be employed. For academic research, on the other hand, DNS can reproduce Chapter 2 Ge era i ed G ver i g Equati s 189 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press experimental resolution for experiments that are either very difficult or impossible. With advancements on parallel computing, one can break up the computational task and set to multiple CPUs to perform computation simultaneously. The scale of the number of grid points required for a DNS to capture the physics of the simplest turbulent flow is (2.529) η where L is the length scale of the largest turbulent eddy and η is the length scale of the viscous effects. The length scale of viscous effect – also referred to as Kolmogorov scale – is given as N ν 3  η =  ε   1/4 L (2.530) where ν is kinematic viscosity and ε is the rate of kinetic energy dissipation. It has been shown that the scale of L / η is N  L / η  Re3/4 (2.531) L where Re L is the turbulent Reynolds number. When DNS is employed to solve a three-dimensional turbulent problem, the number of grid required will be N 3  Re9/4 (2.532) L which shows that the DNS can very easily run into computational limitation of speed and memory for large Reynolds numbers. Nevertheless, significant progresses have been made in DNS of turbulent flow especially turbulent boundary layer structures (Moin and Krishnan, 1998). Refere ces Alazmi, B. and Vafai, K., 2000, “Analysis of Variants within the Porous Media Transport Models,” ASME Journal of Heat Transfer, Vol. 122, pp. 303-326. Alazmi, B. and Vafai, K., 2001, “Analysis of Fluid Flow and Heat Transfer Interfacial Conditions between a Porous Medium and a Fluid Layer,” International Journal of Heat and Mass Transfer, Vol. 44, pp. 1735-1749. Avedisian, C.T., 1997, “Soot Formation in Spherically Symmetric Droplet Combustion,” Physical and Chemical Aspects of Combustion, edited by Irvin Glassman, I., Dryer, F.L., and Sawyer, R. F., pp. 135-160, Gordon and Breach Publishers. Avedisian, C.T., 2000, “Recent Advances in Soot Formation from Spherical Droplet Flames at Atmospheric Pressure,” Journal of Propulsation and Power, Vol. 16, pp. 628-656. Bejan, A., 2004, Convection Heat Transfer, 3rd ed., John Wiley & Sons, New York. 190 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Boysan, F., 1990, A Two-Fluid Model for Fluent, Flow Simulation Consultants, Ltd., Sheffield, England. Cao, Y., and Faghri, A., 1993, “Simulation of the Early Startup Period of High Temperature Heat Pipes From the Frozen State by a Rarefied Vapor SelfDiffusion Model,” ASME Journal of Heat Transfer, Vol. 115, pp. 239-246. Curtiss, C.F., and Bird, R.B., 1999, “Multicomponent Diffusion,” Industrial and Engineering Chemistry Research, Vol. 38, pp. 2115-2522. Curtiss, C.F., and Bird, R.B., 2001, “Errata,” Industrial and Engineering Chemistry Research, Vol. 40, p. 1791. Edwards, D.K., Denny, V.E., and Mills, A.F., 1979, Transfer Process, Hemisphere, New York. Faghri, A., 1995, Heat Pipe Science and Technology, Taylor & Francis, Bristol, PA. Faghri, A., and Zhang, Y., 2006, Transport Phenomena in Multiphase Systems, Elsevier, Burlington, MA. Hadley, G.R., 1986, “Thermal Conductivity of Packed Metal Powders,” International Journal of Heat and Mass Transfer, Vol. 29, pp. 909-202. Hinze, J. O. 1975, Turbulence. 2nd ed., McGraw Hill, New York: Hirschfelder, J.O., Curtiss, C.F., and Bird, R.B., 1966, Molecular Theory of Gases and Liquids, Wiley, New York. Incropera, F.P., and DeWitt, D.P., 2001, Fundamentals of Heat and Mass Transfer, 5th ed., John Wiley & Sons, New York. Kaviany, M., 1995, Principles of Heat Transfer in Porous Media, 2nd ed., Springer Verlag, New York. Kays, W.M., Crawford, M.E., and Weigand, B., 2004, Convective Heat Transfer, 4th ed., McGraw-Hill, New York, NY. Kleijn, C.R., 1991, “A Mathematical Model of the Hydrodynamics and Gas Phase Reaction in Silicon LPCVD in a Single Wafer Reactor,” Journal of Electrochemical Society, Vol. 138, pp. 2190-2200. Lide, D.R. ed., 2009, CRC Handbook of Chemistry and Physics, 89th Ed., CRC Press, Boca Raton, FL. Lock, G.S.H., 1994, Latent Heat Transfer, Oxford Science Publications, Oxford University, Oxford, UK. Mahajan, R.L., 1996, “Transport Phenomena in Chemical Vapor-Deposition Systems,” Advances in Heat Transfer, Academic Press, San Diego, CA. Moin, P., and Krishnan, M., 1998, “Direct Numerical Simulation: A Tool in Turbulence Research,” Annu. Rev. Fluid Mech., Vol. 30, pp. 539-578. Chapter 2 Ge era i ed G ver i g Equati s 191 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Nield, D.A., and Bejan, A., 1999, Convection in Porous Media, 2nd ed., SpringerVerlag, New York. Oosthuizen, P.H., and Naylor, D., 1999, Introduction to Convective Heat Transfer Analysis, WCB/McGraw-Hill, New York. Ward, J. C., 1964, “Turbulent Flow in Porous Media,” ASCE J. Hyd. Div., Vol. 90, HY5, pp. 1-12. Welty, J.R., 1978, Engineering Heat Transfer, John Wiley & Sons, New York. Welty, J.R., Wicks, C.E., Wilson, R. E., 2000, Fundamentals of Momentum Heat and Mass Transfer, 4th ed., John Wiley & Sons, New York, NY. White, F.M., 1991, Viscous Fluid Flow, 2nd ed., McGraw-Hill, New York. Yu, H., Girimaji, S.S., Luo, L., 2005, “DNS or LES of Decaying Isotopic Turbulence with and without Frame Rotation Using Lattice Boltzmann Method,” Journal of Computational Physics, Vol. 209, pp. 599-616. Pr b e s 2.1. The local instance continuity equation can also be obtained by performing a mass balance for the differential control volume shown in Fig. P2.1. Derive the continuity equation and compare your result with eq. (2.54). Derive the local instance continuity equation in a cylindrical coordinate system using the control volume shown in Fig. P2.2. 2.2. Figure P2.1 Figure P2.2 2.3. 2.4. The continuity equation for incompressible flow is eq. (2.55). Is it valid for transient flow? Why or why not? In order to take water without stopping the train, a narrow trough of several thousand feet long can be placed in the midway between the rails of a railroad track. As a locomotive passes over the narrow trough, the water can be forced up into the tender through a scoop by the speed of the 192 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press locomotive (see Fig. P2.3). Assuming the cross-sectional area of the scoop is A, find the force acting on the train by the water using momentum balance in (a) a coordinate system that is attached to and moves with the locomotive, and (b) a fixed coordinate system. Locomotive tender Train Speed, U trough Figure P2.3 2.5. Use the integral form of the energy equation to calculate the steady mass flow rate of water in a pump that delivers 4 horsepower to the fluid, as given in the figure below, for the control volume shown. The pump inlet and outlet diameters are 10 inches and 4 inches, respectively. Control Volume 2 10” A1 u1 P1 1 pump A2 u2 P2 4” 3” Figure P2.4 2.6. Consider a liquid heated in an upward flow configuration in a vertical 10 meter long tube of constant diameter. The inlet average velocity and pressure are 2 m/sec and 350 kPa, respectively. Determine the heat needed to be added to the fluid so that the losses in internal energy equal 180 kJ/kg. Assume the density is 1000 kg/m3 and outlet pressure is 1atm. Water flows through a 3 inch diameter pipe at a rate of 45 gal/min. Frictional forces produce a pressure drop of 8 psi. Determine the 2.7. Chapter 2 Ge era i ed G ver i g Equati s 193 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press temperature change of the water if the heat transfer to the pipe can be neglected. 1 3 u1 A1 ρ1 T1 θ u3 A3 ρ3 T3 u2 A2 ρ2 T2 2 Figure P2.5 2.8. Determine the change in water temperature before (1) and (3) in terms of the inlet parameters A1, A2, u1, u2, cv and θ for the configuration shown in Fig. P2.5. Assume T1 = T2 and p1 = p2. The change in internal energy can be approximated as de=cvdT. A gas enters an insulated tank with total volume V through a pipe with diameter D with an inlet velocity and temperature of ui and Ti respectively as shown in Fig. P2.6. Using the integral form of the energy equation determine a first-order equation for the rate of change in temperature in the tank as a function of (V, D, ui, Ti). 2.9. Volume, V T Gas u i i D Figure P2.6 2.10. Determine viscous dissipation for compressible flow in the Cartesian coordinate system. 2.11. Write down the governing equations for two-dimensional, compressible flow with the effect of viscosity being accounted for. Assume that the 194 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press constant properties assumptions apply and that the only body force is gravity. 2.12. Fourier’s law of conduction [see eq. (1.66)] assumes that heat flux responds immediately to a temperature gradient. However, for cases with a very high temperature gradient or extremely short duration, heat is propagated at a finite speed, which can be accounted for by a modified heat flux model: q′′ + τ ∂q′′ = − k ∇T ∂t where τ is the thermal relaxation time. Derive the energy equation for a pure substance using the above model and eq. (2.92). Internal heat generation and viscous dissipation can be neglected. 2.13. The enthalpy for the ith component is hi is function of temperature, T, and pressure, p. Show that the substantial derivative in terms of the velocity of the ith component, Vi, can be obtained from eq. (2.171). u Inlet Boundary layer region y ,v D δ 0 x, u Hydrodynamic entrance length Figure P2.7 Fully developed flow 2.14. Formulate the hydrodynamics for developing flow in a duct formed by two infinite flat plates. Assume the flow to be laminar, incompressible, 2D and steady. It can also be assumed that the inlet velocity at the entrance of the duct is uniform, as shown in Fig. P2.7. Using the integral method, find out the hydrodynamic entrance and the skin friction coefficient in terms of Reynold’s number. Is the analysis valid for a circular pipe? 2.15. Formulate the hydrodynamics when flow is fully developed in a duct formed by two infinite flat plates and a round pipe. Assume the flow to be laminar, incompressible, two-dimensional and steady. Because flow is fully developed, it can be assumed that the velocity in the duct does not change in the x-direction. Determine the skin friction coefficient. Clearly describe the methods used to solve the momentum and continuity equations. Chapter 2 Ge era i ed G ver i g Equati s 195 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2.16. Obtain the Nusselt number for constant wall temperature and constant wall heat flux cases for fully developed flow in terms of temperature and velocity profiles in a circular pipe. Assume laminar, incompressible, twodimensional and steady state. 2.17. Develop analytical expressions for the Nusselt number for constant wall temperature and constant heat flux cases in a circular pipe with fully developed velocity profile but developing temperature profile. Assume laminar, incompressible, two-dimensional and steady state with uniform inlet velocity. 2.18. Consider steady state, two-dimensional, fully developed thermal and hydrodynamic laminar flow with constant properties in a circular pipe. Let there be heat transfer to or from the fluid at a constant rate per unit of tube length. Determine the Nusselt number if the effect of frictional heating due to viscosity (viscous dissipation) is included in the analysis. How does frictional heating affect the Nusselt number? Identify the pertinent nondimensional parameters appearing in this problem. 2.19. The Couette flow was discussed in Example 2.4. When both surfaces are stationary and flow is caused by an applied pressure gradient in the xdirection, dp/dx, it is called Poiseuille flow (see Fig. P2.8). The fluid is Newtonian and incompressible. Temperature in the upper and lower plates are kept constant at T1 and T2 respectively. Determine the fluid temperature distribution and heat flux with the following assumptions (Poiseuille flow): (1) steady state, (2) constant properties, (3) laminar flow, (4) constant ∂p/∂x, (5) no end effects, (6) no edge effects, and (7) gravity is the only body force. T1 u(y) z,w y,v x, u 2L g T2 Figure P2.8 2.20. Consider the steady flow of an incompressible fluid making an angle θ with the vertical direction with constant liquid film thickness δ along an inclined flat surface with length L and width W as shown in Figure 2.9. Determine the average velocity, volumetric flow rate, film thickness δ, and the force of the fluid on the surface. Neglect the end, exit, and edge effects and assume constant properties. The air is stagnant next to the liquid free surface. 196 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press y u(y) δ x θ g τyx(y) L Figure P2.9 2.21. Develop the energy equation and boundary conditions, for Problem 2.21 for the case of constant wall temperature Tw and uniform inlet temperature T0. Determine the temperature distribution and heat transfer rate for the wall to film for startup period (short contact time) that the fluid temperature changes only in the vicinity of the wall. Liquid free surface u0 h0 Figure P2 10 2.22. Develop formulation for the case of heating with no evaporation for the flow configuration shown below over a rotating disk (Fig. P2.10) Assume that the flow is laminar, incompressible, steady, axisymmetric and the disk has a constant wall temperature. The flow is introduced from a central collar that directs the liquid radially outward with velocity u0 over a gap height of h0. δ is the liquid film thickness and r0 is the initial radius of the disk at the outlet of the collar, h0 is the collar height and δT is the thickness of the thermal boundary layer. Obtain Nusselt number Nu = hr0 / k , and film thickness δ / h0 as a function of the Reynolds number and 2 Roosby numbers [Re = u0 r0 / ν and Ro = u0 / (ω 2r02 )] using the integral method. Chapter 2 Ge era i ed G ver i g Equati s 197 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2.23. The generalized mass balance on a liquid-vapor interface is given in eq. (2.197). Rewrite it in a three-dimensional Cartesian coordinate system. 2.24. Equation (2.197) is valid for an interface between liquid and vapor phases. Derive the conservation of mass for solid-liquid interface that is applicable to the melting problem. 2.25. The jump conditions at an interface can also be obtained by applying the various conservation laws to a control volume that includes the interface. Apply the conservation of mass principle to the control volume shown in Fig. P2.11, and show that when the thickness of the control volume Δx goes to zero, the conservation of mass is ρ1 (u1 − uI ) = ρ2 (u2 − uI ) . Liquid R R0 Tc r Solid Figure P2.11 Figure P2.12 2.26. For the control volume shown in Fig. P2.11, show that the momentum balance at the interface is p1 − p2 = ρ2 (u2 − uI )u2 − ρ1 (u1 − uI )u1 when the thickness of the control volume Δx goes to zero. 2.27. The momentum balance equation at a three-dimensional interface was discussed in Section 2.3.7. Derive the momentum balance equation for a two-dimensional interface between phases 1 and 2. The effect of surface tension can be neglected. 2.28. A tube with radius Ro and zero thickness passes through a liquid phase change material at its melting point, Tm (see Fig. P2.12). A cold fluid with a temperature of Tc, which is below Tm, flows through the inside of the tube. Solidification will occur on the outer surface of the tube. Formulate the problem by giving the governing equations and corresponding boundary conditions in the cylindrical coordinate system. If the problem 198 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press can be assumed to be quasisteady state, solve for the instantaneous location of the solid-liquid interface. 2.29. In a counter-current condenser shown in Fig. P2.13, the liquid flows  downward with a mass flow rate of m , while the vapor flows upward  ′′ with a mass flow rate of mv . The heat flux at the external wall is qw . Derive the jump conditions using mass and energy balances at the liquidvapor interface. Vapor z  m ( z ) ′′ qw  mv ( z ) z  ′′ mδ  m ( z + dz )  m v ( z + dz ) ′′ qw z+dz Liquid Figure P2.13 Liquid 2.30. Redo the Problem 2.29 for the case of co-current condensation, i.e., both liquid and vapor flow downward. 2.31. Redo Problem 2.29 for the case of counter-current evaporation and discuss the effect of the heat transfer direction on the energy balance at the interface. 2.32. The time it takes to completely burn a fuel droplet of specified size as estimated by eq. (2.310) is valid for conduction-controlled combustion processes. For a liquid droplet that has sufficient inertia, the heat transfer between vapor and liquid is dominated by convection, i.e., eq. (2.304) is f replaced by m′′ hv = qI ′′ = hθ∞ , where h is the heat transfer coefficient between the gaseous mixture and liquid droplet, and can be obtained by NuD = hD / k = C Re1/2 . What is the time it takes to completely burn the D liquid fuel droplet? Compare your result with eq. (2.310). 2.33. Redo Example 2.8 using the result from Problem 2.32 and compare your result with the result in Example 2.8. The relative velocity of the fuel Chapter 2 Ge era i ed G ver i g Equati s 199 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press droplet is V p = 5m/s and the heat transfer correlation is NuD = 0.6 Re1/2 . D The viscosity of the gaseous mixture is 4 × 10−5 m 2 /s . 2.34. A fuel droplet with an initial diameter of Dp moves at a velocity vp in an oxidant gas. The gravitational acceleration vector is g and the local oxidant gas velocity vector is vg. The fuel droplet evaporates as it moves in the combustor. If the drag coefficient of the fuel droplet is CD, write the momentum equation for the fuel droplet. 2.35. If the temperature of the fuel droplet in Problem 2.34 can be assumed to be uniform at any time, and the latent heat of evaporation of the fuel is hs , derive the energy equation for the fuel droplet. 2.36. A saturated mixture of water and its vapor at 180 °C flows in a 0.1 m ID tube with a mass flow rate of 0.5 kg/s. The liquid water is dispersed in the vapor phase in the form of 0.1-mm-diameter droplets, and the quality of the mixture is x=0.75. The average velocity of the vapor phase is 30m/s. Find: (a) the average velocity of the liquid droplets, and (b) the interactive force between the liquid and vapor phase. 2.37. If the liquid droplet is subcooled at a temperature of 175 °C, estimate the interphase heat transfer rate between liquid and vapor phases in the above problem. 2.38. The mass flux in the channel is defined as the mass flow rate per unit area. What is the total mass flow rate for the multiphase flow discussed in Example 2.11? Figure P2.14 200 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2.39. A capillary porous structure with partial heating and evaporation on its upper surface is shown in Fig. P2.14. The entire porous structure is saturated with liquid from the bottom surface (y=0), which is connected to ′′ a pool for liquid supply. A constant heat flux, q0 , is applied over part of the upper surface (0<x<Lxf), which is impermeable to the fluid. The rest of the upper surface (Lxf<x<Lx) exposes liquid in the pores of the capillary structure to the vapor space above. Heat is transferred from the left portion of the upper surface to the right portion of the upper surface, where evaporation takes place. The left- and right-side walls are adiabatic and impermeable. At steady-state, the liquid is drawn from the bottom of the porous structure and flows to the right portion of the upper surface due to evaporation. If the liquid velocity at the bottom surface and the vapor velocity at the upper surface are uniform, determine the liquid velocity at the bottom, v , and the vapor velocity at the upper surface, vv . 2.40. The area-averaged velocity components in the x- and y-directions in Problem 2.39 can be estimated by Darcy’s law: u=− K ∂p μ ∂x v=− K ∂p μ ∂y where K is permeability. The continuity equation is ∂u ∂v + =0 ∂x ∂y Determine the analytical solution for the pressure and velocity distributions in the porous structure. 2.41. Specify the energy equation and the corresponding boundary conditions for Problem 2.38; then nondimensionalize them. Write a computer program to obtain dimensionless temperature distribution in the porous media. Quadrant used to model elemental volume y,v z,w x,u Figure P2.15 Chapter 2 Ge era i ed G ver i g Equati s 201 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2.42. For a porous bed made of spherical balls packed in a body-centered-cubic (BCC) formation, use the continuity equation (2.54) and momentum equations (2.74) – (2.76), to develop the boundary conditions and numerically solve the flow field, and determine the permeability of the packed bed by plotting the product of μv vs. dp/dx. Assume that the elemental cube’s sides all have a length of 10-3 m and that the spheres all have a radius of 4.6×10-4 m. (Note: the spheres are slightly compressed against each other, and therefore their radii are slightly larger than if they were only contacting one another at a single point). Model the top right quadrant of an elementary cubic cell, Fig. P2.15. Numerically solve the problem for different fluids and determine the permeability. 2.43. For most heat pipes, the wick structures are very thin, so the liquid flow in the wick can be simplified to one-dimensional axial flow. Since the liquid velocity and its gradient in the axial direction are very small, the liquid  , where K is μ dz permeability. For an axially-grooved wick structure as shown in Fig. P2.16, determine the permeability analytically. flow in the wick satisfies Darcy’s law, i.e., w = − K dp Figure P2.16 Figure P2.17 2.44. The designer of a coffee maker wants to make a filtration system in which the water is heated and filtered at the same time. The filter can be modeled as a porous zone within a tube, and the heat can be modeled as a uniform heat flux applied to the tube’s outer wall. Set up the governing equations and boundary conditions so that the designer can get a temperature and pressure field, and thus know where the liquid may start to boil. 2.45. A conventional heat pipe with multiple heat sources, as shown in Fig. P2.17, is divided into three regions: vapor-flow region, liquid-wick region, and solid wall. The effect of gravity can be neglected so that the fluid flow 202 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press and heat transfer can be treated as two-dimensional. Specify governing equations and corresponding boundary conditions for all three regions. 2.46. The boundary of a porous zone is subjected to a constant heat flux. The porous zone is in local thermal nonequilibrium. Specify the thermal boundary conditions for the fluid phase and the thermal boundary condition for the solid phase. Also discuss the physical significance of your boundary conditions. Inlet 2: Liquid water/air mixture entering velocity of each component known g y, v x, u Outlet: Constant pressure Inlet 1: Liquid water/air mixture entering velocity of each component known Figure P2.18 2.47. Consider flow entering from different sides of a tee-junction in a duct. The fluid entering is a mixture of air and liquid water, with known velocities and volume fraction at each inlet. The outlet is held at a constant pressure. A schematic of this problem is presented in Fig. P2.18. Model the problem using a two-dimensional Cartesian homogeneous model to find the mass averaged velocity. Assume that the flow is in the laminar regime and that gravitational effects are important. Also, assume that the flow field is isothermal. The volume fraction of air is small; therefore, the air exists in the form of bubbles. The diameter of these bubbles is known from experimental observation. 2.48. A parametric study on the effect of the inlet velocities and gas volume fraction from Problem 2.46 is needed. However, before performing this study using the homogeneous model, the inaccuracies of the model must Chapter 2 Ge era i ed G ver i g Equati s 203 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press first be quantified. The inaccuracies should be quantified by comparing the homogeneous results to a volume-averaged multifluid model. Therefore, set up the governing equation using a volume-averaged multifluid model with the corresponding boundary conditions with same assumptions used in Problem 2.47. 2.49. There are many processes where large pressure drops in a liquid occur. If the pressure is low enough at one side, the liquid may locally vaporize. This process is called cavitation. Cavitation may occur in a liquid flowing past a sharp-edged circular orifice. The schematic of this problem is presented in Fig. P2.19. Liquid enters the orifice from the left. At the orifice there is a sharp and large pressure drop, causing vapor bubbles to form. The walls of the orifice are well insulated and can be considered adiabatic. The flow is inherently turbulent; therefore, use the k − ε model to model turbulent kinetic energy and dissipation. Develop the twodimensional steady state governing equations using a homogeneous volume-averaged model. Assume that the vapor and liquid have the same velocities; therefore the mass-averaged velocity of the mixture is identically equal to each intrinsic phase average velocity. Liquid In r, v x, u pin Adiabatic Wall Liquid and Vapor Out Axis pout> pin Figure P2.19 2.50. Consider the Direct Methanol Fuel Cell (DMFC) with a liquid feed presented in Fig. P2.20. Assume that the methanol solution is passing over the fuel cell at a very high velocity, and that solution has a constant concentration at the top of the anode gas diffusion layer. The oxygen supplied to the cathode comes from air. Assume that the catalyst layers are very thin, therefore their thickness can be neglected and the reaction rates can be considered to be a surface reaction rate. The membrane can be considered saturated at all times. Develop the governing equations and boundary conditions for the transport of liquid and gas through the fuel cell as well as the electric potentials of the membrane and gas diffusion layers for a one-dimensional steady state simulation. 204 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press Methanol Solution x Catalyst Layers Anode Gas Diffusion Layer Membrane Cathode Gas Diffusion Layer Air Figure P2.20 Current Collectors 2.51. The multiphase transport model in porous media can also be used to model two-phase flow in microchannel heat exchangers. The transport phenomena in a microchannel heat exchanger can be characterized as unidirectional flow combined with three-dimensional conduction. Set up the governing equations for this system, while leaving the terms caused by phase interaction and solid/fluid interaction unclosed. The energy equation shall be written in terms of internal energy, and the solid, liquid and vapor equations shall be set up separately. z L Gas Liquid layer r Figure P2.21 2.52. Consider a case of annular flow in which the liquid flows as a film on the wall of a pipe while the gas flows down the central cylindrical core. For a Chapter 2 Ge era i ed G ver i g Equati s 205 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press pipe of diameter D let the interfacial and wall shear stresses be τ δ and τ w as given inputs, respectively. Assuming vertical co-current flow configuration as shown in Fig. P2.21, develop the formulation for onedimensional, unsteady, laminar flow with multifluid modeling. 2.53. In a fluidized bed particles are supported by an upward flow of fluid around them as shown in Fig. P2.22. Formulate the problem as multifluid model and deduce the fluid flux and pressure gradient needed to cause fluidization in a bed with void fraction εf by neglecting inertia. Assume one-dimensional, incompressible, steady flow with negligible interparticle forces and no phase change. z Figure P2.22 2.54. Deduce the general one-dimensional governing equations for mass and momentum for stratified turbulent flow of two phases (1) and (2) with evaporation over a flat plate as shown in Fig. P2.23. 2.55. Formulate the analysis for heat transfer for Problem 2.54 using drift flux model. Assume steady, incompressible flow with negligible viscous dissipation. Use the results obtained in Problem 2.54 in carrying out the analysis. z Steam flow (1) Liquid film (2) Figure P2.23 2.56. Formulate steady flow in a one-dimensional nozzle under compressible inviscid conditions. The problem should be modeled as a homogeneous two-component flow with no phase change. Neglect the effects of surface tension. Obtain velocity and volume fraction distribution along the nozzle. 2.57. Derive the equations for mass, momentum, and energy for the flow of compressible gas through a channel containing dispersed particles. The flow occurs at high speed. Neglect the following effects: exchange of 206 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press momentum, heat and mass with the channel walls, phase change, chemical reactions, viscous drag and viscous dissipation and gravity. Assume the flow is one-dimensional, steady, and constant particle phase density. Do the analysis using the area-averaged multifluid models for channel flows. 2.58. Electro-osmosis is flow induced by an external electric field. When surfaces acquire a finite charge density and are in contact with a polar solution, they induce a distribution of electrical charges within the electrically charged solution. Electro-osmosis is used in biochip technology as a pumping mechanism in MEMS devices for chemical and biological analysis and diagnoses. Two-fluid electro-osmotic flow is used in many microfluids. The working fluid needs to be a conducting liquid with significant electrical conductivity to function properly. As shown in Fig. P2.24, a high electro-osmotic mobility liquid is at the bottom of the channel, and a low electro-osmotic mobility liquid is in the upper section of the microchannel. When an electric field is applied across the channel, pressure and electro-osmosis effects drive the liquids. The flow depends on the viscosity ratio of the two fluids, the external electric field, electro-osmosis characteristics of the high mobility fluid, and the interfacial curvature between the two fluids. Develop the steady formulation for velocity and the electric potential due to an electrically applied field caused by charges on the wall. y Liquid 2 u2,e Interface Liquid 1 u1,e δe δ0 H x Ex Figure P2.24 2.59. Derive eq. (2.520) by multiplying temperature to the continuity equation (2.503) and substituting eq. (2.491), (2.494), (2.495), and (2.497) into the resultant equation. Chapter 2 Ge era i ed G ver i g Equati s 207 Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press 2.60. Show that the time-averaged conservation of species mass in a binary system is eq. (2.522). 2.61. Multiply u′2 to the continuity equation (2.503) and time average the resultant equation to show that  ∂u′ ∂v′ ∂w′  u ′2  + + =0  ∂x ∂y ∂z  2.62. In the gravitational field, the effect of gravity must be included. If the direction of the gravity is along the negative z-direction, the momentum equation in the z-direction becomes  ∂2w ∂2w ∂2w  ∂w ∂w ∂w ∂w 1 ∂p +u +v +w =− + ν  2 + 2 + 2  − g[1 − β (T − T0 )]  ∂x ∂t ∂x ∂y ∂z ρ ∂z ∂y ∂z    which will replace eq. (2.506). Rederive the turbulent kinetic energy equation (2.528) for this case. 208 Adva ced Heat a d ass Tra sfer Amir Faghri, Yuwen Zhang, and John Howell Copyright © 2010 Global Digital Press