Back to Transport Phenomena in Multiphase Systems Home
Table of Contents
5 SOLID-LIQUID-VAPOR PHENOMENA AND INTERFACIAL HEAT AND MASS TRANSFER 5.1 Introduction The interfacial region between two homogeneous phases contains matter in a distinct physical state; that is to say, matter in the interfacial state exhibits properties different from those matters in the gaseous, liquid, or solid states. As a result, as soon as interfaces are considered explicitly, new variables – for example, interfacial surface tension – enter into the classical thermodynamic description of equilibrium systems. Interfaces in equilibrium systems need not be considered explicitly unless the surface-to-volume ratio is large, because the contribution of interfacial free energy to the total free energy is usually small. However, interfacial effects on the dynamic behavior of flow systems can be profound, even when the proportion of matter in interfacial regions is extremely small. Furthermore, motion may originate in an interface in systems that are not in thermal or compositional equilibrium. When two adjacent fluids are at rest, their interface ordinarily behaves as if it is in a state of uniform tension. It is both possible and convenient to mathematically represent the interface as a geometric surface in tension. This representation is also appropriate for many flows with free boundaries; indeed, it is the basis of the treatment of capillarity in classical hydrodynamics. Considerations of equilibrium surface tension lead to the conclusion that the normal component of fluid stress, or pressure, is discontinuous at a curved interface, while the shear stress is continuous. Classical hydrodynamics also recognizes – in connection with the calming action of oil on water waves – that extension and contraction of a surface film produces longitudinal variations in surface tension. This in turn gives rise to discontinuities in the tangential components of fluid stress at the interface. A design engineer must have some understanding of these and other phenomena to design effective devices. The field of interfacial phenomena has been the realm of researchers primarily in chemical and mechanical engineering, Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 331 physical chemistry, and material science. Much of the analytical basis for their work comes from a subcontinuum view of the physical world. Models of molecular interaction and the use of statistical mechanics are typical in the literature. The earliest practical work on interfacial phenomena used equilibrium thermodynamics. Today, the ad hoc use of thermodynamics and simple molecular interaction models constitute the most useful treatment of two-phase heat transfer problems. However, a significant effort has been made to develop a unified approach via the extension of continuum mechanics using the conservation equations and some form of the flux laws. Direct solution methods for interfacial phenomena do not always follow from this type of investigation, but one frequently gains valuable insight into the behavior at the interface. The assumption is that one is not interested in the precise details of the interfacial region but rather in how it affects the bulk regions. Section 5.2 discusses capillary pressure and interface shape at equilibrium. Also considered in Section 5.2 are the effects of interfacial tension gradient on the fluid flow and the temperature-dependence of surface tension. Wetting phenomena and contact angles in solid-liquid-vapor systems are taken up in Section 5.3. Section 5.4 address phase equilibrium in a microscale interfacial system; this includes the effects of disjoining pressure on ultra-thin liquid films and the change in saturated vapor pressure over a curved interface due to the combined effects of surface tension and disjoining pressure. Section 5.5 introduces interfacial mass, momentum, and energy balance as well as thermal resistance across ultra-thin films, including the effect on vaporization and condensation. The instability of liquid-vapor interface and waves on liquid film are discussed in Section 5.6. Finally, Section 5.7 discusses numerical simulation of interfaces and free surfaces using both continuum and noncontinuum approaches. 5.2 Surface Tension 5.2.1 Capillary Pressure: The Young-Laplace Equation Since the distance between the molecules in the vapor phase is much greater than that in the liquid phase, the intermolecular force between the molecules in the vapor phase is very weak. The intermolecular attractive force in the liquid phase holds the molecules in the liquid close to each other. For the molecules within the liquid phase, the intermolecular forces from all directions are balanced. Although the forces acting on the molecules at the liquid-vapor interface are balanced along the tangential direction, the attractive force from the molecules in the liquid phase, Fi (in normal direction), tends to pull the molecules at the liquid-vapor interface toward the liquid phase because the attractive force from the vapor phase, Fo, is much weaker. The net inward force 332 Transport Phenomena in Multiphase Systems Vapor Fo Fs Fi Fs Liquid Figure 5.1 Origin of surface tension at liquid-vapor interface. Fi − Fo causes movement of the liquid molecules until the maximum number of molecules is in the interior, which leads to an interface of minimum area (Fig. 5.1). It is generally necessary to specify two radii of curvature to describe an arbitrarily-curved surface, RI and RII, as shown in Fig. 5.2. The surface section is taken to be small enough that RI and RII are approximately constant. If the surface is now displaced outward by a small distance, the change in area is ΔA = ( x + dx )( y + dy ) − xy (5.1) If dxdy ≈ 0 , then ΔA = y dx + x dy (5.2) The work required to displace the surface is obtained from eq. (2.218), i.e., δ W = σ ( x dy + y dx) (5.3) Displacement acting on the area xy over the distance dz also creates a pressure difference Δp across the surface – capillary pressure (pcap). The work attributed to generating this pressure difference is δ W = Δp xy dz = pcap xy dz (5.4) From the geometry of Fig. 5.2, it follows that x + dx x = RI + dz RI or dx = x dz RI y dz RII (5.5) (5.6) Similarly, dy = (5.7) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 333 Figure 5.2 Arbitrarily-curved surface with two radii of curvature RI and RII. For the surface to be in equilibrium across this differential change, the two expressions for the work must be equal: σ ( x dy + y dx) = Δp xy dz i.e., (5.8) σ¨ § xy dz xy dz · + ¸ = Δp xy dz RII ¹ © RI (5.9) The pressure difference between two phases becomes §1 1· pcap = Δp = σ ¨ + ¸ = σ ( K1 + K 2 ) © RI RII ¹ (5.10) where K1 and K2 are curvatures of the surface. This expression is called the Young-Laplace equation, and it is the fundamental equation for capillary pressure. It can be seen that when the two curvature radii are equal, in which case the curved surface is spherical, eq. (5.10) can be reduced to eq. (2.239). 5.2.2 Interface Shapes at Equilibrium Surface tension effects on the shape of liquid-vapor interfaces can be demonstrated by considering a free liquid surface of a completely wetting liquid meeting a planar vertical wall (see Fig. 5.3). Since the interface is twodimensional, the second principal radius of curvature is infinite. For the first principal radius of curvature RI, it follows from analytical geometry that 334 Transport Phenomena in Multiphase Systems . g y Figure 5.3 Shape of the liquid-vapor interface near a vertical wall. 1 d 2 z / dy 2 (5.11) = RI ª1 + (dz / dy ) 2 º 3/ 2 ¬ ¼ The Young-Laplace equation, eq. (5.10), becomes pv − pA = σ / RI (5.12) At a point on the interface that is far away from the wall ( y → ∞, z = 0 ), the liquid and vapor pressures are the same, i.e., pv (∞,0) = pA (∞,0) (5.13) The liquid and vapor pressures near the wall ( z > 0 ) are pv ( z ) = pv (∞,0) − ρv gz (5.14) pA ( z ) = pA (∞,0) − ρA gz (5.15) The relationship between the pressures in two phases can be obtained by combining eqs. (5.13) – (5.15), i.e., pv − pA = g ( ρA − ρv ) z (5.16) Combining eqs. (5.12) and (5.16), the equation for the interface shape becomes ª § dz · 2 º d2z − «1 + ¨ ¸ » =0 σ dy 2 « © dy ¹ » ¬ ¼ Multiplying this equation by dz/dy and integrating gives g ( ρA − ρv ) z −3 / 2 (5.17) g ( ρA − ρv ) z 2 σ ª § dz · 2 º + «1 + ¨ ¸ » « © dy ¹ » ¬ ¼ −1 / 2 = C1 (5.18) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 335 Since at y → ∞ both z and dz/dy equal zero, C1 = 1. The boundary condition for equation (5.18) is (dz / dy ) y =0 → ∞ (5.19) Equations (5.18) and (5.19) can then be solved for z at y=0: ª º 2σ (5.20) z0 = z (0) = « » ¬ ( ρA − ρv ) g ¼ Using eq. (5.20) as a boundary condition, integrating eq. (5.18) yields the following relation for the shape of the interface: 1/ 2 § 2L yª = « cos h ¨ c Lc « © z0 ¬ where 2 § ·º z0 · ª § 2 Lc · º ¸ » − «cos h ¨ ¸» + ¨ 4 + 2 ¸ Lc ¹ © z ¹¼ ¬ ¹» © ¼ 1/ 2 −1 −1 1/ 2 § z2 · −¨4+ 2 ¸ Lc ¹ © 1/ 2 (5.21) ª º σ Lc = « (5.22) » ¬ ( ρA − ρv ) g ¼ is a characteristic length for the capillary scale. For small tubes with R  Lc , the interfacial radius of curvature is approximately constant along the interface and equal to R / cosθ , where is the contact angle. Example 5.1 Two parallel plates are put into bulk water at the bottom end (see Fig. 5.4). Estimate the minimal distance between the plates at which the central point of the liquid-air interface between the plates is not elevated from the bulk water level at equilibrium condition. Assume that water completely wets the material of the plates. The system temperature is 20 °C, σ = 0.07288N/m, and ρA = 999kg / m . 2y z y water Figure 5.4 Two parallel plates in bulk water. 336 Transport Phenomena in Multiphase Systems Solution: Since the density of vapor is much less than that of the liquid – ρv  ρA – eq. (5.22) can be simplified to obtain ª º ªσ º σ ª 0.07288 º −3 Lc = « » « » =« » = 2.728 × 10 m ρA g ¼ ( ρA − ρv ) g ¼ ¬ 999 × 9.8 ¼ ¬ ¬ The rise of the liquid surface at y=0 is obtained from eq. (5.20): 1/ 2 1/ 2 1/ 2 ª º 2σ z0 = z (0) = « » ¬ ( ρA − ρv ) g ¼ 1/ 2 1/ 2 ª 2σ º =« » ¬ ρA g ¼ 1/ 2 ª 2 × 0.07288 º −3 =« » = 3.859 × 10 m ¬ 999 × 9.8 ¼ The value of y for z = 0 can be found by using eq. (5.21), i.e., § § 2L ·º ª z2 · yª § 2L ·º = « cos h ¨ c ¸ » − «cos h ¨ c ¸ » + ¨ 4 + 0 ¸ Lc ¬ L2 ¹ © z ¹¼ « » ¬ © z0 ¹ ¼ c © −1 −1 −1 1/ 2 § z2 · −¨4+ 2 ¸ Lc ¹ © −1 1/ 2 ª ª § 2 × 2.728 × 10−3 · º § 2 × 2.728 × 10−3 · º = « cos h ¨ − «cos h ¨ ¸» −3 ¸ » 0 « « © 3.859 × 10 ¹» © ¹» ¬ ¼ ¬ ¼ 2 § § 3.859 × 10−3 · · ¸ +¨4+¨ −3 ¸ ¨ ¸ © 2.728 × 10 ¹ ¹ © 1/ 2 2 § 0 § ·· −¨4+¨ ¸ −3 ¸ ¨ © 2.728 × 10 ¹ ¸ © ¹ 1/ 2 = 0.909 i.e., y = 0.909 Lc = 0.909 × 2.728 × 10−3 = 2.48 × 10−3 m = 2.48mm This value of y is one-half the minimal distance between the plates. Therefore, the minimal distance is about 5.96 mm. 5.2.3 Effects of Interfacial Tension Gradients Since surface tension depends on temperature, a permanent nonuniformity of temperature or concentration (for a multicomponent system) at a liquid-vapor interface causes a surface tension gradient. The interfacial area with small surface tension expands at the expense of an area with greater surface tension, which in turn establishes a steady flow pattern in the liquid; this flow caused by the surface tension gradient along the liquid-vapor (gas) interface is referred to as the Marangoni effect. The surface tension of a multicomponent liquid that is in equilibrium with the vapor is a function of temperature and composition of the mixture, i.e., σ = σ (T , x1 , x2 ,", xN −1 ) (5.23) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 337 where xi is the molar fraction of the ith component in the liquid phase and N is the total number of components in the liquid phase. The change of surface tension can be caused by either change of temperature or composition, i.e., N −1 § ∂σ · § ∂σ · dσ = ¨ dxi (5.24) ¸ ¸ dT + ¨ © ∂T ¹ xi i =1 © ∂xi ¹T , x ¦ j ≠i Recall from thermodynamics that as the critical temperature for a given fluid is approached, the properties of the liquid and vapor phases of the fluid become identical, i.e., σ vanishes. It therefore follows that the surface tension decreases with temperature, i.e., [ (∂σ / ∂T ) xi < 0 ]. For a pure substance, surface tension is a function of temperature only. The curve-fit equations for surface tension are almost linear for most fluids and can have the form σ = C0 − C1T (5.25) where C0 and C1 are empirical constants that vary between substances. For water, C0 = 75.83 × 10−3 N/m and C1 = 0.1477 × 10−3 N/m- o C . The unit of the temperature T in eq. (5.25) is °C. Since surface tension varies with temperature, the surface tension will not be uniform if the temperature, along the liquid-vapor interface, is nonuniform. Consequently, the liquid in the lower surface tension region near the interface will be pulled toward the region with higher surface tension. Furthermore, because the surface tension usually decreases with increasing temperature, the flow driven by surface tension moves away from the interface with high temperature and towards the interface with low temperature. As noted before, this motion of a liquid caused by a surface tension gradient at the interface is referred to as the Marangoni effect. The most well-known example of surface-tensiondriven flow is Bernard cellular flow, which occurs in a thin horizontal liquid layer heated from below. Figure 5.5 illustrates steady cellular flow driven by the Marangoni effect. Once a steady Marangoni flow is established, the liquid Gas @ TG A hδ h B y, v δ x, u TW TW TA TB TG Figure 5.5 Marangoni effect: cellular flow driven by surface tension gradient. 338 Transport Phenomena in Multiphase Systems velocity is upward at point A and downward at point B. Since the liquid at surface point A comes directly from the hot surface, the temperature at point A is higher than that at point B. Consequently, the surface tension at point A is smaller than that at point B, and the fluid at point A is pulled toward point B. Although a temperature gradient exists in the vertical direction, the actual driving force is the surface tension gradient in the horizontal direction. The liquid surface loses heat to the gas phase at temperature, TG, through convection with a heat transfer coefficient of hδ. An important boundary condition at the interface explains the liquid flow field caused by surface tension variation along the interface: § ∂u · § ∂σ ·§ ∂T · τ yx = μA ¨ ¸ =¨ (5.26) ¸¨ ¸ y =δ © ∂y ¹ y =δ © ∂T ¹© ∂x ¹ y =δ Bernard cellular flow resulting from the Marangoni effect can occur only if certain conditions are met. The conditions for the onset of cellular motion can be predicted using linear stability analysis as presented by Carey (1992). The local temperature is assumed to equal the sum of the basic temperature T – given by the initial linear profile – and a small sinusoidal fluctuation T ′ that represents a Fourier component of random disturbances, i.e., T = T + T ′ = (Tw − ζ y ) + T ′ (5.27) where ζ = (Tw − Tδ ) / δ . Since all base flow velocities are zero, the velocities are v = v + v′ = v′ (5.28) ′ = u′ u =u +u (5.29) Substituting eqs. (5.27) – (5.29) into the continuity, momentum, and energy equations, and subtracting the corresponding base flow equations, the resulting equations are solved assuming the following forms for T ′ and v′ : T ′ = θ ( y )eiα x + β t (5.30) v′ = V ( y )eiα x + β t (5.31) where α = 2π / λ is the wave number and β can be a complex number with its real part representing the amplification factor, and its imaginary part representing the temporal frequency. The resulting stationary wave solutions (for β = 0 ), which are not amplified or dampened with time, are therefore interpreted as stable perturbations. The conditions for which these solutions are obtained are assumed to correspond to the onset of instability, which leads to cellular convection. The case of marginal stability corresponds to ˆˆ ˆ ˆˆ ˆ ˆ 8α (α cosh α + Bi sinh α )(α − sinh α cosh α ) (5.32) Ma = 3 3 5 ˆ ˆ ˆ ˆ ˆ ˆ α cosh α − sinh α − (8Crα cosh α ) /(Bo + α 2 ) where Ma is the Marangoni number: ζ (dσ / dT )δ 2 Ma = (5.33) α A μA Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 339 and α A is the thermal diffusivity of the liquid. ˆ The other dimensionless numbers in eq. (5.32) are the wave number α , the Biot number Bi, the Bond number Bo, and the Crispation number Cr. Their definitions are ˆ α = αδ (5.34) hδ δ Bi = (5.35) kA Bo = ( ρ A − ρ g ) gδ 2 σ μAα A Cr = σδ (5.36) (5.37) The liquid film is stable if the Marangoni number is below that predicted by eq. (5.32). However, the liquid film is not stable if its Marangoni number is above that obtained by eq. (5.32). The marginal stability predicted by eq. (5.32) for some typical combinations of parameters is illustrated in Fig. 5.6. Since the Fourier components of all wavelengths can be contained in a random disturbance, the system becomes unstable when it is unstable at any wavelength. For a system with Cr < 10−4 and Bi → 0, Bo → 0 , Fig 5.6 indicates that the critical Marangoni number is about 80 and the associated dimensionless wave number ˆ α = 2. Figure 5.6 Stability plane for the onset of cellular motion (Carey, 1992; Reproduced by permission of Routledge/Taylor & Francis Group, LLC). 340 Transport Phenomena in Multiphase Systems The Marangoni effect can have an important influence on heat and mass transfer processes, including the evaporation of a falling film, as well as causing vapor bubbles in a liquid with a temperature gradient to move toward the high temperature region during boiling. Example 5.2 A 0.3-mm-thick water film sits on a surface held at a temperature of Tw = 80 °C. The top of the liquid film is exposed to air at a bulk temperature of TG = 20 °C, and the convective heat transfer coefficient between the liquid film and the air is hδ = 10 W/m2K. Determine whether the water film is stable. Solution: Since the liquid film is very thin, the temperature drop across the liquid film will be very small. The properties of the liquid film can be determined at the wall temperature of 80 o C , i.e., kA = 0.67W/m-K, ρA = 971.9kg/m3 , α A = 1.64 × 10−7 m 2 /s, μA = 351.1 × 10−6 N-s/m 2 , σ = 0.0626N/m, and dσ / dT = −1.7 × 10−4 N/m o C. The Biot number is h δ 10 × 0.3 × 10−3 Bi = δ = = 0.00425 kA 0.67 At steady-state, the surface temperature of the interface, T , satisfies T −T kA w δ = hδ (Tδ − TG ) δ i.e., § hδ · § hδ Tδ = ¨ Tw + δ TG ¸ ¨1 + δ kA kA © ¹© = · Tw + BiTG ¸= 1 + Bi ¹ 80 + 4.25 × 10−3 × 20 = 79.75o C 1 + 4.25 × 10−3 which demonstrates that the temperature drop across the liquid film is minimal. The Bond number is obtained by using eq. (5.37), i.e., ( ρA − ρ g ) gδ 2 ρA gδ 2 971.9 × 9.8 × 0.00032 Bo = = = 0.014  σ σ 0.0626 The Crispation number Cr is μ α 351.1 × 10−6 × 1.64 × 10−7 Cr = A A = = 3.07 × 10−6 σδ 0.0626 × 0.0003 which is less than 10-4. The critical Marangoni number, Mac, below which the liquid film is stable, can be obtained from Fig. 5.4; its value is ˆ 80 at α = 2 . The Marangoni number of the system can be obtained from eq. (5.33), i.e., Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 341 Ma = ζ (dσ / dT )δ 2 Tδ − Tw (dσ / dT )δ 2 = α A μA δ α A μA 79.75 − 80 (−1.7 × 10−4 ) × 0.00032 × = 221.4 > Ma c 0.0003 1.64 × 10−7 × 351.1 × 10−6 Therefore, the system is unstable and Marangoni convection will occur. = 5.3 Wetting Phenomena and Contact Angles 5.3.1 Equilibrium and Apparent Contact Angles In addition to the surface tension at a liquid-vapor ( Av ) interface discussed above, surface tensions can also exist at interfaces between solid-liquid interface ( sA ) and solid-vapor interface ( sv ); this can be demonstrated using a liquidvapor-solid system in Fig. 5.7. The contact line is the locus of points where the three phases intersect. The contact angle, θ , is the angle through the liquid between the tangent to the liquid-vapor interface and the tangent to the solid surface. The contact angle is defined for the equilibrium condition. In 1805, Young published the basic equation for the contact angle on a smooth, insoluble, and homogeneous solid: σ − σ sA cos θ = sv (5.38) σ Av which follows from a balance of the horizontal force components (Faghri, 1995), as shown in Fig. 5.7. When there is relative motion of a liquid drop over a solid surface, a different contact angle can be expected. When the relative motion stops, an angle different from the apparent (equilibrium) contact angle is seen; it depends upon the direction of the previous motion, i.e., whether it was a receding or advancing surface, as shown in Fig. 5.8. The minimum wetting contact angles for different solid-liquid combinations obtained by Stepanov et al. (1977) are reproduced in Table 5.1. All contact angle approaches are based on the following Figure 5.7 Drop of liquid on a planar surface. 342 Transport Phenomena in Multiphase Systems (a) (b) (c) Figure 5.8 Schematic of apparent contact angle : (a) stationary liquid, (b) liquid flows upward, (c) liquid flows downward. assumptions (Kwok and Neumann, 1999; Yang et al., 2003): (1) validity of Young’s equation (5.38), (2) pure liquid, (3) constant values of σ sv , σ sA , and σ Av , (4) the value of liquid surface tension which should be higher than the anticipated solid surface tension, and (5) a value of σ sv independent of the liquid used. Table 5.1 Minimum wetting contact angle in arc degree (the upper and lower values are for advancing and receding liquid front, respectively; Stepanov et al. 1977) Acetone Aluminum Beryllium Brass Copper Nickel Silver Steel Titanium 25/11 Water 73/34 63/7 82/35 84/33 79/34 63/38 72/40 73/40 Ethanol 0/0 18/8 15/7 16/7 14/7 19/8 18/8 R-113 16/7 14/6 16/5 For a rough surface, the contact angle θ rough is related to the contact angle on a smooth surface θ by cosθ rough = γ cos θ (5.39) where γ is the ratio of the rough surface area to the smooth surface area. Since γ is always greater than 1, cosθ rough is greater than cosθ . Therefore, the contact angle θ rough is less than the contact angle on a smooth surface, θ . Equilibrium contact angles can vary depending on the motion history of the contact line, particularly for rough surfaces. Equilibrium contact angles may be used in calculations for various heat transfer devices. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 343 5.3.2 Wettability and Adsorption Depending on the contact angle, liquids can be classified as nonwetting ( 90D < θ < 180D ), partially wetting ( 0D < θ < 90D ), or completely wetting ( θ = 0D ). When a small amount of liquid is brought into contact with an initially dry solid surface, the liquid behaves in one of two ways: (1) if the liquid does not wet the solid it may break up into small droplets, or (2) if the liquid wets the solid it may spread over the solid surface and form a thin liquid film. Wettability can be attributed to a strong intermolecular attractive force near the interface between the solid and liquid. The thermodynamic definition of surface tension, eq. (2.212), establishes that there is a significant decrease in the surface free energy per unit area in a wetting liquid. Spreading of a liquid on a solid surface can be described by the spreading coefficient Sp, defined as SpAs = σ sv − σ Av − σ sA (5.40) which is a measure of the ability of a liquid to spread over a solid surface. The spreading coefficient defined by eq. (5.40) is difficult to evaluate, because data for σ sv and σ sA are not available for most substances. Substituting eq. (5.38) into eq. (5.40), one obtains SpAs = −σ Av (1 − cosθ ) (5.41) For a partially wetting liquid ( 0D ≤ θ ≤ 90D ), cos θ ≤ 1 and the spreading coefficient SpAs is always negative, which means that the liquid will partially wet the solid surface and an equilibrium contact angle can be established. The discussion thus far has treated the three phases – solid, liquid, and vapor – as though their boundaries were sharply-delineated lines or surfaces. This idealization, which serves as a useful analytical device at the macroscopic level, does not hold at the microscopic level. At that level, the interfaces between phases appear as regions over which properties vary continuously, rather than as lines or surfaces with discontinuous property changes. Intermolecular forces of both repulsion and attraction influence how material in the various phases is distributed throughout these interfacial regions. Adsorption – which is one of the consequences of this intermolecular action – occurs when a liquid or solid phase adjacent to a second phase (solid, liquid, or gas) retains molecules, atoms, or ions of the second phase at the shared interface. Adsorption affects the wetting process because it alters the interfacial tension of the solid-liquid interface. Introducing the surface pressure of the adsorbed material on the solid surface, ps: ps = σ sv − σ sv , a (5.42) where σ sv ,a is the interfacial tension with the absorbed substance present. Young’s equation, eq. (5.38), can be rewritten as σ Av cosθ = (σ sv − ps ) − σ sA (5.43) which indicates that adsorption changes equilibrium contact angles and surface tension. 344 Transport Phenomena in Multiphase Systems (a) wetting liquid (b) nonwetting liquid Figure 5.9 Capillary phenomenon in an open tube (Faghri, 1995; Reproduced by permission of Routledge/Taylor & Francis Group, LLC). Surface tension, equilibrium contact angles, and capillary pressure determine the behavior of liquids in small-diameter tubes, slots, and porous media. The best example is the capillary rise of a wetting liquid in a small tube (see Fig. 5.9), in which case the capillary force is balanced by gravitational force. The pressure difference across the liquid-vapor interface in such tube is given by eq. (5.10). For tubes with a very small radius, the two radii of curvature RI and RII are the same, i.e., r (5.44) RI = RII = cos θ Substituting eq. (5.44) into eq. (5.10), one obtains the pressure difference between liquid and vapor at point C 2σ (5.45) ( pv − pA )C = pcap ,C = cos θ r The pressure at the flat surface A [see Fig. 5.9 (a)] is related to the vapor pressure by p A = pv + ρv gH (5.46) where H is the height of the meniscus above the flat surface A. The pressure at point B inside the tube must be equal to that at point A because they are on the same horizontal surface, i.e., pB = pv + ρA gH − pcap ,C = p A (5.47) Combining eqs. (5.45) – (5.47), one obtains 2σ cos θ = ( ρA − ρv ) gH (5.48) r Capillary rise phenomena can be observed when the liquid wets the tube wall ( θ < 90D ). If the liquid cannot wet the tube wall ( θ > 90D ), the capillary rise H Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 345 obtained by eq. (5.48) is negative, which indicates that there is a capillary depression, as shown in Fig. 5.9b. In general, solid materials have only two types of behavior when they interact with water. They are either hydrophobic or hydrophilic. Hydrophobic materials have little or no tendency to absorb water, while hydrophilic materials have an affinity for water and readily absorb it. The criteria for hydrophobic or hydrophilic properties of a material is based on its contact angle 2 2 Hydrophobic materials can be observed as beading of water on a surface, such as a freshly waxed car surface. Hydrophilic materials allow water to wet its surface forming a film or coating. Hydrophilic materials are usually charged or have polar side groups to their structure that attract water. There are many cases of hydrophobic surfaces in nature, including some plant leaves, butterfly wings, duck feathers and some insects’ exoskeletons. There are many synthetic hydrophobic materials available including waxes, alkanes, oils, Teflon, and Gortex. There are numerous applications for using these materials such as protection of stone, wood, and cement from the effects of rain, waterproofing fabrics and the removal of water from glass surfaces, such as a windshield, to increase transparency. Hydrophobic materials are also used for cleaning up oil spills, removal of oil from water and for chemical separation processes to remove nonpolar from polar compounds. A hydrophilic material’s ability to absorb and transport water gives it numerous applications in cleaners, housings, cables, tubes and hoses, waterproofing, catheters, surgical garments, etc. A hydrophilic coating on a tube or hose eliminates the need for other lubricants, which is useful to prevent crosscontamination. Hydrophilic coatings on plugs and o-rings increase their ability to stop leaks; this is the basis for water-stop and sealants. Another use for hydrophobic and hydrophilic materials is the storage and distribution of water and methanol in miniature direct methanol fuel cells (DMFCs). For the distribution of methanol, a material is chosen that wets methanol but is hydrophobic to water. This type of a material is a preferential wicking material. This allows neat methanol to be stored and distributed in a DMFC without water diffusing into the methanol storage layers. The water storage layer at the anode of the fuel cell is a hydrophilic material. It is not preferential to either water or methanol and provides a layer in which they can mix. θ Hydrophilic < π ,θ Hydrophobic ≥ π Example 5.3 A 0.25-mm-diameter tube is placed vertically in a pool of water as shown in Fig. 5.7. The density of water is 1000 kg/m3 and its surface tension is 0.06 N/m. It is assumed that the water can completely wet the tube. Find the capillary rise. 346 Transport Phenomena in Multiphase Systems Solution: Since the water can completely wet the tube, the contact angle is θ = 0. Considering that the density of the vapor is much less than that of the liquid, eq. (5.48) can be simplified as 2σ  ρ A gH r Therefore, the capillary rise is 2σ 2 × 0.06 H= = = 0.098m = 98mm ρA gr 1000 × 9.8 × 0.125 × 10−3 5.4 Phase Equilibrium in Microscale Interfacial Systems 5.4.1 Ultra-Thin Liquid Films, Disjoining Pressure When a liquid film on a solid surface becomes very thin, the intermolecular attractive force between the molecules in the liquid and those in the solid surface turns to pull the liquid back to the liquid film. For a flat liquid film in which capillary pressure is absent, the pressure in the liquid film, pA , is changed by an amount, pd ( pd < 0 ), referred to as disjoining pressure (see Section 2.6.4) which causes the liquid pressure to be less than the vapor pressure ( pA < pv ). At equilibrium, the liquid pressure becomes pA = pv + pd ; in this case, the disjoining pressure is negative when the solid draws the liquid into the film. The disjoining pressure is a product of long-range intermolecular forces composed mainly of molecular (dispersion) and electrostatic interactions. They can be effective as long-range forces when the intermolecular spacing in the liquid film is about 0.2 nm to 10 nm. These forces are present even in nonpolar liquids because they are quantum-mechanical in origin. A pressure gradient is generated within the thin layer of liquid that covers the solid section in contact with the vapor. The disjoining pressure becomes more important for ultra-thin liquid films, because the liquid-vapor interface is closer to the solid surface. As a result, the attractive forces from the solid to liquid will have a greater effect on the liquid-vapor interface. As an example, Fig. 5.10 illustrates the variation of disjoining pressure with the liquid film thickness for CCl4 on glass at 77 ºC. The dependence of disjoining pressure, pd, on liquid film thickness, δ, can be obtained by analyzing the molecular forces and the reversible work required to establish a liquid microlayer (see Section 2.6.4): 2 AH δ 3− n (5.49) pd (δ ) = − π ( n − 2)(n − 3) where AH is a constant, and n is the exponent in the long-range interaction potential φ ( r ) = −Cφ / r n that characterizes the force interaction between the Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 347 Figure 5.10 Disjoining pressure of CCl4 on glass at 77 °C (Potash and Wayner, 1972). liquid and solid molecules. This long-range interaction potential can be represented by the second term of the Lennard-Jones 6-12 potential. This results in the exponent in the long-range potential, n, equaling 6. Thus eq. (5.49) becomes A (5.50) pd (δ ) = − H δ −3 6π A more generalized relationship between the disjoining pressure and the film thickness can be expressed as (Potash and Wayner, 1972) pd = − A′δ − B (5.51) where A′ and B are constants that characterize the molecular and electrostatic interactions. Several components of the disjoining pressure are distinguished: dispersion (molecular), adsorption, electrostatic, and structural (Derjaguin, 1955; 1989). Values of dispersion constants strongly depend on liquid types. For ammonia, for example, A′ = 1 × 10−21 J. The nature of this phenomenon produces increasing (absolute value) pressure with decreasing film thickness. The pressure is considered to be positive for repulsion and negative for attraction of the surface film. Due to the extremely high values of pd in ultra-thin films, the liquid transport can be significant; thus the role of disjoining pressure in evaporation can be essential, particularly for low-temperature fluids. Disjoining pressure is one of the fundamental phenomena that affect the formation of thin evaporating films and the magnitude of the contact angles. When a wetting liquid is in contact with a solid, the liquid forms a curved liquid/vapor interface, as shown in Fig. 5.11. When the liquid has a high wetting capability, the liquid spreads on the solid wall to form an extended meniscus. The extended meniscus can be divided into three regions: the equilibrium thin film, microfilm region, and intrinsic region. In the equilibrium thin film, the liquid molecules adhere to the solid, and the interfacial temperature is equal to the wall 348 Transport Phenomena in Multiphase Systems Non-evaporating equilibrium thin film region (Tw = Tδ) Thin film Extended meniscus Vapor Microfilm region ( Tsat ≤ Tδ ≤ Tw ) Intrinsic region Tw Meniscus region ( Tδ = Tsat ) Liquid Figure 5.11 Thin liquid film on a solid surface. temperature. No evaporation occurs in this region because the liquid/vapor interfacial equilibrium temperature is elevated to the wall temperature due to the disjoining pressure effects. Virtually all of the evaporation occurs in the microfilm region. In this region the disjoining and capillary pressure significantly affect its shape. When the wall temperature is greater than the saturation temperature for a given vapor pressure, the interfacial temperature lies between these two values, i.e., Tsat ( pv ) < Tδ < Tw . In the intrinsic region, the effect of disjoining pressure is negligible and surface tension effect is dominating. The liquid flow that continually feeds the evaporating thin film is driven by a pressure gradient. The pressure gradient near the equilibrium region is primarily caused by changes in disjoining pressure as the film thickens. While the film thickens, the curvature of the surface increases to a maximum. Once the surface curvature is at its maximum, it starts to decrease, and both the disjoining pressure gradient and the change in curvature drive the flow. As the film thickens further, the disjoining pressure effects become negligible and curvature change alone drives the flow. As the film thickness increases further, the evaporation rate drops to zero, and the curvature stays at a constant value. This region of the interface can be referred to as the meniscus, and is the fourth region labeled in Fig. 5.11. The ultra-thin film phenomena and the effects of disjoining pressure are very important in two-phase micro/miniature devices such as micro heat pipes (Khrustalev and Faghri, 1994), rotating miniature heat pipes (Lin and Faghri, 1999), and pulsating heat pipes (Zhang and Faghri, 2002). Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 349 5.4.2 Change in Saturated Vapor Pressure over a Curved Interface In very thin films, as noted above, the attractive force from the solid surface to the liquid produces a pressure difference (disjoining pressure) across the liquidvapor interface, in addition to the capillary effect. These two effects reduce the saturated vapor pressure over a thin film with curvature in comparison with the normal saturated condition. Consider a thin liquid film with liquid thickness δ over a substrate with liquid interface temperature Tδ and normal saturation vapor at constant temperature from the normal saturated pressure psat (Tδ ) to an arbitrary pressure gives pressure corresponding to Tδ of psat (Tδ ) . At equilibrium, the chemical potential in the two phases must be equal (see Chapter 2): μA = μv (5.52) Integrating the Gibbs-Duhem equation, d μ = − sdT + vdp (5.53) μ − μ sat = ³ p psat (T δ ) vdp (5.54) Using the ideal gas law vv = Rg Tδ / Pv ( ) for the vapor phase, and assuming the liquid phase is incompressible (v = vA ) , one obtains the following relations upon integration of eq. (5.54) for the vapor and liquid chemical potentials, respectively: pv ,δ μv ,δ = μ sat ,v + Rg Tδ ln (5.55) psat (Tδ ) μA ,δ = μ sat ,A + vA [ pA − psat (Tδ )] (5.56) Since μ sat ,A = μ sat ,v , substituting eqs. (5.55) and (5.56) into eq. (5.53) yields ­ v [ p − psat (Tδ )] ½ ° ° (5.57) pv ,δ = psat (Tδ ) exp ® A A ¾ Rg Tδ ° ° ¯ ¿ The pressure difference in the vapor phase pv ,δ and the liquid phase pA due to capillary and disjoining effects are related as follows: pv ,δ − pA = pcap − pd eq. (5.57). (5.58) where pcap is capillary pressure. Equation (5.58) can be used to eliminate pA in ­[ pv ,δ − psat (Tδ ) − pcap + pd ] ½ ° ° (5.59) pv ,δ = psat (Tδ )exp ® ¾ ρA Rg Tδ ° ° ¯ ¿ When the interface is flat and pd = 0, pv ,δ = psat (Tδ ) . For a curved interface and pd = 0, eq. (5.59) is reduced to the following Kelvin equation: 350 Transport Phenomena in Multiphase Systems ­[ pv ,δ − psat (Tδ ) − pcap ] ½ ° ° pv ,δ = psat (Tδ ) exp ® ¾ ρA Rg Tδ ° ° ¯ ¿ (5.60) 5.5 Transport Effects at the Interface 5.5.1 Interfacial Mass, Momentum, Energy, and Species Balances The conservation laws for transport phenomena can be reduced to local partial differential equations if they are considered at a point that does not belong to a surface of discontinuity, such as an interface. When considering a discontinuous point, appropriate jump conditions relating the values of the fundamental quantities on both sides of the interface should be considered. Jump conditions at an interface were discussed in Section 3.3.6, but the effects of surface tension and disjoining pressure associated with a nonflat liquid-vapor interface were not taken into account. It is the objective of this subsection to specify mass, momentum, and energy balance at a nonflat liquid-vapor interface, as well as species balance in solid-liquid-vapor interfaces. Mass Balance At a liquid-vapor interface, the mass balance is ′′ mδ = ρA (VA − VI ) ⋅ n = ρv ( Vv − VI ) ⋅ n (5.61) ′′ where mδ is mass flux at the interface due to phase change and VI is the velocity of the interface. For a three-dimensional interface, there are three 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 (5.62) V ⋅ t1 = Vt1 (5.63) V ⋅ t 2 = Vt 2 The interfacial mass balance can be rewritten in these terms by: ′′ mδ = ρA (VA ,n − VI ,n ) = ρv (Vv ,n − VI ,n ) (5.64) (5.65) Momentum Balance For a situation where the effects of surface tension and disjoining pressure are negligible, the momentum balance at the liquid-vapor interface can be described by eq. (3.173), i.e., Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 351 ′ − ′ ) ⋅ n = mδ (VA − Vv ) ′′ (5.66) A v where n is the normal direction of the interface and points toward the vapor. The left-hand side of eq. (5.66) represents the force per unit area acting on the interface, while the right-hand side represents the change of momentum across the interface. For applications involving thin film evaporation and condensation, the effects of surface tension and disjoining pressure will create additional forces on the interface; in these cases, the left-hand side of eq. (5.66) should be modified. The force per unit area created by the surface tension as indicated by eq. (5.10) is σ (T )(1/ RI + 1/ RII )n . The force per unit area created by the disjoining pressure is − pd n. For the situations where the interfacial temperature is not a constant, the contribution of the Marangoni effect, (dσ / dT )∇Tδ , should also be included. Therefore, the momentum balance at the interface becomes §1 1· ( ′A − ′v ) ⋅ n + σ (T ) ¨ + ¸ n − pd n © RI RII ¹ (5.67) § dσ · ′′ −¨ ¸ (∇Tδ ⋅ t )t = mδ (VA − Vv ) © dT ¹ 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. 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: 2 τ ′ = − pI + 2 μ D − μ ( ∇ ⋅ V ) I (5.68) 3 The deformation tensor can be written for a reference frame that is adjusted to the interface: ª ∂Vn 1 § ∂Vn ∂Vt1 · 1 § ∂Vn ∂Vt 2 · º + + « ¨ ¸ ¨ ¸» ¨ ¸ ¨ ¸ ∂xn 2 © ∂xt1 ∂xn ¹ 2 © ∂xt 2 ∂xn ¹ » « « » ∂Vt1 « 1 § ∂Vn ∂Vt1 · 1ª 1 § ∂Vt1 ∂Vt 2 · » T + + D = ∇ V + ( ∇V ) º = « ¨ ¸ ¨ ¸» ¼ ¨ ¸ 2¬ 2 ¨ ∂xt 2 ∂xt1 ¸ » ∂xt1 « 2 © ∂xt1 ∂xn ¹ © ¹ « » ∂Vt 2 « 1 § ∂Vn + ∂Vt 2 · 1 § ∂Vt1 + ∂Vt 2 · » ¸ ¨ ¸ « 2 ¨ ∂x » ¨t ∂xn ¸ 2 ¨ ∂xt 2 ∂xt1 ¸ ∂xt 2 ¹ © ¹ ¬© 2 ¼ (5.69) 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, τ ′ ⋅ n = [ − p 0 0] + ª ∂V 2μ « n « ∂xn ¬ 1 § ∂Vn ∂Vt1 + ¨ 2 ¨ ∂xt1 ∂xn © · ¸ ¸ ¹ 1 § ∂Vn ∂Vt 2 + ¨ 2 ¨ ∂xt 2 ∂xn © ·º 2 (5.70) ¸ » − μ [ ∇ ⋅ V 0 0] ¸» 3 ¹¼ ( 352 Transport Phenomena in Multiphase Systems This can be reduced to the three components to obtain: ∂Vn 2 4 ∂V 2 § ∂Vt ∂Vt 2 · − μ∇ ⋅ V = − p + μ n − μ ¨ 1 + ¸ ∂xn 3 3 ∂xn 3 ¨ ∂xt1 ∂xt 2 ¸ © ¹ (5.71) § ∂Vn ∂Vt1 · τ ′ ⋅ n ⋅ t1 = μ ¨ (5.72) + ¸ ¨ ∂xt ¸ © 1 ∂xn ¹ § ∂V ∂Vt · τ ′ ⋅ n ⋅ t2 = μ ¨ n + 2 ¸ (5.73) ¨ ∂xt ∂xn ¸ ©2 ¹ The momentum equation balance at the interface is then broken into its three components, as follows: τ ′ ⋅ n ⋅ n = − p + 2μ Normal Direction ∂V 4 § ∂V − pA + pv + ¨ μA A ,n − μv v ,n 3 ¨ ∂xA ,n ∂xv ,n © 2 ª § ∂VA ,t1 ∂VA ,t 2 − « μA ¨ + 3 « ¨ ∂xt1 ∂xt 2 © ¬ · ¸ ¸ ¹ ·º ¸» ¸» ¹¼ (5.74) · § ∂Vv ,t1 ∂Vv ,t 2 + ¸ − μv ¨ ¸ ¨ ∂xt ∂xt 2 1 ¹ © §1 1· ′′ +σ ¨ + ¸ − pd = mδ (VA ,n − Vv ,n ) © RI RII ¹ Tangential 1 § ∂V § ∂Vv ,n ∂Vv ,t1 · § dσ ∂VA ,t1 · μ A ¨ A ,n + + ¸ − μv ¨ ¸−¨ ¨ ∂xt ¨ ∂xt ∂xn ¸ ∂xn ¸ © dT 1 1 © ¹ © ¹ Tangential 2 § ∂V ∂VA ,t 2 μ A ¨ A ,n + ¨ ∂xt ∂xn 2 © · § ∂Tδ ¸¨ ¹ ¨ ∂xt1 © · § ∂Tδ ¸¨ ¹ ¨ ∂xt 2 © · ′′ ¸ = mδ VA ,t1 − Vv ,t1 ¸ ¹ (5.75) ( ) · ′′ ¸ = mδ VA ,t 2 − Vv ,t 2 ¸ ¹ (5.76) The non-slip condition at the liquid-vapor interface requires that VA ,t1 = Vv ,t1 · § ∂Vv ,n ∂Vv ,t 2 + ¸ − μv ¨ ¸ ¨ ∂xt ∂xn 2 ¹ © · § dσ ¸−¨ ¸ © dT ¹ ( ) and VA ,t 2 = Vv ,t 2 . The momentum balance at the tangential directions becomes § ∂VA ,n ∂VA ,t1 · § ∂Vv ,n ∂Vv ,t1 · § dσ · § ∂Tδ · + + (5.77) ¸ = μv ¨ ¸+¨ ¸ ¸¨ ¨ ∂xt ¸ ¨ ∂xt ∂xn ¹ ∂xn ¸ © dT ¹ ¨ ∂xt1 ¸ 1 1 © © ¹ © ¹ § ∂V § ∂Vv ,n ∂Vv ,t 2 · § dσ · § ∂Tδ · ∂VA ,t 2 · + (5.78) μ A ¨ A ,n + ¸ = μv ¨ ¸+¨ ¸ ¸¨ ¨ ∂xt ¨ ∂xt ∂xn ¸ ∂xn ¸ © dT ¹ ¨ ∂xt 2 ¸ 2 2 © ¹ © ¹ © ¹ ′′ For most applications, the evaporation or condensation rate – mδ – is not ′′ very high; therefore, it can be assumed that mδ (VA − Vv )  0 . If the liquid and μA ¨ Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 353 vapor phases are further assumed to be inviscid ( equation at the interface can be reduced to §1 1· + pv − pA = σ (T ) ¨ ¸ − pd © RI RII ¹ A = v  0 ), the momentum (5.79) Energy The energy balance at the interface can be obtained from eq. (3.180), i.e., ª§ V2 · § V 2 ·º ′′ q′′ − q′′ ) ⋅ n − (n ⋅ ′ ) ⋅ VA ,rel + (n ⋅ ′ ) ⋅ Vv ,rel = mδ «¨ ev + v , rel ¸ − ¨ eA + A ,rel ¸ » (A v v A 2¸¨ 2 ¸» «¨ ¹© ¹¼ ¬© (5.80) where the work done by surface tension and disjoining pressure is neglected. If the velocity of the reference frame is taken as the interfacial velocity VI, eq. (5.80) can be rewritten as ( kv ∇Tv − kA ∇TA ) ⋅ n − (n ⋅ ′A ) ⋅ (VA − VI ) + (n ⋅ ′v ) ⋅ (Vv − VI ) (5.81) ­ª (V − VI ) 2 º ª (VA − VI ) 2 º ½ ° ° ′′ ® « ev + v = mδ » − «eA + »¾ 2 2 °¬ ¼° ¼¬ ¯ ¿ where the heat fluxes in the liquid and vapor phases have been determined by Fourier’s law of conduction. Equation (5.81) can be rewritten in terms of enthalpy ( h = e + p / ρ ), i.e., ( kv∇Tv − kA ∇TA ) ⋅ n − (n ⋅ ′ ) ⋅ (VA − VI ) + (n ⋅ ′ ) ⋅ (Vv − VI ) A v (5.82) ª º §p p· 1 1 ′′ = mδ « hAv − ¨ v − A ¸ + Vv2 − Vv ⋅ VI − VA2 + VA ⋅ VI » 2 « » © ρv ρA ¹ 2 ¬ ¼ where hAv is the difference between the enthalpy of vapor and liquid at saturation, i.e., the latent heat of vaporization. The stress tensor in eq. (5.82) can be expressed as 2 τ ′ = − pI + 2μ D − μ ( ∇ ⋅ V ) I = − pI + τ (5.83) 3 ′′ Since the relative velocities at the interface satisfy ( VA − VI ) ⋅ n = mδ / ρA and ′′ (Vv − VI ) ⋅ n = mδ / ρv , we have ª §p p ·º ′′ pA ( VA − VI ) ⋅ n − pv ( Vv − VI ) ⋅ n = mδ « − ¨ v − A ¸ » (5.84) « © ρv ρA ¹ » ¬ ¼ Substituting eqs. (5.83) and (5.84) into eq. (5.82), the energy balance can be written as 354 Transport Phenomena in Multiphase Systems ( kv∇Tv − kA ∇TA ) ⋅ n − ( n ⋅τ A )( VA − VI ) + ( n ⋅τ v )( Vv − VI ) (5.85) 1 1 º ′′ ª = mδ « hAv + Vv2 − Vv ⋅ VI − VA2 + VA ⋅ VI » 2 2 ¬ ¼ To simplify the energy equation, the kinetic energy terms are considered negligible and no-slip conditions are assumed at the interface, VA,t = Vv ,t = VI ,t (5.86) Therefore, the energy equation can be rewritten as § ∂Tv ∂T · − kA A ¸ ¨ kv ∂xn ¹ © ∂xn ª 4 ∂V 2 § ∂VA ,t1 ∂VA ,t 2 · ′′ (5.87) = mδ « hAv + ν A A ,n − ν A ¨ + ¸ ¨ ¸ 3 3 © ∂xt1 ∂xn ∂xt 2 ¹ « ¬ 4 ∂V 2 § ∂Vv ,t1 ∂Vv ,t 2 · º − ν v v ,n + ν v ¨ + ¸» 3 3 ¨ ∂xt1 ∂xn ∂xt 2 ¸ » © ¹¼ where ν A 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., (5.88) ( kv∇Tv − kA ∇TA ) ⋅ n = mδ′′hAv Equations (5.61), (5.79), and (5.88) are widely used in the analysis of evaporation and condensation. Species Balance In general, to obtain a detailed solution one needs to solve the continuity, species, momentum, and energy equations in each phase (presented in Chapters 3 and 4) and use the above interfacial balances to couple the conditions in each phase. This provides information for dependent variables, such as pressure, density, mass fraction, velocity, and temperature, as functions of time and space. Such an approach requires detailed numerical simulation in all phases. One may often be interested in a restricted solution for making an order of magnitude analysis or doing an analytical solution for a limiting case, in contrast to a detailed numerical solution. In some applications, dependent variables are also coupled, preventing the solution of one dependent variable independent of the others. In practice, one often must make assumptions to simplify the interface conditions or complexity of the multiphase problem. In other cases one can neglect, based on physical significance, the resistance to mass transfer and/or heat transfer in one of the phases; this permits simplification of the solution techniques and the solution of conservation equations in only one phase. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 355 Some dependent variables, such as temperature, are continuous across the phases; other variables, such as concentration, are discontinuous. 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 ( ) (5.89) 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 (5.90) (5.91) ρ j ,i V j ,i = J j ,i + ω j ,i ρ j V j 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 jk − VI (5.92) ( ) Remembering the overall mass conservation at the interface, m′′ = ρ k ( Vk − VI ) ⋅ n = ρ j V j − VI ⋅ n ( ) (5.93) (5.94) The interfacial species mass flux is: mi′′ = J k ,i ⋅ n + ωk ,i m′′ = J j ,i ⋅ n + ω j ,i m′′ 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 (5.95) 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. (5.94): J j ,i − J k ,i ⋅ n m′′ = (5.96) ωk , i − ω j , i ( ) 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 (5.97) ω j ,i = ω j , i () (ω ) k ,i (5.98) 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 (5.99) 356 Transport Phenomena in Multiphase Systems J j ,1 ⋅ n = − ρ j D j ,12∇ω j ,1 ⋅ n (5.100) 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 ( ∂xA / ∂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. 5.12, where x A,A , x A, s , x A, g are molar fractions in liquid, solid, and gas, respectively, and m′′ is mass flux for A species A. There are two cases shown in Fig. 5.12. 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 in the liquid/solid phases. The mass transfer rate is limited by both phases. Also note that the concentration gradient in the solid [Fig. 5.12(b)] is steeper and does not penetrate as deeply when compared to the liquid concentration gradient [Fig. 5.12(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. 1 I. x A,A I. x A, s II. x A ,A xA Liquid II. x A , s m′′ Gas A I. x A, g II. x A, g y Solid m′′ A Gas 0 y (a) Liquid to Gas I. x A, g II. x A, g (b) Solid to Gas Figure 5.12 Species concentration and mass transfer from solid and or liquid to a gas mixture. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 357 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: (5.101) p A y =0− = x A,A + ( p A, sat ) Liquid to gas y =0 pA y = 0− = x A, s y = 0+ ( pA,sat ) y = 0+ Solid to gas y = 0+ (5.102) is the mole where pA is the partial vapor pressure of A in the gas mixture, x A,A fraction of species A in liquid, x A, s 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,A = x A, s =1 y =0 y =0 and eqs. (5.101) and (5.102) reduce to p A y =0− = p A, sat (5.103) 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 B. 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 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 p A y = 0− x A, g = p (5.104) (5.105) ω A, g = x A, g M A 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. 5.12(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. 5.12(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. ∂x J * = −cDAB A (5.106) A ∂y where J * can be a function of concentration and not necessarily constant. One A application is related to catalytic surface reaction. For this case, m′′ = 0 because 358 Transport Phenomena in Multiphase Systems the mass flux of reactants equals the mass flux of products. Therefore, from eq. (5.94), mi′′ = J k ,i ⋅ n . Catalytic surfaces are used to promote heterogeneous reactions (see Chapter 3), 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. ( y =0 ) n (5.107) b. Mass transfer from gases to liquids or solids There are two major forms of mass transfer from gas to liquid, shown in Fig. 5.13. 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,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. 5.13(a)] p A y = 0− (5.108) x A, A + = y =0 H 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 B.77 provides H for selected aqueous solutions. For binary soluble gases, Henry’s law is not appropriate, and solubility data are usually presented in term 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 Liquid 1 Gas II. x A, g I. x A, g Solid II. x A,A I. x A,A I. x A, s Gas II. x A, g xA I. x A, g II. x A, s y 0 a. Gas to Liquid b. Gas to Solid Figure 5.13 Species concentration and mass transfer from gas mixture to liquid or solid Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 359 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 (5.109) c A, s + = S p A, g − 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 B (Tables B.76, B.78 and B.79). Special care should be exercised when using eq. (5.109) for the variation in form and units – two different versions of S are presented in Appendix B. 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. 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. 5.14). 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, 2004) and some are reproduced in Appendix B. A conventional example for this case is the dissolution of salt into water. 1 I. x A, s Liquid Solution Solid I. x A,A II. x A,A Liquid Solution xA II. x A, s II. x A, s Solid II. x A,A I. x A,A I. x A, s 0 (a) Solid to Liquid (b) Liquid to Solid Figure 5.14 Species concentration and mass transfer from solid to liquid and liquid to solid 360 Transport Phenomena in Multiphase Systems 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. c K p = A, s (5.110) c A, A 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 5.4 Write the continuity, momentum, energy, and species equations and the necessary boundary conditions for the horizontal evaporating capillary tube in Fig. 5.15, which is open at one end and closed at the other. The capillary tube is on the order of 100 m. The evaporation is driven by the concentration gradient of vapor in the air. The evaporation cools the interface while the wall heats the fluid, causing a temperature gradient along the interface. Since the surface tension is a function of temperature, Maragoni stresses are created due to the temperature gradient along the interface. The wall is at a constant temperature, which is the ambient air temperature. The meniscus of a volatile liquid in ambient conditions with no forced heating recedes into the tube because of evaporation. The fluid density and viscosity are constant. Gravitational effects are negligible because the tube lies horizontally and the pore diameter is small enough that free convection can be considered negligible. The disjoining pressure effects are put into the contact angle, because the thin film region of the capillary diameter has minimum effects on the total heat transfer and evaporation. The gas in the vapor region is assumed to be an ideal gas. Assume quasi-steady state and that the specific heats of the vapor and air are the same. Solution: Figure 5.16 shows the control volume, which were chosen due to symmetry, and the boundaries identified with numerical indices of the problem. The incompressible Navier-Stokes equations are solved, and the fluid is considered to be Newtonian. The continuity and momentum equations are ∇⋅V = 0 (5.111) Tw = Tref Axis Vapor/Air Liquid Figure 5.15 Physical model of a horizontal evaporating capillary tube. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 361 6 11 10 8-b 8-a 9 1 3 2 4 5 7 Figure 5.16 Control volumes and boundaries for the horizontal evaporating capillary tube. DV = −∇p + ∇ ⋅ ( μ∇V ) Dt The energy equation in the liquid region is ∂ρ h + V ⋅ ∇ ( ρ h ) = ∇ ⋅ ( k ∇T ) ∂t If the solid wall is being modeled, only the conduction in modeled, since only the steady-state solution is of interest. ∇ 2T = 0 The species equation for vapor phase is ∂ρω + V ⋅ ∇ ( ρω ) = ∇ ⋅ ( ρ D∇ω ) ∂t In the energy equation, the enthalpy, h, is h = c p (T − Tref ) ρ (5.112) (5.113) the solid is (5.114) (5.115) (5.116) This condition is used for both the liquid and the gaseous regions, because the vapor has a specific heat similar to air. The density in the liquid is constant, and the density of the gas is determined from the ideal gas law. pref ρ= (5.117) § ω 1− ω · + ¨ ¸ RuT © M v M air ¹ The diffusion coefficient, D, in the species equation is T 3/ 2 (5.118) D=B p The viscosity and thermal conductivity are considered to be constant in the liquid region, while the mass weighted average of these properties is used in the gaseous region. φ = ωφv + (1 − ω ) φair (5.119) 362 Transport Phenomena in Multiphase Systems At the liquid/vapor interface (boundary 1 in Fig. 5.13), the interfacial boundary conditions are Conservation of Mass ∇ω ⋅ n m′′ = ρ D = ρA ( VA − VI ) ⋅ n = ρ g ( Vg − VI ) ⋅ n 1−ω The subscript g refers to the mixture of vapor and air. Conservation of Normal Momentum mI′′ ( VA − Vg ) ⋅ n + pA − pg = σ ( K1 + K 2 ) Conservation of Tangential Momentum − μA ∇ ( VA ⋅ t ) + μ g ∇ ( Vg ⋅ t ) ⋅ n = γ∇T ⋅ t (5.120) (5.121) (5.122) ( ) where is the surface temperature change with temperature dσ dT . Conservation of Energy −kA ∇T ⋅ n + k g ∇T ⋅ n = hAv m′′ (5.123) Saturation Mass Fraction ω = ωsat No-Slip Velocity in the Tangential Direction VA = Vg Continuity of Temperature TA = Tg (5.124) (5.125) (5.126) The surface tension is a function of temperature and is approximated by the following equation: σ T = σ Tref + γ ⋅ (T − Tref ) (5.127) The mass fraction at the liquid/vapor interface was found using the molar concentration of an ideal gas p xg = sat (5.128) p ωsat = xg M g + (1 − xg ) M air xg M g (5.129) The partial pressure of the vapor is a very crucial parameter in the evaporation process. The function of pressure related to temperature from Yaws (1992) is used. p b log (5.130) =a− 132.953 (T + c ) Since the free surface is always a function of the radial location, f, the two curvatures, K1 and K2, are: Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 363 K1 = − f ′′ ( f ′2 + 1 f′ ) 3 (5.131) 2 (5.132) r f ′2 + 1 The first and second derivates of f should be taken by differentiating the face center locations of the interface. The pressure at the face closest to the wall is calculated with a central differences scheme, but the point of intersection between the wall and the liquid/vapor interface is calculated from the prescribed contact angle, . 1§ a2 · a2 a − ¸ + xa − 2 xb ¨ tan θ © b¹ b (5.133) xw = 2 a 1− 2 b where a is the radial distance from the wall to its adjacent face and b is the radial distance from the wall to the second face from the wall. If the interface is close to the tube mouth, and xw exceeds xm, then the prescribed location of intersection between the interface and the wall to calculate the pressure is xm. The conditions at the inner wall (the boundaries labeled 3, 5, and 8-b in Fig. 5.16) are T = Tref (5.134) V=0 (5.135) ∇ω ⋅ n = 0 (5.136) Far from the tube mouth, ambient air conditions of pressure, temperature, and mass fraction are constants, These conditions apply to a quarter circle with a radius of 6r0, as seen in Fig. 5.16 (boundaries 6 and 7): p = pref (5.137) K2 = − If V ⋅ n ≤ 0 → T = Tref else ∇T ⋅ n = 0 (5.138) (5.139) At the tube liquid entrance (boundary 10) the boundary conditions are: ′′ mI , j A j V ⋅n = ω =0 ¦ j ρπ r02 (5.140) ∇T ⋅ n = 0 (5.141) At the axis of axi-symmetric geometry (boundaries 2, 4, and 8-a), the boundary conditions are ∇(V ⋅ t) ⋅ n = 0 (5.142) V ⋅n = 0 (5.143) (5.144) ∇T ⋅ n = 0 364 Transport Phenomena in Multiphase Systems 1.4 1.3 Distance from tube mouth (mm) 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 10 20 30 40 50 θ= π/24, Rice and Faghri (2005a) = /24, numerical (present) θ= π/12, Rice and Faghri (2005a) = /12, numerical (present) θ= π/6, Rice and Faghri (2005a) = /6, numerical (present) Experimental (Buffone and Experimental (Buffone et al., 2004) Sefiane 2004) time (s) Figure 5.17 Distance of meniscus center distance from mouth vs. time for acetone, for a tube with a diameter of 600 μm (Rice and Faghri, 2005a). Rice and Faghri (2005a) developed a computational liquid/vapor interface tracking technique (see Section 5.7.1) to model an interface between a liquid and a vapor, including mass transfer, for this problem. The effect of contact angle on the distance the meniscus recedes inside the capillary over time is presented in Fig. 5.17. The rate at which the meniscus recedes inside the tube increases with decreasing contact angle. The distance the meniscus recedes over time asymptotically approaches the solution with a contact angle of 0Û. The time it takes the center of the meniscus to be approximately 1.1 mm away from the tube mouth is more than 100% different for the contact angles of /6.0 and /12.0. However the difference reduces to about 6% between contact angles of /12.0 and /24.0. This result verifies the assumption of including all of the solid/liquid/vapor interaction into the contact angle. The evaporation rate increases with decreasing contact angle, because the diffusion length scale decreases with decreasing contact angle. With a contact angle of zero, if the meniscus has a spherical shape, the meniscus spans r0 in the axial direction. The experimental results of Buffone et al. (2004) are also plotted in Fig. 5.17 for comparison. With the smaller contact angles, the numerical work has a slightly higher evaporation rate than Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 365 experimentally observed. The increased evaporation rate can be attributed to the constant temperature boundary condition, whereas the experimental condition is natural convection with the ambient air. Therefore, the numerical work has a slightly higher temperature, and therefore a slightly higher vapor pressure at the interface, leading to a higher evaporation rate. 5.5.2 Interfacial Resistance in Vaporization and Condensation High-heat transfer coefficients, typically associated with evaporation and condensation processes in heat transfer devices, are restricted by interfacial resistance. When condensation occurs at an interface, the flux of vapor molecules into the liquid must exceed the flux of liquid molecules escaping to the vapor phase. When evaporation occurs, on the other hand, the flux of liquid molecules escaping to the vapor phase must exceed the flux of vapor molecules into the liquid. Schrage (1953) used the kinetic theory of gases to describe condensation and evaporation processes and considered separately the fluxes of condensing and vaporizing molecules in each direction. From kinetic theory, the mass flow rate (of molecules) passing in either direction (right or left) through an imagined plane is given by § Mv · p mδ = ¨ (5.145) ¸ 1/ 2 © 2π Ru ¹ T ′′ where mδ is the flux of molecules, Ru is the universal gas constant, and Mv is the molecular mass of the vapor. The net molecular flux through an interface is ′′ ′′ ′′ (5.146) mδ = mδ + − mδ − Actually, the Boltzmann transport equation should be solved with appropriate boundary conditions of thermal equilibrium at several mean free path distances from the interface. However, a reasonable approximation can be obtained by means of correction factors if it is assumed that the interaction between molecules leaving the interface and those approaching was at equilibrium. This obtains the following relation for the net mass flux at the interface: q′′ M v § Γpv p· ′′ − A¸ mδ = δ = α (5.147) ¨ 2π Ru ¨ Tv hAv TA ¸ © ¹ where α is the accommodation coefficient ( α ≤ 1 ), and the function Γ is given by Schrage (1953): Γ ( a ) = exp a 2 + a π ª1 + erf ( a ) ¼ º (5.148) ¬ 1/ 2 () Γ ( − a ) = exp ( a ) − a 2 π ¬1 − erf ( a ) ¼ ª º (5.149) 366 Transport Phenomena in Multiphase Systems where a= and erf ( a ) = ′′ qδ ρv hAv Mv 2 RuTv a (5.150) π³ 2 0 e − x dx 2 (5.151) is the Gaussian error function. The heat flux to the interface is equal to the net mass flux multiplied by the ′′ ′′ ′′ latent heat ( qδ = mδ hAv ) . Since Γ is a function of qδ , eq. (5.147) does not provide an explicit relation for the interfacial heat flux. Assuming that pA and pv are the saturation pressures corresponding to TA and Tv, eq. (5.147) can be represented in the following form: M v ª Γpsat (Tv ) psat (TA ) º ′′ − qδ = α hAv (5.152) « » 2π Ru « Tv TA » ¬ ¼ For evaporation and condensation of working fluids at moderate and high temperatures, a is usually very small according to its definition [see eq. (5.150)]. In such a case, eq. (5.146) can be approximated by Γ =1+ a π (5.153) ′′ ′′ An explicit relation for qδ and mδ was obtained by Silver and Simpson (1961) by substituting eq. (5.153) into eq. (5.147) and using ρv = pv M v / RuTv : § pv p· (5.154) − A¸ ¨ ¨T TA ¸ ©v ¹ which is referred to as the Kucherov-Rikenglaz equation (Kucherov and Rikenglaz, 1960) in the Soviet literature. One can develop an alternative form of eq. (5.152) for small a by assuming that ( pv − pA ) / pv << 1 , (Tv − TA ) / Tv << 1 , and by using the Clausius-Clapeyron relation. Mv § pv vAv · 2α · § hA2v · ′′ § qδ = ¨ (5.155) ¸ ¨1 − ¸ (Tv − TA ) ¸¨ 2hAv ¹ © 2 − α ¹ © Tv vAv ¹ 2π RuTv © ′′ mδ = ′′ qδ § 2α · M v =¨ ¸ hAv © 2 − α ¹ 2π Ru Using the above relation, the heat transfer coefficient at the interface hδ is obtained from the following equation: 2 ′′ qδ Mv § pv vAv · § 2α · § hAv · hδ = =¨ (5.156) ¸ ¨1 − ¸ ¸¨ (Tv − TA ) © 2 − α ¹ © Tv vAv ¹ 2π RuTv © 2hAv ¹ If hδ is of the same order of magnitude as the other h values, the effects of interfacial resistances should be accounted for. It should also be noted that the ′′ above equations relating qδ , h δ , and (Tv − TA ) apply equally well to Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 367 ′′ evaporation and to condensation, with the convention that qδ is positive for condensation and negative for evaporation. It is clear that predicting the interfacial resistances using of any of the above equations depends on the value of the accommodation coefficient α , which varies widely in the literature. Paul (1962) compiled the accommodation coefficients for evaporation for a large number of working fluids. Mills (1965) recommended that α should be less than unity when the working fluids or the interface is contaminated. 5.5.3 Formation of and Heat Transfer through Thin Liquid Films The next two sections describe the heat and mass transfer analysis of an axiallygrooved heat pipe (AGHP) to illustrate an application of the ideas presented thus far. The operating principles of the heat pipe are outlined in Section 1.6.4. A partial cross-section of the particular axially-grooved heat pipe configuration considered here is shown in Figures 5.18 (a) and (b) and 5.19. Slots with walls of angle γ relative to the radial direction, and of mouth width W, are cut along the inner surface length (z-axis) of a circular pipe; this forms a series of fins through which heat is transferred between the working fluid and the heat pipe structure. For this analysis, the heat transfer processes in the heat pipe container and the working fluid are considered to be one-dimensional in the radial direction; as a result, axial heat conduction (into or out of the page in Fig. 5.18) is neglected. The radial direction is represented by the x-axis in Fig 5.18. The working fluid forms a meniscus of angle men with the slot wall as shown in Fig. 5.18(a), and the fluid may form a thin film of thickness δ over the fin top surface, particularly during the condensing portion of the heat pipe cycle. Since the thermal resistance of a low-temperature AGHP depends primarily on the film thickness in the condenser and evaporator sections, evaluation of the thin liquid film heat transfer will be the focus of this analysis. Moreover, since the heat transfer and fluid dynamic processes in a thin film are similar in both the evaporator and condenser sections of the heat pipe, the formation of the films may be described by the same equations as long as one takes into account the different directions of the temperature potential. During the condensation process, liquid in the subcooled thin film flows toward the meniscus region along the s-coordinate. During evaporation, liquid in the superheated thin film flows from the meniscus region in the opposite direction. The model presented here assumes a given meniscus contact angle and accounts for the effects of interfacial thermal resistance, disjoining pressure, and surface roughness. In the evaporating section of the heat pipe, the thin liquid film flows over the crown of the fin and undergoes evaporation caused by heat transfer from the heat-loaded fin surface, which has the curvature Kw. The local heat flux through the liquid film due to heat conduction is (see Fig. 5.19) 368 Transport Phenomena in Multiphase Systems Figure 5.18 Cross-section of the characteristic element of (a) an axially-grooved condenser, and (b) an evaporator (Khrustalev and Faghri, 1995; A and B are shown in Fig. 5.19). q′′ = kA Tw − Tδ where the local thickness of the liquid layer δ and the temperature of the free liquid film surface T are functions of the s-coordinate. For small Reynolds δ (5.157) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 369 numbers (less than 10), it is valid to assume a fully-developed laminar liquid flow velocity profile (Khrustalev and Faghri, 1995): 1 dpA uA = − 2ηδ − η 2 (5.158) 2 μA ds where η is the coordinate normal to the solid-liquid interface. In arriving at eq. (5.158), it is assumed that the shear stresses at the interface are negligible. It is further assumed that the vapor pressure is constant along the s-coordinate, and that the liquid flow is driven mainly by the surface tension and the adhesion forces (see Problem 5.10). §1 dpA dK dpd dσ dTδ d 1· 2 = −σ + −K + (5.159) ρ v2 vv ,δ ¨ − ¸ ds ds ds dTδ ds ds © ρv ρA ¹ ( ) ( ) 1 dª 3d º kA (Tw − Tδ ) (5.161) «δ ds ( pd − σ K ) » = h ρ δ 3μA ds ¬ ¼ Av A The film surface curvature K expressed in terms of the solid surface curvature Kw and film thickness is d 2δ K = Kw + 2 ds −3 / 2 where K is the local interface curvature and pd is the disjoining pressure. The impact of the last two terms on the results was found to be negligible in the present analysis, and therefore they are omitted in the following equations. The continuity and energy equations for the evaporating liquid layer lead to (see Problem 5.11) dδ q′′ (5.160) uA dη = ds 0 hAv ρA Substituting eqs. (5.157) – (5.159) into eq. (5.160) gives the following relation for the thickness of the evaporating film, δ ( s ) : ³ ª § dδ ·2 º (5.162) «1 + ¨ ¸» « © ds ¹ » ¬ ¼ It is assumed that the absolute value of the vapor core pressure at any z-location along the groove is related to vapor temperature by the saturation conditions pv = psat (Tv ) (5.163) and therefore can be defined for a given Tv using the saturation tables. The interfacial temperature T , which is affected by the disjoining and capillary pressure, also depends on the value of the interfacial resistance. For the case of a comparatively small heat flux, interfacial resistance is defined by the following relation [see eq. (5.154)]: § 2α · hAv ª pv ( psat )δ º ′′ qδ = − ¨ (5.164) − « » ¸ Tδ » © 2 − α ¹ 2π Rg « Tv ¬ ¼ where pv and ( psat )δ are the saturation pressures corresponding to Tv (assumed to be a constant along s) in the bulk vapor and at the thin liquid film interface (varies with s), respectively. 370 Transport Phenomena in Multiphase Systems The relation between the vapor pressure over the thin evaporating film, ( psat )δ –which is affected by the disjoining pressure – and the saturation pressure corresponding to Tδ , psat (Tδ ) , is given by eq. (5.59): ª ( p ) − psat (Tδ ) + pd − σ K º (5.165) ( psat )δ = psat (Tδ ) exp « sat δ » ρA Rg Tδ « » ¬ ¼ This reflects the fact that under the influence of disjoining and capillary pressures, the liquid free surface saturation pressure ( psat )δ differs from the normal saturation pressure psat (Tδ ) and varies along the thin film (i.e., scoordinate); pv and Tv are, however the same for any value of s at a given zlocation. ( psat )δ also varies due to the fact that Tδ changes along s. Under steady-state conditions, the heat flux obtained by eqs. (5.157) and (5.164) are the same. Combining these two expressions yields δ § 2α · hAv ª pv ( psat )δ º − (5.166) Tδ = Tw + ¨ « » ¸ kA © 2 − α ¹ 2π Rg « Tv Tδ » ¬ ¼ Equations (5.165) and (5.166) determine the interfacial temperature Tδ and pressure psat (Tδ ) . The wall temperature Tw , which results from solving heat conduction between grooves (i.e., heat conduction by the fins), must be provided as an input to the solution procedure. As the liquid film thins, the disjoining pressure pd and the interfacial temperature T increase. A nonevaporating film thickness is present under specific conditions, as demonstrated in Section 5.4.1, that gives the equality of the liquid-vapor interface and solid surface temperatures: Tδ = Tw (5.167) Substituting eq. (5.167) into eq. (5.166), one obtains T (5.168) ( psat )δ = pv w Tv Substituting eq. (5.168) into eq. (5.165), one obtains the disjoining pressure for the nonevaporating film thickness as § T pv Tw · +σ K pd = − pv w + psat (Tw ) + ρA Rg Tw ln ¨ (5.169) ¨ p (T ) T ¸ ¸ Tv v¹ © sat w At the nonevaporating film thickness, the disjoining pressure can also be obtained by eq. (5.51), i.e., pd = − A′δ 0− B (5.170) Combining eqs. (5.169) and (5.170) yields the nonevaporating film thickness: ­1 ª º½ § T pv Tw · ° ° δ 0 = ® « pv w − psat (Tw ) − ρA Rg Tw ln ¨ ¸ − σ K »¾ ¨ p (T ) T ¸ Tv »° v¹ ° A' « © sat w ¼¿ ¯¬ −1/ B (5.171) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 371 which is applicable to nonpolar liquids. For water, however, the following logarithmic dependence of disjoining pressure on the liquid film thickness is preferable (Holm and Goplen, 1979): ª § δ ·b º pd = ρA Rg Tδ ln « a ¨ (5.172) ¸» « © 3.3 ¹ » ¬ ¼ where a = 1.5336 and b = 0.0243. Example 5.5 Derive the expression for the equilibrium film thickness of water on a wall. If the wall temperature Tw = 101 °C and the saturated vapor temperature Tv = 100 °C, find the non-evaporating film thickness 0. The density of the liquid water can be assumed to be ρA = 1000kg/m3 . Solution: The disjoining pressure for water at nonevaporating film thickness is expressed by eq. (5.172). The nonevaporating film thickness δ 0 can be obtained by rearranging eq. (5.172), i.e., ª1 § pd · º δ 0 = 3.3 « exp ¨ » ¨ ρA Rg Tδ b ¸ » ¸ «a © ¹¼ ¬ Substituting eq. (5.169) into the above equation, one can obtain the non-evaporating film thickness for water as follows: ­1 ª p (T ) − pv Tw / Tv + σ K § pv Tw · º ½ ° ° + ln ¨ δ 0 = 3.3 ® exp « sat w ¨ p (T ) T ¸ » ¾ ¸ ρA Rg Tw « v ¹» ° °a © sat w ¬ ¼¿ ¯ (5.173) Since the molecular weight of water is 18, the gas constant for the water vapor is Rg = 8314 /18 = 461.9J/kgK . The wall and vapor temperatures are respectively Tw = 101 + 273 = 374K and Tv = 100 + 273 = 373 K. The curvature of the liquid film is assumed to be negligible, i.e., K = 0. The vapor pressure is the saturation pressure corresponding to the vapor temperature, i.e., pv = psat (Tv ) = 1.013 × 105 Pa The saturation pressure corresponding to the wall temperature can be obtained from the Table B.48 as psat (Tw ) = 1.0547 × 105 Pa. Therefore, the nonevaporating film thickness is ­1 ª p (T ) − pv Tw / Tv + σ K § pv Tw · º ½ ° ° δ 0 = 3.3 ® exp « sat w + ln ¨ ¨ p (T ) T ¸ » ¾ ¸ ρA Rg Tw « v ¹» ° °a © sat w ¬ ¼¿ ¯ 1/ b 1/ b 1/ b 372 Transport Phenomena in Multiphase Systems ­1 ª1.0547 × 105 − 1.013 × 105 × 374 / 373 ° = 3.3 × ® exp « 1000 × 461.9 × 374 °1.5336 ¬ ¯ § 1.013 × 105 374 · º ½ ° + ln ¨ = 1.51 × 10−8 m = 0.0151 m ¨ 1.0547 × 105 373 ¸ » ¾ ¸ © ¹» ° ¼¿ It can be seen that the non-evaporating film is extremely thin. 1/ 0.0243 5.5.4 Heat Transfer in the Thin-Film Region of an Axially-Grooved Structure Khrustalev and Faghri (1995) modeled an evaporating extended meniscus in a capillary groove, shown in Fig. 5.18(b). Most of the heat is transferred through the region where the liquid layer is extremely thin. The significance of the temperature difference between the saturated vapor core and the interface has been stressed in the model. Manufacturing processes always leave some degree of roughness on metallic surfaces. Alloys of copper, brass, steel, and aluminum invariably have some distinct grain structure resulting from processing. Corrosion and deposition of some substances on the surface can further influence its microrelief. This means that the solid surface of a working material is totally covered with microroughnesses, where the characteristic size Rr may vary from 10-8 to 10-6 m, for example. Apparently, thin liquid film formation can be affected by some of these microroughnesses. It can be assumed that at least some part of a single roughness fragment on which thin film formation takes place has a circular crosssection and is extended in the z-direction due to manufacturing of the axial grooves (see Fig. 5.19). As can be seen in Fig. 5.19, the free liquid surface is divided into four regions. The first region is the equilibrium nonevaporating film. The second (microfilm) region ranges in the interval δ 0 ≤ δ ≤ δ1 , where the liquid film thickness increases up to the value δ1 as described by eqs. (5.161) and (5.162). In this region, the generalized capillary pressure pcap ≡ σ K − pd (5.174) (defined so that its value is positive) changes drastically along the s-coordinate. The capillary pressure changes from the initial value up to an almost-constant value at point s1, where the film thickness δ1 is large enough to warrant neglecting the capillary pressure gradient. It is useful to mention that some investigators have denoted this microfilm region as the interline region. The third (transition) region, where the liquid-vapor interface curvature is constant, is bounded by δ1 < δ < Rr + δ 0 , and the local film thickness is determined by the geometry of the solid surface relief and the value of the meniscus radius Rmen. In the fourth (meniscus) region, where δ > Rr + δ 0 by definition, the local film Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 373 Figure 5.19 Thin evaporating film on a fragment of the rough solid surface (Khrustalev and Faghri, 1995; A and B corresponding to points shown on Fig. 5.18). thickness can be considered independent of the solid surface microrelief. In the third and fourth regions, heat transfer is determined by heat conduction in the meniscus liquid film and the metallic fin between the grooves. In the second region, however, the temperature gradient in the solid body can be neglected in comparison to that in the liquid due to the extremely small size of this region. The total heat flow rate per unit groove length in the microfilm region is defined as s1 T − T s1 w δ ′ qmic ( s1 ) = ds ≡ q′′ds (5.175) 0 δ /k 0 A Equations (5.161) – (5.163), (5.165) and (5.166) must be solved for four ′ variables: δ , d δ / ds , pcap and qmic ( s ) in the interval from s = 0 to the point s = s1, where pcap can be considered to be constant. Instead of considering the two ³ ³ 374 Transport Phenomena in Multiphase Systems second-order equations (5.161) and (5.162), the following four first-order equations should be considered with their respective boundary conditions: dδ =δ ' (5.176) ds −B 3/ 2 § pcap − A ' δ dδ ' 1· = 1 + δ '2 +¸ (5.177) ¨ ¨ ds σ Rr ¸ © ¹ dpcap 3ν ′ = − A 3 qmic ( s ) (5.178) ds hAvδ ′ dqmic Tw − Tδ = (5.179) δ / kA ds ( ) δ s =0 = δ0 (5.180) δ ' s =0 = 0 pcap s =0 (5.181) + A 'δ 0− B =− σ Rr + δ 0 s =0 (5.182) (5.183) ′ qmic =0 where the value of δ 0 can be found from eq. (5.171) with K = −1 / Rr . Though the initial-value problem described by eqs. (5.176) – (5.183) is completely determined, its solution must satisfy one more condition – namely, (5.184) Rmen Since the only parameter in this problem that is not fixed is connected with the surface roughness characteristics, the boundary condition (5.184) can be satisfied by the choice of Rr. Physically, this means that the beginning of the evaporating film shifts along the rough surface depending on the situation so as to satisfy the conservation laws. However, in a smooth surface model ( Rr → ∞ ) the solution will probably not satisfy eq. (5.184). As a result of this problem, the ′ values of δ1 and qmic ( s1 ) can be determined and the transition region can be considered provided that δ1 < Rr , where the free liquid surface curvature is constant and its radius Rmen is many times larger than Rr. Based on the geometry shown in Fig. 5.18, s = s1 pcap = σ where θ f is the contact angle in the microfilm region. Equation (5.185) is valid for the rough surface model ( θ f can be set equal to zero for very small Rr) as well as the smooth surface model in the meniscus region ( Rr → ∞ and θ f is given as a result of the microfilm problem solution). The heat flow rate per unit groove length in the transition region is δ = δ 0 + Rr − Rr2 − x 2 − Rmen + ( Rmen 2 + x 2 + 2 Rmen x sin θ f ) 1/ 2 (5.185) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 375 Tw − Tδ dx (5.186) δ / kA where xf and xtr can be obtained from eq. (5.185), provided δ = δ1 and δ = Rr + δ 0 , respectively. At this point, the point of connection between the transition and meniscus regions must be considered. Here the film thickness, free surface curvature, and liquid surface slope angle must coincide on both sides. In the rough surface model, the last condition is always satisfied because the length of the microfilm region is smaller than Rr and the rough fragment with the film can be “turned” around its center in the needed direction. In other words, because of the circular geometry of the rough fragment and the constant temperature of the solid surface in the microfilm region, the slope of the film free surface is not fixed in the mathematical model. On the contrary, in the smooth surface model the numerical results give θ f , which is generally not equal to θ men as determined by the fluid ′ qtr = ³ xtr xf flow along the groove. Note that in a situation where θ f ≠ θ men , the smooth surface model can be used along with the rounded fin corner, where the radius is Rfin. In this type of case eq. (5.185) can also be used, provided Rr is changed to Rfin. It is useful to note that the values of Rmen and men are connected by the geometric relation θ men = arccos (W / 2 Rmen ) − γ . These values should be given as a result of the solution of the problem for fluid transport along the groove. The fin top temperature Tw should be defined based on consideration of the heat conduction problem in the fin between the grooves and in the meniscus liquid film, as discussed below. The free liquid surface curvature K in the microfilm region varies from the initial value to that in the meniscus region. Its variation is described by eqs. (5.176) – (5.184) with respect to the pcap and pd definitions. In spite of a sharp maximum, which the K function has in the microfilm region, its variation only slightly affects the total heat transfer coefficient. To check this hypothesis numerically, a simplified version of the heat transfer model of the microfilm region was developed by Khrustalev and Faghri (1995), wherein it was assumed that the microfilm free surface curvature is equal to that in the meniscus region. Therefore, instead of solving eqs. (5.176) – (5.184), the microfilm thickness in this region (and in the transition region) can be given by eq. (5.185) for the interval 0 ≤ x ≤ xtr . In this case, the heat flow rate per unit groove length in both the microfilm and transition regions is xtr T − T w δ ′ ′ qmic + qtr = dx (5.187) 0 δ / kA Khrustalev and Faghri (1995) solved eqs. (5.165) and (5.166) simultaneously for Tδ (absolute error Δ a = 0.0001 K) and ( psat )δ ( Δ a = 1 Pa) for every point ³ on s. The system of the four first-order ordinary differential equations with four 376 Transport Phenomena in Multiphase Systems initial conditions and one constitutive condition describing the evaporating microfilm region, eqs. (5.176) – (5.184), were solved using the fourth-order Runge-Kutta method and the shooting method (on parameter Rr). The controlled relative error was less than 0.001 for each of the variables. The results obtained for comparatively small temperature drops through the thin film were compared with those from the simplified model. The results presented in this section mostly refer to the AGHP with a total length of Lt = 0.914 m, a condenser length of Lc = 0.152 m, and an evaporator length of 0.15 ≤ Le ≤ 0.343 m; other geometric parameters are: W = 0.61 mm, tg Figure 5.20 Characteristics of the evaporating film along the solid-liquid interface (ammonia, Tv=250 K): (a) Free liquid surface temperature; (b) Thickness of film; (c) Generalized capillary pressure (Khrustalev and Faghri, 1995). Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 377 = 1.02 mm, L1 = 0.43 mm, Rv = 4.43 mm, Ro = 7.95 mm, γ = 0o , and N = 27. The working fluids were ammonia and ethane, and the thermal conductivity of the casing material was assumed to be kw = 170 W/m-K, with dispersion constants A′ = 10−21 J and B = 3. The data in Figs. 5.20 and 5.21 were obtained for ammonia in the evaporator with a vapor temperature of Tv = 250 K and α = 1 . The solid surface superheat was ΔT = Tw − Tv , and results obtained using the simplified model for evaporating film are denoted as SIMPL. Figure 5.20(a) shows the variation of the free liquid surface temperature along the evaporating film for ΔT = 0.047 K, 0.070 K, and 0.120 K, which values are from the solutions of eqs. (5.51), (5.165), (5.166) and (5.176) – (5.184) in the microfilm region. These results are compared to those obtained by the simplified model, which is based upon the assumption that the curvature of the microfilm free surface is equal to that of the meniscus region. In the simplified model for the case of a smooth surface, the value of the contact angle in the microfilm region was θ f = 7 o ; this value was obtained by the numerical solution of eqs. (5.176) – (5.184). The corresponding variations in the film thickness δ and generalized capillary pressure pcap are shown in Figs. 5.20 (b) and (c). The results obtained by the simplified model have been artificially shifted along the s-coordinate in these figures (and also in Fig. 5.19) to make the comparison more understandable. Also, it should be noted that there is some difference between the s-coordinate and the x-coordinate used in the simplified model. The following relation has been used in the present analysis: s = Rr arcsin( x / Rr ). In Fig. 5.20(a), the interval of Tδ variation along the evaporating film from the value of Tw to approximately Tv is very prolonged, and the interfacial thermal resistance is significant even when the film thickness is larger than 0.1 m. For a smaller characteristic size Rr, the film thickness increased more sharply along the solid surface, as shown in Fig. 5.20(b), which is in agreement with eq. (5.185). It should be mentioned that unlike for the simplified model, the problem defined by eqs. (5.176) – (5.184), Rr is not a parameter but the result of the numerical solution. The values of the maximum heat flux in the microfilm region were extremely high in comparison to those in the meniscus region (see Fig. 5.21). For ΔT = 0.120 K, the generalized capillary pressure pcap decreased from the initial value to an almost-constant value by approximately 5000 times [see Fig. 5.20(c)]. For a larger ΔT , this sharp decrease can cause some difficulties in the numerical treatment while solving eqs. (5.176) – (5.184); that is why the simplified model is useful. The simplified model has given the variation of pcap along the film, which is even more drastic because the surface tension term is absent in the capillary pressure gradient [see Fig. 5.20(c)]. However, this assumption caused a comparatively small decrease in total heat flow rate in the microfilm region, as illustrated by Fig. 5.18(a). The distribution of the heat flux in the microfilm, transition and beginning of meniscus regions for different 378 Transport Phenomena in Multiphase Systems Figure 5.21 Heat flux through the evaporating film (ammonia, Tv=250 K, =1): (a) along the solid-liquid interface (microfilm region); (b) along the fin axis (Rr=33 m, T = 1 K) (Khrustalev and Faghri, 1995). meniscus contact angles θ men – as predicted by the simplified model – are presented in Fig. 5.21(b). The total heat flow through the meniscus region was significantly larger than that through the microfilm region. This means that the simplified model should provide the needed accuracy when estimating the heat transfer coefficient for an evaporator element, as shown in Fig. 5.18. Khrustalev and Faghri (1995) showed that the simplified model underestimated the overall heat transfer coefficient by only 5%, which permits its use as a means of avoiding the numerical difficulties described above. Example 5.6 A schematic of a pulsating heat pipe (PHP) with one end sealed and the other end open is shown in Figs. 5.22 and 5.23 (Zhang and Faghri, 2002). The evaporator section of the heat pipe is near the closed end, and its length is Lh. The condenser section with a length of Lc is near the open end of the heat pipe, and the adiabatic section with length La is located between the evaporator and condenser sections. When the evaporator is heated, the vapor pressure increases. The liquid plug moves toward the open end because the vapor pressure, pv, is higher than the environment pressure, pe. As a result, the volume of the vapor slug is increased and part of the vapor slug is exposed to the cooled section of the heat pipe, where vapor is condensed to liquid. When the rate of condensation Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 379 exceeds the rate of evaporation, the vapor pressure will decrease. When the vapor pressure, pv, is decreased to a value below the environment pressure, pe, the liquid plug will be pushed back to the closed end. At this point, the rate of evaporation again exceeds the rate of condensation, which enables the vapor pressure to increase and pushes the liquid plug to the open end. This process is repeated, and the oscillation of the liquid plug can be maintained. As the liquid slug moves toward the open end of the pulsating heat pipe, the trailing edge of the liquid plug leaves a thin liquid film on the pipe wall. Evaporation and condensation over this thin liquid film are the driving forces of pulsation flow in a PHP with an open end. Obtain the equation that governs the liquid film thickness in 0 < x < L1 . Solution: Figure 5.23 shows the physical model of the evaporator section of a PHP with an open end. The continuity equation of the liquid film is δ 1§ q· uA dη = (5.188) ¨ mA − ¸ 0 ′¹ hAv ,e ¸ 2π R ρ A ¨ © ³ Figure 5.22 Pulsating heat pipe with an open end. Figure 5.23 Heat transfer model for evaporator. 380 Transport Phenomena in Multiphase Systems where q is the rate of heat through a given cross section due to phase change and is defined as follows s k (T − T ) c Ah q = 2π R ds (5.189) 0 δ The revised latent heat of evaporation is defined as ′ hAv ,e = hAv + 0.68c pA (Th − Tv ) (5.190) ³ to account for the contribution of the sensible heat. For small Reynolds numbers (less than unity), a fully-developed laminar liquid velocity profile can be assumed. 1 dpA uA = − (5.191) 2ηδ − η 2 2 μA dx Substituting eq. (5.191) into eq. (5.188), the continuity equation in the liquid film becomes dp 3μ A § q· (5.192) − A= ¨ mA − ¸ ′¹ dx 2π R ρAδ 3 ¨ hAv ,e ¸ © The pressures in the vapor and liquid phases have the following relationship: pv − pA = σ K − pd (5.193) where the curvature is ( ) ª § dδ ·2 º ª 1 § dδ · º − (5.194) cos «arctan ¨ «1 + ¨ ¸» ¸» R −δ © dx ¹ ¼ ¬ « © dx ¹ » ¬ ¼ and the disjoining pressure is ª § δ ·b º pd = ρA Rg Tv ln « a ¨ (5.195) ¸» « © 3.3 ¹ » ¬ ¼ where a = 1.5336 and b = 0.0243. It is assumed that the pressure in the vapor phase is constant. By differentiating eq. (5.193), we have dp d − A = (σ K − pd ) (5.196) dx dx Combining eqs. (5.192) and (5.196), one obtains 3μ A § d q· (5.197) (σ K − pd ) = m− ¸ 3¨ A ′¹ dx hAv ,e ¸ 2π R ρAδ ¨ © The boundary conditions of eq. (5.197) are dδ = 0, x = 0 (5.198) dx d 2δ K= 2 dx −3 / 2 Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 381 ­1 , x p < Lh d 2δ ° x=0 (5.199) = ® R − δ tr 2 dx °0, x p < Lh ¯ δ = δ 0 , x = L1 (5.200) dδ = tan α , x = L1 (5.201) dx where δ 0 is nonevaporating liquid film thickness. The liquid film thickness can be found by numerical solution or by an approximate solution of eq. (5.197) (Zhang and Faghri, 2002). 5.6 Dynamic Behaviors of Interfaces The instability of a horizontal or vertical interface between a liquid and vapor at different velocities will be discussed in Section 5.6.1, followed by a discussion of waves at liquid-vapor interface in Section 5.6.2. 5.6.1 Rayleigh-Taylor and Kelvin-Helmholtz Instabilities A given physical state is said to be stable if it can withstand a disturbance and still return to its original state. Otherwise, the particular state is unstable. The objective of stability analysis is to analyze the effect of a particular disturbance on the physical state. If φ (it can be velocity, pressure, or temperature) represents a basic solution, a disturbance φ ′ is added to this basic solution and φ + φ ′ will be substituted into the governing equations. The governing equations with φ as a dependent variable are then subtracted from the governing equation y, v uA g Liquid x, u δ ( x, t ) uv Vapor Figure 5.24 Rayleigh-Taylor instability for dense liquid overlay less dense vapor. 382 Transport Phenomena in Multiphase Systems Liquid Vapor δ ( x, t ) g uA x, u uv y, v Figure 5.25 Kelvin-Helmholtz instability in vertical liquid-vapor interface for concurrent flow. with φ + φ ′ as dependent variables to yield disturbance equation for φ ′ . If the disturbance φ ′ damps out, φ is stable; otherwise, if the disturbance φ ′ grows with increasing time, φ is unstable. Interface morphology is important for heat and mass transfer. In a horizontal co-current flow system where a dense liquid phase overlays a less dense vapor phase (see Fig. 5.24), both phases are incompressible, inviscid, and immiscible; the interface may become unstable if there is a disturbance δ ( x, t ) . This instability is referred to as Rayleigh-Taylor instability. On the other hand, if the gravity is parallel to the directions of the liquid-vapor co-current flow (see Fig. 5.25), the instability is referred to as the Kelvin-Helmholtz instability. The conditions for which the interfaces are stable with respect to an arbitrary perturbation δ ( x, t ) will be presented in this subsection (Carey, 1992). Assuming that the liquid and vapor flows are two-dimensional, the governing equations for configuration in Fig. 5.24 are ∂u ∂v + =0 (5.202) ∂x ∂y ρ« ρ« ª ∂u ∂u ∂u º ∂p +u +v » =− ∂x ∂y ¼ ∂x ¬ ∂t (5.203) ª ∂v ∂v ∂v º ∂p + u + v » = − − ρg (5.204) ∂x ∂y ¼ ∂y ¬ ∂t The velocities and the pressure are decomposed as follows into base flow and perturbed components: Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 383 u = u + u ′ , v = v + v′ p = p + p ′ (5.205) Substituting eq. (5.205) into eqs. (5.202) – (5.204) and considering the fact that the base flow should also satisfy eqs. (5.202) – (5.204), the equations for the flow simplify to ∂u ′ ∂v′ + =0 (5.206) ∂x ∂y ρ« ª ∂u ′ ∂p′ § ∂u ′ · º +u¨ ¸» = − ∂x © ∂x ¹ ¼ ¬ ∂t (5.207) ª ∂v′ ∂p′ § ∂v′ · º +u¨ (5.208) ¸» = − ∂y © ∂x ¹ ¼ ¬ ∂t In arriving at eqs. (5.206) – (5.208), the products of perturbation (primed) terms are neglected, and we recognize that ∂u / ∂x = ∂u / ∂y = v = 0 . Differentiating eqs. (5.207) and (5.208) with respect to x and y, respectively, then summing them and substituting the continuity equation, yields the Laplace equation for the pressure perturbation field: ∂ 2 p′ ∂ 2 p′ + 2 =0 (5.209) ∂x 2 ∂y The shape of the interface at time t can be described by δ ( x, t ) = Aeiα z + β t (5.210) and the perturbation quantities v′ and p′ can be postulated to have the following forms: ˆ v′( x, y, t ) = veiα x + β t (5.211) ρ« ˆ p′( x, y, t ) = peiα x + β t (5.212) ˆ ˆ where v and p are the magnitudes of the perturbation. Employing the Young-Laplace equation and the equation for the curvature of the liquid film, pv − pA = σ / R, Carey (1992) used the perturbation analysis to obtain the following condition for an unstable interface ªσα + ( ρA − ρv ) g / α º ( ρA + ρv ) 2 ¼ (5.213) uA − uv > ¬ where α = 2π / λ is the wave number. Surface tension and gravity tend to stabilize the interface (for this configuration only). The right side of this inequality has a minimum when the wave number is equal to a critical wave number α crit : ρA ρv ª ( ρ − ρv ) g º α crit = « A » σ ¬ ¼ It follows from eq. (5.213) that 1/ 2 (5.214) 384 Transport Phenomena in Multiphase Systems ª 2 ( ρA − ρv ) º ª σ ( ρA − ρv ) g º (5.215) uA − uv crit = « »« » ρA ρv2 ¬ ¼¬ ¼ For motionless liquid over motionless vapor ( uA = uv = 0 ), we obtain from eq. (5.215) 1/ 2 1/ 4 ª ( ρA − ρv ) g º (5.216) » σ ¬ ¼ The critical wavelength corresponding to the critical wave number is 1/ 2 α > α crit = « ª º σ λc = 2π « (5.217) » « ( ρA − ρv ) g » ¬ ¼ A perturbation with a wavelength greater than λc will grow and result in instability. If the length of the interface in the x-direction is less than λc , the interface is stable because a perturbation of wavelength greater than λc cannot arise. A specific value of α exists where β in eqs. (5.210) – (5.212) is at its maximum. Its value is ª ( ρ − ρv ) g º (5.218) α max = « A » 3σ ¬ ¼ The disturbance wavelength corresponding to α max is referred to as the most dangerous wavelength, λD – 1/ 2 1/ 2 ª º 3σ λD = 2π « (5.219) » = 3λc « ( ρA − ρv ) g » ¬ ¼ which has many applications, including derivation of the critical heat flux for pool boiling. When the direction of the gravity is parallel to the direction of the liquidvapor co-current flow as shown in Fig. 5.24, the gravity will not have a significant effect on the pressures in the liquid and vapor phases. Equation (5.213) becomes σα ( ρA + ρv ) 2 uA − uv > (5.220) 1/ 2 ρA ρv which is the condition for an unstable vertical interface, and this condition is termed the Kelvin-Helmholtz instability. Example 5.7 In a stratified horizontal two-phase flow system, saturated water vapor flows above saturated liquid water. The pressure of the two-phase system is at one atmosphere. Determine the critical velocity, the Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 385 corresponding critical wavelength wavelength λD . λc , and the most dangerous Solution: The properties of saturated vapor and liquid water at one ρv = 0.5974kg / m3 , ρA = 958.77kg / m3 , and atmosphere are σ = 0.05891N / m. The critical velocity is obtained by eq. (5.215), i.e., uA − uv crit ª 2 ( ρA − ρv ) º =« » ρA ¬ ¼ 1/ 2 ª σ ( ρA − ρv ) g º « » ρv2 ¬ ¼ 1/ 4 ª 2 × ( 958.77 − 0.5974 ) º ª 0.05891 × ( 958.77 − 0.5974 ) × 9.8 º =« »« » 958.77 0.59742 ¬ ¼¬ ¼ = 8.87m/s The critical wavelength is obtained by eq. (5.217), i.e., 1/ 2 1/ 4 ª º σ λc = 2π « » « ( ρA − ρv ) g » ¬ ¼ 1/ 2 ª º 0.05891 = 2π « » = 0.0157m = 15.7mm « ( 958.77 − 0.5974 ) × 9.8 » ¬ ¼ The most dangerous wavelength can be obtained by eq. (5.219), i.e., λD = 3λc = 3 × 0.0157 = 0.0273m = 27.3mm 1/ 2 5.6.2 Surface Waves on Liquid Film Flow As will become evident in later chapters, condensation or evaporation at the interface of a thin liquid film flowing over a solid surface is often encountered in engineering applications. This type of process is commonly found in chemical and mechanical engineering equipment such as condensers, long-tube evaporators, wetted wall columns, and cooling towers. The fact that current literature shows continued interest in this general area affirms its continued importance. The problem of the transport from a liquid layer running down a vertical wall has been considered in literature with particular attention to heating, evaporation, condensation, and gas absorption for various applications. The flow in such a liquid layer may be either (a) laminar with a constant thickness, or (b) laminar with variable thickness due to waves, or turbulent with waves moving down the interface. Kapitza (1964) predicted theoretically that the mean film thickness is reduced by approximately 8% due to ripples on the interface. Many experiments have focused on determining the values of the wave properties in both laminar and turbulent flows. Experimental data on the mean film thickness of wavy laminar 386 Transport Phenomena in Multiphase Systems flow falls between the predictions of Kapitza and Nusselt (smooth interface) theories. Although the wave effect does not have any significant effect on the mean film thickness in wavy laminar regions, it plays a major role in heat and mass transfer processes. The structure of the liquid film running down a plane usually shows a randomly distributed wave on the surface, except at low flow rate at the entry region, where a periodic motion exists. Figure 5.26, which is taken originally from Rogovan et al. (1969) shows such a variation as obtained at a Reynolds number of 54 for water. Comparisons of the experimental and theoretical values of the transport in falling film should be supported by a specification of the nature of the waves, for it is through the interaction of the normal velocities produced by the waves and the distribution of the transported quantity that the additional transport arises. But this is complicated by the existence of an initial wave-free length, depending upon the way in which the film was formed by the liquid supply, and by a region of wave development that ostensibly culminates in a steady regime at a distance far enough down the height of the film. Brauer (1956), for example, made measurements on water, and mixtures of water and diethylene glycol, at a distance of 1.3 m from the point of film initiation, and found sinusoidal waves to begin at a Reynolds number, Re = 4Γ / μ = 1.2 / Ka1/10 , where Ka = ν 4 ρ 3 g / σ 3 is the Kapitza number. The ratio of the crest height of the waves to the mean film thickness increased with the Reynolds number to about twice for this Reynolds number than for higher Reynolds number. This ratio remained the same up to a Reynolds number of 140 / Ka1/10 . Between these Reynolds numbers the waves become distorted and no longer sinusoidal; above the higher Reynolds number there were secondary capillary waves on the surface. Throughout, up to 4Γ / μ = 1600, the average film thickness remains essentially that given by the Nusselt 0.4 0.2 0 a Thickness of film, mm 0.4 0.2 0 0.4 0.2 0 0.4 0.2 0 0 0.1 0.2 b c d 0.3 0.4 Time, sec. 0.5 Figure 5.26 Change in the local liquid film thickness at Re=54 at a distance from a sprayer of (a) 130, (b) 310, (c) 445, and (d) 670 mm (Faghri and Payvar, 1979). Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 387 analysis. But in the range where the ratio of the crest height to the mean film thickness is almost constant, the wave frequency increases with the Reynolds number to indicate a basis for a proportional increase in heat and mass transfer. Koizumi et al. (2005) studied the behavior of a liquid film flowing down the inner surface of a pipe with countercurrent gas flow. In the experiment, the gas used was air, and silicone oils of 500, 1000, and 3000 cSt as well as water, were used for the liquid phase. Three primary parameters were correlated: the film thickness, wave velocity, and wavelength. The film thickness was given three different values: minimum, maximum and mean film thickness. The film height was determined by taking photos of the film and measuring its thickness. The wave velocity was determined by seeing how far a wave peak traveled during the time interval between picture frames. The wavelength was determined in a similar manner; that is, photographs were taken, and the distance between peaks of the wave were measured. Figure 5.27 presents experimental results for water and silicone oil with no airflow for dimensionless δ * versus the Reynolds number, where δ * = ( y /ν ) τ w / ρA , is the film thickness, w is the wall shear stress and ρA is the liquid density. When there is no airflow, the maximum film thickness of the 500 cSt and the 1000 cSt silicone and water films is much greater. To correlate the wave properties using a non-dimensional analysis, the nondimensional wave velocity, N w = ww /ν 1/ 3 , and maximum film thickness, g 102 Max. thickness All average Min. thickness Open: Water Solid: Silicone500cSt Double: Silicone100cSt Half: Silicone3000cSt 101 (δ /ν ) τ w / ρA 100 Nusselt Universal Velocity Profile 10-1 -2 10 103 Liquid film Reynolds number, 4Γ / μ 10-1 100 101 102 104 Figure 5.27 Dimensionless liquid film thickness as a function of Reynolds number for falling liquid film (adopted from Koizumi et al., 2005). 388 Transport Phenomena in Multiphase Systems + δ max = δ max ( g /ν 2 )1/ 3 , are found according to the Buckingham pi theorem. They are functions of the Reynolds number ( Re = 4Γ / μ ), the Morton number ( K F = ρ 3ν 4 g / σ 3 ) and the nondimensional wavelength ( N λ = λ ( g /ν 2 )1/ 3 ), where ww is the wave velocity and is the wavelength. Nosoko et al. (1996) developed the following correlations for wave velocity and maximum film thickness: 0 0.31 N w = 0.68K F.02 N λ Re0.37 (5.221) (5.222) Equations (5.221) and (5.222) require information about the wavelength; therefore, the correlations are not closed. The wave velocity and maximum film thickness should be evaluated only from the film flow rate and physical dimensions, if possible. The wavelength correlation was incorporated into the Nosoko correlations by Koizumi et al. (2005), and the constants and exponents were modified for better accuracy. The final forms are N λ = 14.9Fr1.29 Re −0.133 (5.223) 0 N w = 0.88K F.008 Fr 0.977 Re0.214 + 0 0.39 δ max = 0.26 K F.044 N λ Re0.46 (5.224) (5.225) δ + max 0 = 1.09 K F.021Fr 0.316 Re 0.424 where Fr = ww / gδ and δ is the mean film thickness. The thickness of these liquid films is small enough that boundary layer assumptions are valid. However, these films are not thin enough for disjoining pressure to play any significant role. Figure 5.28 shows a slow liquid flow on an inclined surface with negligible shear stress at the free surface. The gravitational and liquid viscous forces, as well as surface tension effects, are dominant. As the mass flow rate of the liquid increases, waves on the liquid film interface can be observed, as indicated by Fig. 5.28. The flow regimes of the liquid film include smooth laminar, wavy laminar, and turbulent, which correspond to different film Reynolds numbers. The Reynolds number is defined as 4Γ 4wδ (5.226) Reδ = = where Γ is the mass flow rate per unit width and w is the mean velocity of the liquid film. The liquid film is laminar if the film Reynolds number is below 30. A wavy flow regime can be observed when the film Reynolds number is between 30 and 1600. The liquid flow becomes turbulent when the film Reynolds number is above 1600. For a fully-developed steady-state downward laminar flow, the momentum equation reads ( ρA − ρv ) g sin θ d 2w + vA 2 = 0 (5.227) ρA dy μA μA Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 389 Vapor δ w Figure 5.28 Liquid film on an inclined surface. Integrating eq. (5.227) twice and considering the boundary conditions w = 0 at y = 0 and dw / dy = 0 at y = δ , the velocity distribution in the liquid film is obtained. ( ρ − ρv )(2δ y − y 2 ) g sin θ w= A (5.228) 2 μA The mean velocity in the liquid film is then ( ρ − ρv ) g sin θδ 2 1δ (5.229) w= w ( y ) dy = A 3μA δ0 which is obtained under assumption that the liquid flow is steady-state. Under these assumptions (Nusselt analysis) and assuming ρA  ρv , the dimensionless film thickness, and film velocity at the surface δ , wδ and the mean velocity w ³ are given for θ = 90D below: δ+ =δ ¨ §g· 2¸ ©ν ¹ 1/ 3 = 0.908Re1/ 3 (5.230) g2 δ (5.231) 2ν 2 w = wδ (5.232) 3 When waves are present, the liquid film flow will be unsteady. However, it is assumed that eq. (5.229) is also valid for unsteady film flow. Consideration of hydrostatics and the Young-Laplace equation at the interface yields: ∂pA ∂pv ∂δ ∂ 3δ (5.233) = + ( ρA − ρv ) g cos θ −σ 3 ∂z ∂z ∂z ∂z Since there is no motion in the vapor phase, the pressure gradient in the vapor phase satisfies wδ = 390 Transport Phenomena in Multiphase Systems ∂pv = ρv g sin θ (5.234) ∂z The momentum equation for an unsteady-state flow in the liquid film is § ∂w ∂p ∂w ∂w · ∂2w +v + w ¸ = − A + ρA g sin θ + μA 2 ρA ¨ (5.235) ∂y ∂z ¹ ∂z ∂y © ∂t Substituting eqs. (5.233) – (5.234) into eq. (5.235), the momentum equation becomes ∂w ∂w ∂w ( ρA − ρv ) g § ∂δ · σ ∂ 3δ ∂2w +v +w = sin θ − cosθ + + vA 2 ¨ ¸ ∂t ∂y ∂z ∂z ¹ ρA ∂z 3 ρA ∂y © (5.236) Integrating eq. (5.236) and considering the continuity equation, two governing equations for δ ( z , t ) can be obtained: 3 3w2 ∂δ 1 ∂δ 3w § · σ ∂δ = ( ρA − ρv ) g ¨ cos θ − sin θ ¸ − = vA 2 (5.237) 3 δ ∂z ρA ∂z δ © ¹ ρA ∂z ∂δ ∂δ = 3w (5.238) ∂t ∂z These equations can be solved along with the energy equation. Our purpose, however, is to find the conditions of stability of the film. Assuming that δ = δ 0 + δ ′ and w = w0 , and noticing from eq. (5.229) that 3w v 1 (5.239) ( ρA − ρv ) g sin θ = 02 A ρA δ0 the following is obtained: 2 3w0 ∂δ ′ 1 ∂δ ′ σ ∂ 3δ ′ = ( ρA − ρv ) g cos θ − (5.240) δ 0 ∂z ρA ∂z ρA ∂z 3 The fluctuating component δ ′ is assumed to be a sinusoidal wave of the form δ ′ = δ *eiα ( z −ct ) (5.241) where α is the wave number and c is the wave velocity. For α > 0, which is the condition under which waves exist, it follows from eqs. (5.240) and (5.241) that ª º 3μA2 cos θ δ >« 2» « ( ρA − ρv ) ρA g sin θ » ¬ ¼ 1/ 3 (5.242) which can also be written in terms of film Reynolds number, Reδ = 4 ρA w0δ / μA , i.e., Reδ > 4cot θ (5.243) Equation (5.243) indicates that the critical Reynolds number at which waves appear equals cot θ . For a vertical surface where θ = 90D , eq. (5.243) become Reδ > 0 , which means that waves can be present no matter how small the film Reynolds number is. This result contradicts experimental observation, because waves are present only when the film Reynolds number reaches a certain value. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 391 This contradiction can be solved by considering the amplification rate of wave amplitude. If amplification over the time interval required to travel 100 times the liquid film thickness is considered, the amplification factor is (Benjamin, 1957; Carey, 1992) −1 ­ª ½ º δ ° 4 / 3 1/ 3 § σ · 8° (5.244) ν = exp ® «0.31 A g ¨ ¸ » Reδ / 3 ¾ δ t =0 © ρA ¹ » °« ° ¼ ¯¬ ¿ For saturated water at one atmosphere, ν A = 2.91 × 10−7 m 2 /s, σ = 0.0589N/m, and ρA = 958.77kg/m3 , and eq. (5.244) becomes δ 8 = exp ( 2.08 × 10−5 Reδ / 3 ) δ t =0 (5.245) The dependence of disturbance amplification versus film Reynolds number is shown in the following table: Reδ 10 1.010 20 1.063 40 1.476 80 11.85 120 1466.2 160 6.588×106 δδ t =0 It can be seen from the above table that the wave becomes significant after Reδ is greater than 40. Waves increase the heat transfer coefficient during evaporation (also condensation) of wavy films as compared to a smooth surface due to an increase in interfacial surface area and mixing action (Faghri and Seban, 1985). Heat transfer during condensation and vaporization over a wavy liquid film will be discussed in detail in Chapter 8 and 9, respectively. 5.7 Numerical Simulation of Interfaces and Free Surfaces Many engineering applications involve interfacial phenomena, because of the high heat transfer rates that can be achieved. There are many complexities that need to be addressed when modeling an interface, since an interface is generally irregular, involves mass transfer, is three-dimensional, and is not at a fixed location. There are two different approaches when considering an interface: a continuum and a noncontinuum approach. In many problems, the scale of the entire computational domain is much larger than the scale of the interfacial thickness, therefore the interface can be considered a planar surface and the continuum approach will give highly accurate results. With the development of nanotechnology, the scales of systems are becoming very small. Therefore, neglecting the thickness of the interface is no longer a valid assumption, and a noncontinuum approach is needed. This approach can be achieved through molecular dynamic simulations. A discussion tracking an interface using the continuum approach is given in Section 5.7.1, while the noncontinuum approach is given in Section 5.7.2. 392 Transport Phenomena in Multiphase Systems 5.7.1 Continuum Approach: Interface Tracking Techniques There are many engineering problems that involve separated phases, in which the interface between phases is clearly defined but not at a fixed location. The interface can be assumed to be a planar surface in these situations, and therefore a continuum approach is valid. In these problems, it is necessary to solve for both the phases as well as the interfacial location. These problems can be solved numerically on an Eulerian or Lagrangian mesh. An Eulerian mesh is stationary and defined prior to the start of a solution. When using an Eulerian mesh, the interface is tracked by solving an additional scalar equation. In the Lagrangian approach, a boundary of the mesh is aligned with the interface, and this boundary moves with the interface. When thinking of a multiphase system from a continuum approach, in the bulk region a phase is continuous and is discontinuous at an interface between different phases. In general, the interface is free to deform based on the nature of the flow; therefore, it is difficult to efficiently capture an interface between phases with just one model. Consequently, there have been strong efforts to track an interface based on several different techniques, each with its own pros and cons. An actual interface, represented by the different numerical techniques, is presented in Fig. 5.29. The techniques are as follows: 1. Lagrangian Techniques: grid and fluid move together, interface is directly captured 2. Stationary Grid Approach: standard CFD modeling in cells that contain only a single phase, special consideration taken in cells in the vicinity of an interface 3. Phase Interface Fitted Grid: interface directly tracked as boundary, rest of grid moves as a function of interfacial movement, “semi-Lagrangian” 4. Front Tracking Methods: stationary grid is used and modified near the interface, so that the grid is aligned with the interface; combination of 2 and 3. Each of the above methods will be briefly decribed below. Lagrangian Techniques In a Lagrangian technique, fluid particles are tracked directly. Therefore, a particle on the interface will also be tracked directly, so the interface will automatically be resolved. The Lagrangian approach is used by Shopov et al. (1990) to model a droplet interacting with a wall. Even though a Lagrangian approach directly captures an interface, an Eulerian approach is often preferable to a Lagrangian approach to solve the bulk flow of a fluid. An Eulerian approach is generally preferable, because fluid flow is usually complex, and fluid particle path lines often intersect and/or disperse, making fluid interaction very difficult Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 393 Actual Interface Interface Modeled as a Boundary (a) Actual interface Interface approximated by piecewise interpolation from cell volume fraction (c) Phase interface fitted grid Interface captured by a separate interfacial grid (b) Volume of Fluid model (d) Front tracking/capturing Figure 5.29 Interfacial representations for different interface tracking techniques. to manage with a Lagrangian method. Techniques 2, 3, and 4 all use an Eulerian point of view to solve the governing equations in the bulk fluids with special consideration to track the interface. Stationary Grid Approach The second class of numerical techniques for multiphase fluids is based on a stationary grid where the fluid interface is captured directly. The first of such approaches is the marker-and-cell approach originated by Harlow and Welch 394 Transport Phenomena in Multiphase Systems (1965). In this approach, massless particles are introduced into the flow field, and the locations are projected from their interpolated velocities. Cells with a particle are considered to have one phase in them, and cells without a particle do not have that phase in them. An interface is considered where cells with particles are neighbored with cells without particles. For efficiency, this method was extended to only track particles on the surface, (Nichols and Hirt, 1975). Further development of this class of models lead to the Volume of Fluid method (VOF) by Hirt and Nichols (1981); they used a donor-acceptor method to effectively eliminate numerical diffusion at an interface. The VOF method is one of the most widely used routines to solve a free surface problem on an Eulerian mesh. In this method, there is one velocity, pressure, and temperature field, and it is shared by all of the phases modeled. The interface between the phases is tracked by the volume fraction of phase k, ε k . The volume fraction equation is the continuity equation of phase k divided by the density of that phase. Π º 1 ª∂ m′′′ » (5.246) « (ε k ρk ) + ∇ ⋅ (ε k ρk V ) = jk ρ k « ∂t » j =1( j ≠ k ) ¬ ¼ ¦ When the volume fraction is between 0 and 1 in a computational cell, that cell is considered to be an interfacial cell. For a two-phase system, the phases are k and j. When the volume fraction is 1, that cell is occupied by only phase k, and when it is 0, that cell is occupied by only phase j. The sum of the volume fractions is unity. ¦ε k =1 Π k =1 (5.247) Therefore, (k-1) volume fraction equations need to be solved. The fluid properties, such as density, viscosity, and thermal conductivity, are calculated by their volume weighted average. Φ eff = ¦ε Φ k k =1 Π k (5.248) The overall continuity equation in conjunction with the VOF method is: ∂ (5.249) ( ρeff ) + ∇ ⋅ ( ρeff V ) = 0 ∂t The continuity equation is the same as eq. (3.38), except it uses the effective density and the velocity field is shared by both phases; therefore, subscript k is dropped. The momentum equation is Π Π ∂ ′ m′′′V jk (5.250) ( ρeff V ) + ∇ ⋅ ( ρeff VV ) = ∇ ⋅τ eff + ρeff X + ∂t k =1 j =1( j ≠ k ) ¦¦ The momentum equation is also the same as eq. (3.52), with the same exceptions as the continuity equation and the momentum transfer due to phase change. In the energy equation, the enthalpy is mass averaged instead of volume averaged. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 395 n (a) (b) (c) Figure 5.30 (a) An actual interface between two phases, (b) an interfacial representation using the Donor-Acceptor scheme, and (c) a piecewise linear reconstruction scheme with the VOF method heff = 1 ρeff ¦ε k =1 Π k ρ k hk (5.251) Therefore the energy equation, neglecting pressure effects and viscous dissipation, is: Π Π ∂ ′′′ ρeff heff ) + ∇ ⋅ ( ρeff Vheff ) = −∇ ⋅ q′′ + m′′′ hk + qeff (5.252) ( eff jk ∂t k =1 j =1( j ≠ k ) ¦¦ The energy equation has the same form as eq. (3.79), except effective properties are incorporated and a latent heat term due to phase change is added. The continuity, momentum, and energy equations can be solved by standard solution procedures, such as the SIMPLE class of algorithms in Patankar (1980). However, if the volume fraction equation is solved using a standard implicit or explicit scheme, the interface will quickly lose resolution due to numerical diffusion. The numerical diffusion will lead to inaccurate solutions and/or a solution that will not converge. Therefore, special consideration must be taken on interfacial cells to construct the interface so that the advection transport of fluid is representative of the physical problem. An actual interface between two phases, as well as the corresponding interface represented by two special interpolation schemes, is presented in Fig. 5.30. 396 Transport Phenomena in Multiphase Systems As noted before, the first of such methods is the donor-acceptor scheme proposed by Hirt and Nichols (1981). If a cell is an interfacial cell, 0 < ε k < 1 , the fluid will be rearranged in the cell to be on one face, as show in Fig. 5.30(b). The face on which the fluid will be rearranged will depend on the normal direction of the interface. The normal direction of the interface with respect to phase k can be calculated by the gradient of the volume fraction: ∇ε k nk = − (5.253) ∇ε k The component in which the normal is the greatest will occur where the fluid is perpendicular to the interface while the interface is either horizontal or vertical. Once this is done, one cell is designated as a donor, and its neighbor is designated as an acceptor. The amount of fluid leaving the donor cell is exactly equal to the amount of fluid entering the acceptor through each computational face. Also, the maximum amount of fluid that can leave a cell is equal to the amount of fluid in that cell or the amount that would make another cell fill with fluid. A more refined interface interpolation scheme was developed by Youngs (1982), in which the interface was approximated as piecewise linear in each cell. The normal of the reconstructed interface in each cell is the same as the normal calculated in eq. (5.253), as shown in Fig. 5.30(c). The advection of fluid through each face is calculated in a manner similar to the donor-acceptor scheme. Both the donor-acceptor and piecewise linear interpolation can only be run in transient simulations, and the time step must be kept small enough so that the fluid near the interface will only advance by one cell at a time. Also, the resolution of the interface is limited to the grid spacing of the computational mesh. Any surface waves that are smaller than the spacing of the mesh will be smoothed out, as seen in Fig. 5.30 on the left side of the interface. Surface tension effects are important in many free surface problems. The surface tension effects can be applied as a body force in the momentum equations. This method is called the continuum surface force (CSF) model and was proposed by Brackbill et al. (1992). The curvature is defined as the divergence of the surface normal vector. Kk = ∇ ⋅ nk (5.254) The volume force due to surface tension, Fσ , is Π k −1 ε j ρ j K k ∇ε j + ε k ρ k K j ∇ε k −σ jk Fσ = (5.255) 1 k =1 j =1 ( ρ j + ρk ) 2 If only two phases are present, the body force reduces to ρeff K k ∇ε j (5.256) Fσ = σ jk 1 ρk + ρ j ) ( 2 ¦¦ Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 397 It can be seen that the body force is proportional to the cell density. It is important to note that when surface tension forces are large compared to other flow characteristics, numerical inaccuracies of the surface tension forces in the CSF model create artificial currents called parasitic currents. Much work has been done to eliminate parasitic currents, such as the second gradient method proposed by Jamet et al. (2002). Despite the adverse effects of parasitic currents, the VOF method has been widely used and can give reasonably accurate results for a wide range of applications. It is robust in its handling of free surface problems with large interface distortion and can easily handle problems in which the free surface breaks apart, such as droplet formation. One other drawback of the VOF method is that the interface resolution is limited to the grid spacing. Therefore, a refined mesh is needed anywhere the interface is going to travel. This refinement can lead to many computational cells in regions of the mesh resided in by the interface for a short period of time, which will increase the total computational time of the solution. Therefore, advanced remeshing algorithms are needed for problems of this type. Phase Interface Fitted Grid The third method, in which a phase interface fitted grid is employed, can be very useful for multiphase problems in which the interface does not deform greatly or break apart (although it can be handled with complex algorithms). The benefit of these models is that they directly capture an interface; therefore, a phenomenon occurring at a phase interface is not subject to interpolation error, such as the interface location itself. The third class of numerical method used to capture a free surface is the semi-Lagrangian approach. Since a Lagrangian approach is not feasible for most fluid mechanics problems, only the interface is tracked in a Lagrangian manner, while the rest of the mesh can be manipulated to maintain good mesh qualities. Therefore, this type of technique can be deemed “semi-Lagrangian,” and is a front tracking technique. Rice and Faghri (2005a) developed a semi-Lagrangian approach that works effectively on the interface between phases with phase change. They used a transient method in which the present interfacial velocity, VI , is calculated based on the error in mass flux through each computational face from the previous time step. Therefore the interfacial velocity is lagging by one time-step. The conservation of mass between phase k and phase j is ′′ m′′ = ρ k ( Vk − VI ) ⋅ n = ρ j ( V j − VI ) ⋅ n = mactual (5.257) The interfacial continuity equation represents the actual mass flux at the interface. Instead of implicitly calculating the interfacial velocity of the next time step (n+1), the interfacial velocity from the previous time step (n) is used. By lagging the interfacial velocity, there is an inherent error in mass flux at the interface. 398 Transport Phenomena in Multiphase Systems ρ VIn ⋅ nΔtA merror Δt min +1Δt ρ VIn +1 ⋅ nΔtA = mI Figure 5.31 Error in interfacial location due to the lag in time. m′′n +1 = ρ k Vkn +1 − VIn ⋅ n − m′′n +1 = ρ k VIn +1 − VIn ⋅ n error actual ( ) ( ) (5.258) The normal velocity of phase j at the interface is fixed by continuity. However, the normal velocity of phase k at the interface is not fixed, and fluid is allowed to pass through this interface without regard to the interfacial mass balance. Figure 5.31 demonstrates how the interfaces moves and shows that there is an error in the mass of the system due to the interfacial movement. In order to compensate for the error in mass through each computational face, the new velocity can be found by rearranging eq. (5.258) and integrating over time and area. Δt ³ AI ρ k VIn +1 ⋅ n dA = Δt ³ AI ( m′′ n +1 error + ρ k VIn ⋅ n dA ) (5.259) In discretized form, eq. (5.259) becomes ª Δt ρ k VIn +1 ⋅ n Aº = ª Δt m′′n +1 + ρ k VIn ⋅ n Aº = mi (5.260) error ¬ ¼i ¬ ¼i The subscript i represents a single face on the interface. Rice and Faghri (2005a) suggest moving the interface at the beginning of each time step based on the error in mass flux and interfacial velocity of the previous time step, so that the mass the interface gains by movement is equal to mi. They specified a direction in which the interface location (δ I ) will always be a function of δ I = f (η ) , and moved the nodes on the interface normal to this direction. Specifying the coordinate system can prevent the interfacial face area of a single cell from becoming zero or negative. Some knowledge of the problem is necessary in order to specify the direction in which the interface can move. The pressure of phase k at the interface should be defined from the interfacial momentum equation. When a second phase also influences the interface, this interface will move with the interface of phase k. However, the mass flux that ′′ passes through this interface should be identically equal to mactual . This side of ( ) ( ) Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 399 the interface affects the interfacial movement through the momentum equations. When surface tension effects are important, they are directly applied as an interfacial boundary condition in the interfacial momentum equation. This technique was shown by Rice and Faghri (2005b) to eliminate the parasitic currents caused by surface tension effects that are experienced in the VOF method. Also, this technique directly tracks the interface; therefore, the interface is always directly resolved, no matter what the mesh spacing is. For problems with large fluid distortion, it is very difficult to define a coordinate system of which the interface will always be a function; therefore, this coordinate system will need to be able to move over time. If the interface is largely distorted, advanced remeshing techniques will be needed for the computational domain away from the interface to maintain a good quality of computational cells. Rice and Faghri’s (2005a, b) technique used a transient solver that allowed a small amount of liquid mass to pass through the interface. During each new time-step, the liquid that passed through the interface during the previous timestep was recaptured in a manner in which mass is directly conserved. Even though this technique is generally not suitable for interfaces undergoing large deformations, such as jet breakup, it has the advantage of directly capturing the interface. Therefore, this technique offers a great advantage in problems involving multiple components, whose interfacial values are not continuous. Front Tracking Methods The fourth class is a combination of the second and third class of multiphase and interface modeling techniques, and is called the front tracking technique. One form was developed by Glimm et al. (2001). This method uses a fixed grid that is modified only near the interface. The interface is modeled as a separate moving grid. The interface can move, and it separates the stationary grid cells into two cells near the interface, representing each fluid. Therefore, each fluid is treated separately. Another front tracking technique developed by Tryggvason et al. (2001) also uses a separate grid to track the interface; however, the phases are considered together, and a single set of governing equations is solved for the whole field. Therefore, the stationary cells that lie below the interface have a shared density and velocity field. Also, since the interfacial area is free to expand or contract, the subgrid that tracks the interface must have the ability to add or subtract nodes as needed. Further consideration of this method is discussed in the following section. The main concept of the front tracking method is to identify a fluid with a Heaviside function, H. This function is one at a location where a particular fluid exists and zero elsewhere. The interface is located where the H function changes from zero to one, which is the location where the delta function, , is nonzero. H can be expressed as a function of the function. (5.261) H ( x, t ) = δ ( x − x I ) dV ³ ΔV 400 Transport Phenomena in Multiphase Systems The integration is performed over an entire volume, V, with δ being nonzero only on the interface between the phases, AI, at location xI. The calculation of the fluid properties, Φ , is straightforward with the Heaviside function. Φ ( x, t ) = Φ k H ( x, t ) + Φ j (1 − H ( x, t ) ) (5.262) Note, that when these values are integrated over a finite volume, the properties used in the calculation are volume averaged. The calculation of the interface is done in a Lagrangian fashion. dx I ⋅ n = VI ⋅ n (5.263) dt To write the equivalent expression using a kinematic equation, the gradient of H must be defined ∇H = ³ AI n δ ( x − x I ) ds (5.264) where s is a surface that is integrated over AI. Therefore, the kinematic equation using H to represent the interface is: ∂H = − VI ⋅ ∇H = − VI ⋅ nδ ( x − x I ) ds (5.265) AI ∂t To calculate the interfacial velocity, the overall continuity equation can be considered. ∂ (5.266) ( ρ ) = −∇ ⋅ ( ρ V ) ∂t The evaluation of fluid properties and flow values uses the H function as shown in eq. (5.262). The overall continuity equation can be written as: ∂ (5.267) ( ρk H + ρ j (1 − H ) ) + ∇ ⋅ ( ρk Vk H + ρ j V j (1 − H ) ) = 0 ∂t Using the chain rule, the continuity equation is broken into three regions: the bulk regions for phase k and phase j, and the interface. The interfacial continuity equation is: ∂ (5.268) ( ρk − ρ j ) ∂t ( H ) = − ( ρk Vk + ρ j V j ) ⋅ ∇ ( H ) Applying eqs. (5.264) and (5.265) to eq. (5.268) yields: ³ ³ (ρ AI k − ρ j ) VI ⋅ nδ ( x − x I ) ds = ³ (ρ V AI k k + ρ j V j ) ⋅ nδ ( x − x I ) ds (5.269) The interface velocity can be calculated using eq. (5.263). The momentum equation is: ∂ρ V + ∇ ⋅ ( ρ VV ) = ∇ ⋅ τ + ρ g + σκ nδ ( x − x s ) ds (5.270) AI ∂t The last term represents the surface tension effects where the twice the mean curvature is κ . Also, the energy equation can be written as: ∂ ( ρ c pT ) + ∇ ⋅ ( ρ Vc pT ) = −∇ ⋅ q′′ − AI m′′ ªhkj + ( c p,k − c p, j )δ ( x − xs )ºds ¬ ¼ ∂t (5.271) ³ ³ Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 401 Note that pressure effects and viscous dissipation are neglected. The last term represents the latent heat release in a phase change process. The interface is made of a separate mesh that goes through the stationary mesh. It generally does not lie on top of the points in which the conservation laws are applied. Therefore, the interfacial values are calculated by grid interpolation, and mass and momentum of the advection front is not necessarily conserved. The error in mass and momentum conservation at the interface can be reduced to an acceptable level by refining the mesh. 5.7.2 Noncontinuum Approach: Molecular Dynamic Simulation More detail about what is actually happening at an interface requires a different approach, because the continuum assumption breaks down. For these situations, a molecular dynamic simulation is needed. Research at the molecular level on vapor formation, condensation, and other interface phenomena is fairly scarce due to the limitations of experimental instruments, because the usual macroscopic continuum numerical models do not describe the phenomena at such a small scale. A promising approach to analyzing phenomena at the molecular level is molecular dynamic simulation (MDS). By solving Newton’s equation of motion for every molecule in the system numerically, information on the macroscopic system can be found. The basic idea behind MDS is that as every substance is made of particles (atoms or molecules), if the dynamics of these particles, such as velocity, position, and interaction force, can be determined, then physical macroscopic properties such as temperature, pressure, volume, etc. may also be found. The scale treated by MDS is only about 10 nm in space and 10 ns in time with the latest computers. The system needs to be decomposed into elements that can be examined using MDS. Since MDS gives only the transient position and velocity of the molecules in the system, macroscopic quantities such as temperature, internal energy, pressure, heat flux, and momentum flux are obtained by averaging the MDS system to interpret the phenomena. In MDS, each atom is treated as a point mass and the interactions between atoms are described by Newton’s second law of motion. The equation of motion for each atom or molecule under consideration may be expressed as N dφij (r ) rij d 2r Fi = mi ai = mi 2i = − (5.272) dr rij dt j ≠i ¦ where Fi is the total force vector acting on the ith molecule; N is the number of molecules; φij is the pair potential between the ith molecule and the jth molecule; rij represents the position vector and rij is the distance between the ith molecule and the jth molecule; mj is the mass of the ith molecule; and ai is the acceleration vector of the ith molecule. The pair potential approximation does not give the correct potential energy between two molecules when they are isolated in space. 402 Transport Phenomena in Multiphase Systems Therefore, more accurate models are required to give the correct overall potential of the system. The most basic model is the Lennard-Jones potential model for spherical molecules, which is presented in Chapter 1, eq. (1.28). MDS has been used to determine the evaporation and condensation coefficients of fluids, such as water and argon, because experimental determination of these values is very difficult. This is due to the large temperature increase in the interface region, which is made up of only a few molecule layers and cannot be directly measured. Some examples of MDS in recent research are: Yasuoka and Matsumoto (1995) used MDS to obtain the condensation coefficient of argon, methanol, and water at certain temperatures; Matsumoto (1998) used MDS to explore evaporation-condensation dynamics of pure fluids under both equilibrium and non-equilibrium conditions; Wu and Pan (2003) investigated the evaporation of a thin argon layer into a vacuum using MDS with the Lennard-Jones potential; Yang and Pan (2005) examined the evaporation of a thin layer of water and determined an evaporation coefficient. Yang and Pan (2005) used the TIP4P model as the effective pair potential between water molecules. The model consisted of four sites: two positively charged hydrogen atoms, one negatively charged hydrogen atom, and one Lennard-Jones site for oxygen. A total of 1000 water molecules were used for the simulation with dimensions of 25 × 25 Å. A periodic boundary condition was used for all six boundaries. The simulation shows that the hydrogen bond in water has a significant effect on the molecules at the interface and may reduce the evaporation coefficient. This study showed that a combination of MDS and the classic Schrage model might provide a methodology to find the evaporation coefficient at different temperatures and pressures. Figure 5.32(a) shows the density distribution for several different temperatures in the normal direction to the interface at equilibrium. The figure shows three zones: liquid, interface, and vapor. The density profile may be fitted to a hyperbolic-tangent function, as shown in Fig. 5.32 (b): 1 1 ª z − z0 º (5.273) ρ ( z ) = ( ρA + ρv ) − ( ρA − ρv ) × tanh « » 2 2 ¬d¼ where ρA and v are the liquid and vapor densities, respectively; z0 is the position of the Gibbs dividing surface d is a measure of the interfacial thickness and is related to the interfacial thickness, δ I , by δ I = 2.1792d (5.274) The interfacial thickness ranges from 4 to 8 Å and increases with temperature by δ I = −10.95 + 0.04068T (5.275) where δ I and T are in Å and K, respectively. The interfacial thickness for water is much smaller than that for argon, which is ranged from 16 to 20Å. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 403 Figure 5.32 (a) Density profile at various different temperatures, (b) density profile at 423.15K and its fitted curve of hyperbolic tangent form (Yang and Pan, 2005; Reproduced by permission of Routledge/Taylor & Francis Group, LLC). Molecular dynamics simulation can apply not only to liquid-vapor interfaces, but also to solid-liquid interfaces, to provide detailed information at the microscopic level. For example, conventional disjoining pressure theory, presented in Section 2.6.4, provides a base for determining the effects of wallfluid attractive forces but ignores many details of the variation in thin liquid film. Molecular dynamics simulation was used by Carey and Wemhoff (2005) to examine the detailed features of the structure and thermophysics of thin liquid films on solid surfaces. The film features and equilibrium vapor pressure were determined for films of varying thickness and were compared to the conventional disjoining pressure results. The molecular dynamics simulation used combines traditional NVE-type MD simulation with a stochastic boundary condition to represent the equilibrium of pressure. This method allows monatomic fluids that interact via the Lennard- 404 Transport Phenomena in Multiphase Systems Jones potential ª§ σ LJ ·12 § σ LJ ·6 º (5.276) φLJ ( r ) = 4ε LJ «¨ ¸ −¨ r ¸ » © ¹» «© r ¹ ¬ ¼ where LJ and LJ are the Lennard-Jones length and energy parameters and r is the distance between two molecules. The system obtains equilibrium by raising the temperature from an initially low value to the desired equilibrium temperature by velocity rescaling, so that the system kinetic energy follows the relation 1 N23 m ci = Nk BT (5.277) 2 i =1 2 where ci is the speed of molecule i and N is the number of molecules in the system. A flux boundary at z = L is used to adjust the pressure to an equilibrium value. For ideal gases, the vapor boundary molecular flux J is related to pressure p and temperature T by kinetic theory relation k BT § p · J= (5.278) ¨ ¸ 2π m © k BT ¹ ¦ where m is the molecular mass of the injected species. Figure 5.33 shows the variation of pv, /psat with film thickness, which was predicted by the MDS for 2.0 pv ,δ psat 1.0 0.8 conventional 0.6 disjoining pressure 0.4 Molecular Dynamic Simulation 0.2 0 1 2 3 4 f 5 (nm) 6 Liquid film thickness, Figure 5.33 Comparison of molecular dynamic simulation results with conventional theory (Carey and Wemhoff, 2005). Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 405 thickness less than 6 nm. Also shown is the variation predicted by the conventional disjoining pressure theory, presented in Section 2.6.4. It can be seen that the MDS results generally agree with tabulated bulk saturation pressure within about 10% for film thickness greater than 1.5 nm, with the exception of 4.9 nm. The thermodynamic model only offers a prediction of the relationship between vapor pressure and film thickness at a given temperature. The molecular dynamics simulation, however, provides detailed information on the density and pressure variation in the thin liquid film and the variation of the equilibrium vapor pressure with temperature and thickness. Shi et al. (2005) used molecular dynamic simulations to analyze the contact angles of water droplets on a platinum surface over a range of temperature. To obtain high accuracy for the calculation of long-range forces, the P3M and Ewald summation methods were used. They found that the contact angle decays as the system temperature or the solid wall potential increases. In high temperature regions the contact angle decays very rapidly. The trends of the contact angle for water-platinum, argon-virtual solid wall and water-aluminum cases are similar. References Benjamin, T.B., 1957, “Wave Formation in a Laminar Flow Down an Inclined Plane,” Journal of Fluid Mechanics, Vol. 16, pp. 554-574. Brackbill, J.U., Kothe, D.B., and Zemach, C., 1992, “A Continuum Method for Modeling Surface Tension,” Journal of Computational Physics, Vol. 100, pp. 335-354. Brauer, H., 1956, “Stromung and Warmenbergang bei Rieselfilmen,” VDI Forschunsheft, Vol. 457, Dusseldorf. Buffone, C., Sefiane, K., and Christy, J.R.E., 2004, “Experimental investigation of the hydrodynamics and stability of an evaporating wetting film placed in a temperature gradient,” Applied Thermal Engineering, Vol. 24, pp. 1157-1170. Carey, V.P., 1992, Liquid-Vapor Phase-Change Phenomena: An Introduction to the Thermophysics of Vaporization and Condensation Processes in Heat Transfer Equipment, Hemisphere Publishing Corp., Washington, D. C. Carey, V.P., and Wemhoff, A. P., 2005, “Disjoining Pressure Effects in UltraThin Liquid Films in Micropassages – Comparison of Thermodynamic Theory with Predictions of Molecular Dynamics Simulations,” IMECE2005-80234, Proceedings of 2005 ASME International Mechanical Engineering Congress and Exposition, Orlando, FL. Derjaguin, B.V., 1955, “Definition of the Concept of and Magnitude of the Disjoining Pressure and its Role in the Statics and Kinetics of Thin Layers of Liquid,” Kolloidny Zhurnal, Vol. 17, pp. 191-197. 406 Transport Phenomena in Multiphase Systems Derjaguin, B.V., 1989, Theory of Stability of Colloids and Thin Films, Plenum, New York, NY. Faghri, A., 1995, Heat Pipe Science and Technology, Taylor & Francis, Washington, DC. Faghri, A., and Seban, R.A., 1985, “Heat Transfer in Wavy Liquid Films,” International Journal of Heat and Mass Transfer, Vol. 28, pp. 506-508. Faghri, A. and Payvar, P., 1979, “Transport to Thin Falling Liquid Films,” Reg. Journal of Energy Heat and Mass Transfer, Vol. 1, pp. 153-173. Glimm, J., Li, X.L., Liu, Y., and Zhao, N., 2001, “Conservative front tracking and level set algorithms,” Proceedings of National Acedemic of Science: Applied Mathemitics, Vol. 98, pp. 14198-14201. Harlow, F.H., and Welch, J.E., 1965, “Numerical Calculation of Time-Dependent Viscous Incompressible Flow,” Physics of Fluids, Vol. 8, pp. 2182-2189. Hirt, C.W. and Nichols, B.D., 1981, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” Journal of Computational Physics, Vol. 39, pp. 201-225. Holm, F.W., and Goplen, S.P., 1979, “Heat Transfer in the Meniscus Thin-Film Transition Region,” ASME Journal of Heat Transfer, Vol. 101, No. 3, pp. 543547. Jamet, D., Torres, D. and Brackbill, J.U., 2002, “On the Theory and Computation of Surface Tension: The Elimination of Parasitic Currents through Energy Conservation in the Second-Gradient Method,” Journal of Computational Physics, Vol. 182, pp. 262-276. Kapitza, P.L., 1964, “Wave Flow of Thin Layers of a Viscous Fluid,” Collected Papers of P.L. Kapitza, MacMillan, NY. Khrustalev, D., and Faghri, A., 1994, “Thermal Analysis of a Micro Heat Pipe,” ASME Journal of Heat Transfer, Vol. 116, No. 1, pp. 189-198. Khrustalev, D. K., and Faghri, A., 1995, “Heat Transfer during Evaporation and Condensation on Capillary-Grooved Structures of Heat Pipes,” ASME Journal of Heat Transfer Vol. 117, pp. 740-747. Koizumi, Y., Enari, R., and Ohtake, H., 2005, “Correlations of Characteristics of Waves on a Film Falling Down on a Vertical Wall,” Proceeding of the 2005 ASME International Mechanical Engineering Congress and Exposition, Nov. 511, Orlando, FL. Kucherov, R.Y., and Rikenglaz, L.E., 1960, “The Problem of Measuring the Condensation Coefficient,” Doklady Akad. Nauk. SSSR, Vol. 133, pp. 11301131. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 407 Kwok, D.Y., and Neumann, A. W., 1999, “Contact Angle Measurement and Contact Angle Interpretation,” Advances in Colloid and Interfaces Science, Vol. 81, pp. 167-249. Lide, D.R. ed., 2004, CRC Handbook of Chemistry and Physics, 85th Ed., CRC Press, Boca Raton, FL. Lin, L., and Faghri, A., 1999, “Heat Transfer in the Micro Region of a Rotating Miniature Heat Pipe,” International Journal of Heat and Mass Transfer, Vol. 42, pp. 1363-1369. Matsumoto, M., 1998, “Molecular Dynamics of Fluid Phase Change,” Fluid Phase Equilibria, Vol. 144, pp. 307-314. Mills, A.F. 1965, The Condensation of Steam at Low Pressures, Report No. NSF GP-2520, Series No. 6, Issue No. 39, Space Sciences Laboratory, University of California at Berkeley. Nichols, B.D., and Hirt, C.W., 1975, “Methods for Calculating Multidimensional, Transient Free Surface Flows Past Bodies,” Proc. of the First International Conf. On Num. Ship Hydrodynamics, Gaithersburg, MD, Oct. 2023. Nosoko, T., Yoshimura, P.N., Nagata, T., and Oyakawa, K., 1996, “Characteristics of Two-Dimensional Waves on a Falling Film,” Chemical Engineering Science, Vol. 51, pp. 725-732. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC. Paul, B., 1962, “Compilation of Evaporation Coefficients,” ARSJ, Vol. 32, pp. 1321-1328. Potash, M., and Wayner, P.C., 1972, “Evaporation from a Two-Dimensional Extended Meniscus,” International Journal of Heat and Mass Transfer, Vol. 15, pp. 1851-1863. Rice, J., and Faghri, A., 2005a, “A New Computational Method to Track a Liquid/Vapor Interface with Mass Transfer, Demonstrated on the Concentration Driven Evaporation from a Capillary Tube, and the Marangoni Effect,” Proceeding of the 2005 ASME International Mechanical Engineering Congress and Exposition, Nov. 5-11, Orlando, FL. Rice, J., and Faghri, A., 2005b, “A New Computational Method for Free Surface Problems,” Proceeding of the 2005 ASME Summer Heat Transfer Conference, July 17-22, San Francisco, CA. Rogovan, I.A., Olevskii, V.M., and Runova, N.G., 1969, “Measurement of the Parameters of Film Type Wavy Flow on a Vertical Plate,” Theoretical Foundations of Chemical Engineering, Vol. 3, pp. 164-171. 408 Transport Phenomena in Multiphase Systems Schrage, R.W., 1953, A Thermal Study of Interface Mass Transfer, Columbia University Press, New York. Shi, B., Sinha, S., Dhir, V.K., 2005, “Molecular Simulation of the Contact Angle of Water Droplet on a Platinum Surface,” IMECE2005-80136, Proceedings of 2005 ASME International Mechanical Engineering Congress and Exposition, Orlando, FL. Shopov, P.J., Minev, P.D., Bazhekov, I.B., and Zapryanov, Z.D., 1990, “Interaction of a deformable bubble with a rigid wall at moderate Reynolds numbers,” Journal of Fluid Mechanics, Vol. 219, pp. 241-271. Stepanov, V.G., Volyak, L.D., and Tarlakov, Y.V., 1977, “Wetting Contact Angles for Some Systems,” Journal of Engineering Physics and Thermophysics, Vol. 32, No. 6, pp. 1000-1003. Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S., and Jan, Y., 2001, “A Front-Tracking Method for the Computations of Multiphase Flow,” Journal of Computational Physics, Vol. 169, pp. 707-759. Wu, Y. W., and Pan, C., 2003, “A Molecular Dynamics Simulation of Bubble Nucleation in Homogenous Liquid under Heating with Constant Mean Negative Pressure,” Microscale Thermophysical Engineering, Vol. 7, pp. 137-151. Yang, J., Han, J., Isaacson, K., and Kwok, D. Y., 2003, “Effects of Surface Defects, Polycrystallinity, and Nanostructure of Self-Assembled Monolayers for Octadecanethiol Adsorbed on to Au on Wetting and its Surface Energetic Interpretation,” Langmuir, Vol. 19, pp. 9231-9238. Yang, T.H., and Pan, C., 2005, “Molecular Dynamics Simulation of a Thin Water Layer Evaporation and Evaporation Coefficient,” International Journal of Heat and Mass Transfer, Vol. 48, pp. 3516-3526. Yasuoka, K., and Matsumoto, M., 1995, “Molecular simulation of evaporation and condensation. I. Self condensation and molecular exchange,” Proceedings of ASME/JSME Thermal Engineering Conference, pp. 459-464. Yaws, C.L., 1992, Thermodynamic and Physical Property Data, Gulf Pub. Co., Houston, Texas. Youngs, D.L., 1982, “Time-Dependent Multimaterial Flow with Large Fluid Distortion,” Numerical Methods for Fluid Dynamics, K.W. Morton and M.J. Baines, eds. Academic Press, pp. 273-285. Zhang, Y., and Faghri, A., 2002, “Heat Transfer in a Pulsating Heat Pipe with Open End,” International Journal of Heat Mass Transfer, Vol. 45, No. 4, pp. 755-764. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 409 Problems 5.1. Liquid droplets that condense to form fog in the atmosphere can be as small as 2 m in diameter when they first form. Determine the pressure in a droplet of this size at 20 °C (101.3 kPa). Note that for water in contact with air, the surface tension σ is 0.0728 N/m. Describe any assumption you need to make in order to get the final answer. 5.2. A 0.2-mm-thick water film sits on a surface held at a temperature of 60 °C. The liquid film is exposed to air at a bulk temperature of TG=20 °C, and the convective heat transfer coefficient between the liquid film and the air is h ˆ = 15 W/m2K. Suppose the wave number is α = 2 . What is the critical Marangoni number, Mac ? Is the liquid film stable? 5.3. A 0.2-mm-diameter tube is vertically placed in a pool of water at 20 °C. Find the capillary rise in tubes of the following materials: (a) aluminum, (b) brass, (c) copper, and (d) steel. 5.4. Two nearly vertical but nonparallel plates touch a pool of water with their parallel bottom edges. The angle between the plates is 2γ so that these plates would intersect somewhere under the pool surface. The distance between the plates at the pool surface level, 2W, is small. The apparent contact angle is denoted by θ . The liquid-vapor meniscus between the plates is elevated from the pool surface level due to capillary pressure. Assume that the meniscus curvature does not change along the liquid-vapor interface. Derive an algebraic equation for the height of the capillary rise H for this situation (for θ ≠ 0 and γ ≠ 0 ). Using the derived equation, estimate H for the following data: W = 0.5 mm, θ = 0 , γ = 0 , ρA = 1000 kg/m3 , σ = 0.06 N/m. 5.5. A small amount of liquid metal resides between two high melting point powder particles separated by a distance of A (A > 2r0 ) (see Fig. P5.1). If the effect of gravity on the liquid metal is negligible, estimate the force required to hold the two particles together. Figure P5.1 410 Transport Phenomena in Multiphase Systems 5.6. Write the continuity, momentum, energy, and species equations and the necessary boundary conditions for a very long capillary tube that is open at both ends, as shown in Fig P5.2. The evaporation is driven by the concentration gradient of vapor in the air. The evaporation cools the interface while the wall heats the fluid, causing a temperature gradient along the interface. The assumptions are the same as Example 5.4. ′′ qconv Vapor/Air Axis Liquid Figure P5.2 5.7. Consider an ultrathin evaporating liquid film on a hot solid surface of temperature Tw in the presence of a saturated vapor of temperature Tv. The film thickness δ is small and the disjoining pressure pd is significant. The film is flat. Does the disjoining pressure increase or decrease the evaporation rate in comparison to the hypothetical case where pd = 0 (for the same film thickness and temperatures)? Assuming that Tw − Tv is extremely small, explain the effect of pd on the evaporation rate by referring to the extended Kelvin equation and interfacial resistance. 5.8. A liquid film of ammonia with thickness δ = 0.375 × 10−8 m is evaporating from a smooth solid surface into bulk vapor. The solid surface temperature is Tw = 250.12 K. The bulk vapor temperature is Tv = 250 K. The constants in the correlation for disjoining pressure pd = − A′δ − B are A′ = 10−21 and B = 3. The accommodation coefficient is equal to unity. In this particular case, the perfect gas law can be used instead of the saturation tables, for convenience. Estimate the heat flux at the liquid interface. 5.9. A slug of water 4.5 mm in length sits in a vertical copper tube with an inner radius of 2 mm at room temperature. The temperature of the whole system is slowly increasing. At what temperature will the water slug start to move? 5.10. Derive eq. (5.159) from the momentum balance at the liquid interface. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 411 5.11. Starting from the appropriate equations of Chapter 3, derive eq. (5.160). Show the required control volumes. 5.12. Nusselt-type condensation takes place on a flat vertical plate made from aluminum. The thickness of the plate is 5 mm and there are 1 mm ID micro channels inside it containing cooling liquid. You must increase the rate of condensation, which is restricted primarily by the conduction through the condensate film. What would you do? 5.13. The effective heat transfer during evaporation of liquid is more extensive from a solid surface with small grooves than from a plain surface. How would you explain this? 5.14. A vertical capillary tube is connected to bulk liquid at the bottom end and to a heat source at the upper end. The liquid evaporates from the liquidvapor meniscus, so a quasi-steady state exists. What is the heat load that would initiate dry-out of the upper end and lead to unbounded increase of its temperature? Write corresponding general formulas or equations. 5.15. The temperature at the meniscus of Problem 5.14 equals the saturation temperature, and the bulk liquid is at a lower temperature Tb. Repeat Problem 5.14 accounting for the conduction heat loss from the meniscus to the bulk liquid. 5.16. A capillary tube of inner radius r is vertically dipped into a liquid. The depth from the free surface to the lower end of the tube is a. When the liquid rises in the capillary tube due to capillary force, it can be viewed as a variable mass system with capillary force, gravity, and viscous force acting upon it. What is the governing equation for the liquid slug in the capillary tube? The shape of the meniscus in the tube can be assumed to be the same at all times. 5.17. Heat transfer in the thin-film region of an axially-grooved evaporator was analyzed in Section 5.5.4. For the condenser region, the effect of interfacial resistance becomes negligible so that the temperature of the liquid vapor interface is equal to the vapor temperature ( Tδ = Tv ). It is further assumed that (1) the surface of the liquid film is smooth and the film thickness variation along the s-coordinate is weak (see Fig. 5.15), i.e., (d δ / ds ) 2  1 , and (2) the disjoining pressure gradient along the film can be neglected in comparison to that of the capillary pressure because of the large film thickness. Analyze the heat transfer on the condenser region based on the above assumptions. 412 Transport Phenomena in Multiphase Systems (a) (c) Figure P5.3 (b) 5.18. A typical horizontally oriented rotating miniature heat pipe is shown in Fig. P5.3(a). The fluid flow in the rotating miniature heat pipe is characterized by axial fluid flows and liquid film flows transverse to the axial direction. Axial fluid flows include an axial liquid flow in the groove and a countercurrent vapor flow with varied mass rate due to evaporation and condensation. Liquid film flows transverse to the axial direction occur both in the evaporator and in the condenser. In the evaporator, an evaporating thin liquid film flowing up the groove side wall is driven by the disjoining pressure and the surface tension. In the condenser, a condensate film flowing into the groove is sustained by centrifugal force. A change in radius of curvature of the liquid-vapor interface causes a capillary pressure difference between the condenser and evaporator, which promotes the flow of condensate back to the evaporator. In addition, a varying liquid depth allows the centrifugal acceleration to produce a hydrostatic pressure change along the groove that pumps the condensate back to the evaporator. In order to analyze heat transfer in the micro region, the radius of curvature of the liquid-vapor interface in the meniscus region, as shown in Fig. P5.3 (b) is assumed to be constant. The coordinate system for the evaporating thin film flow is shown in Fig. P5.3(c). The assumptions made are as follows: The vapor temperature, Tv, and the wall temperature, Tw, are constant; the radius of curvature of the liquid-vapor interface in the meniscus region and the rotational radius of the micro region, rmic, are constant, and the influences of the Coriolis force, gravitational force, and Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 413 vapor drag on the evaporating thin film flow are negligible. Analyze heat transfer in the evaporator of the rotating miniature heat pipe. 5.19. The cross-sectional view and the coordinate system of the condenser section of a rotating miniature heat pipe are shown in Fig. P5.4. It is assumed that the vapor and liquid flows are incompressible and that the vapor and liquid are in a saturated state. The contact angle of the meniscus, , is constant. The effect of gravitational force on the flow is negligible, and the influences of the Coriolis force and surface tension on the circumferential condensate film flow are negligibly small. Thus, the circumferential film flow is symmetrical about the rotational coordinate of y. Assuming that the condensate film thickness is very small compared to sb, the momentum equations in the y and x directions for the condensate film flow and boundary conditions are similar to those used for rotating heat pipes with a large diameter. Derive the ordinary differential equation that governs the liquid film thickness. Figure P5.4 5.20. A Newtonian liquid is in an open cylindrical container with a radius of R and height of H (see Fig. P5.5). A stationary continuous axisymmetric laser beam with an intensity of I = I 0 exp(−r 2 / r02 ) irradiates the top of the liquid surface. The side of the cylindrical container is maintained at a constant temperature, T0 , and the bottom of the container is adiabatic. Supposing the surface tension is a linear function of temperature σ = C0 − C1T , write the governing equations and the corresponding boundary conditions of the problem. The effect of gravity is negligible, and the liquid flow induced by the laser beam can be assumed to be axisymmetric. 414 Transport Phenomena in Multiphase Systems Figure P5.5 5.21. Introduce appropriate nondimensional variables and nondimensionalize the governing equations and boundary conditions in Problem 5.20. Identify the dominant dimensionless variables in the problem. 5.22. Saturated vapor enters a miniature channel with a radius of R and a length of L (see Fig. P5.6). Condensation takes place on the wall of the channel, since the wall temperature Tw is below the saturation temperature of the working fluid Tsat. The condensate fluid flows in the positive x-direction due to the effects of shear force and surface tension. The problem is condensation on the inner surface of the miniature channel with cocurrent vapor flow. The average vapor velocity is not constant along the xdirection because condensation occurring on the wall reduces the amount of vapor flow in the core of the tube. Capillary blocking occurs when capillary force causes the liquid to block the tube cross-section at some distance, Lδ , from the condenser entrance. Supposing the effect of gravity on the two-phase flow in the miniature tube is negligible, write the governing equations and boundary conditions for the liquid and vapor flows. Figure P5.6 Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 415 5.23. Heat is applied to the outer surface of a microchannel with inner and outer radii of ri and ro , respectively, as shown in Fig. P5.7. Specify the governing equations that describe evaporation in the microchannel. The effects of disjoining pressure and interfacial thermal resistance must be taken into account. Figure P5.7 Figure P5.8 5.24. Evaporation in the evaporator section of a pulsating heat pipe (PHP) with an open end is analyzed in Example 5.6. The physical model of the condenser section of the PHP with an open end is shown in Fig. P5.8. It is assumed that the liquid film in the condenser section is divided into two regions: a thin film region and a meniscus region. These two regions are separated by a transition point where the film thickness is δ = δ tr . The curvature of the liquid film in the meniscus region is assumed to be constant. Assuming the effect of disjoining pressure is negligible, specify the governing equation for the liquid film thickness in 0 < x < L2 . 5.25. The fluid mechanics of capillary tubes and microchannels is of great interest for modeling transport phenomena in micro devices. Develop the physical formulation, including continuity and momentum equations with appropriate boundary and initial conditions, for a typical miniature channel with a free surface as depicted in Fig. P5.9 for both submerged and blocked end configuration. Neglect mass transfer and Marangoni effects. Describe the physical/mathematical formulation to model the problem related to instability growth of a free-flowing cylindrical column of fluid as shown in Fig. P5.10. 5.26. 416 Transport Phenomena in Multiphase Systems Open to exterior environment Vapor /gas mixture Parallel plates or cylindrical tube g (if applicable) 2y 0 Liquid x,u Submerged or blocked end y, v Figure P5.9 Vapor or vapor /gas mixture 2y 0 y,v x, Location of initial disturbance Liquid Initial column of liquid Figure P5.10 Disturbed column of liquid, f(x) 5.27. Obtain an analytical solution for the shape of the meniscus in a capillary tube in 1 g environment for mechanical equilibrium, as shown in Fig. P5.11 for both parallel plates and circular tubes. Assume the effect of the gaseous region is negligible with vapor pressure constant. You should get a second order ordinary differential equation for f ( y ) which can be solved by the Runge-Kutta method. Obtain numerical results for the working fluid as acetone with y0 = 0.001 m and contact angle θ = 0.0 . Obtain a closed analytical solution for the situation with no gravity. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 417 Vapor Interface, f(y) (meniscus) g Parallel plates or cylindrical tube 2y0 Liquid x,u y,v Figure P5.11 Vapor (meniscus) g 2y0 Liquid H x,u Submerged in liquid y,v Liquid Figure P5.12 5.28. Obtain an analytical solution for the rise of a meniscus in a capillary tube with both ends open, as shown in Fig. P5.12, for both parallel plates and a circular tube. If u is the bulk motion or the average velocity, obtain your solution in terms of u = dH / dt where H is the mean height of the meniscus. 418 Transport Phenomena in Multiphase Systems z y w(y) yz g Figure P5.13 5.29. Consider the steady, laminar flow of a liquid film down an inclined flat plate (see Fig. P5.13). Assume the flow is fully developed. However, the presence of a small temperature gradient in the axial direction dT/dz produces a constant surface tension gradient d /dz = C. Assume the surface tension effect is only at the interface, which produces a nonzero shear stress and no other effect on physical properties of the system. Solve the momentum equation to determine the shear stress variation and velocity distribution as a function of distance perpendicular to the plate. Consider the absorption of gas A (O2) to a laminar, steady, fully developed, isothermal falling film of liquid B (water) (see Fig. 5.14). Since O2 is slightly soluble in water, one can assume the properties of water are not affected. Since diffusion of O2 is taking place slowly and thermal diffusion is in the immediate vicinity of liquid, O2 does not penetrate very far into the liquid. Obtain the total mole flow rate of gas A to liquid B. Consider the solid dissolution of a component A for a solid wall to a steady, laminar, isothermal falling liquid film B (see Fig. 5.15). Assume that species A is slightly soluble in B. Obtain the physical formulation of the simplified model. Saturated vapor at 100 °C flows inside a horizontal tube cooled at the external surface. A pool of the liquid water exists at the bottom of the tube and flows with a velocity of 0.5 m/s. Determine the critical velocity of the vapor when waves appear on the liquid-vapor interface. 5.30. 5.31. 5.32. Chapter 5 Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer 419 wall y z water g O2 y z liquid g A B ′′ nO2 Figure P5.14 Figure P5.15 5.33. In a stratified, horizontal, two-phase flow system, the saturated water vapor flows above the saturated liquid water. The pressure of the two-phase system is at 5 bar. What is the most dangerous wavelength, λD ? 420 Transport Phenomena in Multiphase Systems