- Preface
- Introduction to Transport Phenomena
- Thermodynamics of Multiphase Systems
- Generalized Governing Equations in Mutliphase Systems: Local Instance Formulations
- Generalized Governing Equations for Multiphase Systems : Averaging Formulations
- Solid-Liquid-Vapor Phenomena and Interfacial Heat and Mass Transfer
- Melting and Solidification
- Sublimation and Vapor Deposition
- Condensation
- Evaporation
- Boiling
- Two-Phase Flow and Heat Filter
- Appendix A : Constants, Units and Conversion Factors
- Appendix B: Transport Properties
- Appendix C: Vectors and Tensors
- Index

4
GENERALIZED GOVERNING EQUATIONS FOR MULTIPHASE SYSTEMS: AVERAGING FORMULATIONS
4.1 Introduction
In order to use the governing equations and the jump conditions derived in Chapter 3, one must explicitly track the interfaces in a multiphase system. However, this is not always possible, especially for cases with multiple interfaces. Many multiphase-flow problems encountered in engineering – such as those with dispersed and mixed phases in Table 1.10 – have extremely complicated and deformable interfaces. It is not always possible to solve the local instance fluid flow, because the difficulty associated with interface tracking exceeds present computational capability. Fortunately, information about the discontinuity of properties at the interfaces and the exact locations of the interfaces are not always of interest to practical engineers. The macroscopic aspects of multiphase flow are more important to the design and operation of a multiphase system. Appropriate averaging can obtain the mean values of flow and thermal properties and eliminate the need to explicitly track interfaces and/or the local instance fluctuations of properties. Either multifluid or homogeneous models based on spatial- and/or timeaveraging techniques can be employed (Ishii, 1975); therefore, microscopically detailed interfaces between phases need not be explicitly identified. Various averaging techniques, including time, spatial, and area averaging, are used to obtain nondimensional parameters that correlate the experimental data as well as flow maps for two-phase flow. The local instance formulation is the fundamental basis for all two-phase flow models, regardless of the averaging method employed. When each subregion is bounded by interfaces that can be considered continuous, the local instance formulation is mathematically rigorous. Consequently, all two-phase flow models should start from this formulation and proceed by applying appropriate averaging methods. For instance, if spatial
238 Transport Phenomena in Multiphase Systems
averaging is performed for each individual phase within a multiphase control volume, the multifluid model is obtained. If averaging is performed over all phases in the control volume, the homogeneous model is obtained. If averaging is performed over the cross-section of a channel, the area-averaging model is obtained and the governing equations are reduced to a one-dimensional formulation. Section 4.2 presents a thorough overview of various averaging methods used to describe multiphase flow and heat and mass transfer problems; this is followed by the governing equations for the multifluid and homogeneous models in Sections 4.3 and 4.4, respectively. The area-averaged governing equations for one-dimensional flows are discussed in Section 4.5. Single- and multiphase transport phenomena in porous media are presented in Section 4.6. Chapter 4 closes with a discussion of the Lattice Boltzmann model.
4.2 Overview of Averaging Approaches
The objectives of the various averaging methods are twofold: (1) to define the average properties for the multiphase system and correlate the experimental data, and (2) to obtain solvable governing equations that can be used to predict the macroscopic properties of the multiphase system. The application of averaging methods to the definition of two-phase flow properties will be discussed in Chapter 11. This chapter will address the application of averaging methods to the governing equations. Based on the physical concepts used to formulate multiphase transport phenomena, the averaging methods can be classified into three major groups: (1) Eulerian averaging, (2) Lagrangian averaging, and (3) Molecular statistical averaging. These averaging techniques are briefly reviewed below.
4.2.1 Eulerian Averaging
Eulerian averaging is the most important and widely-used method of averaging, because it is consistent with the control volume analysis that we used to develop the governing equations in the preceding section. It is also applicable to the most common techniques of experimental observations. Eulerian averaging is based on time-space description of physical phenomena. In the Eulerian description, changes in the various dependent variables, such as velocity, temperature, and pressure, are expressed as functions of time and space coordinates, which are considered to be independent variables. One can average these independent variables over both space and time. The integral operations associated with these averages smooth out the local spatial or instant variations of the properties within the domain of integration. For a generalized function Φ = Φ ( x, y, z , t ) , the most widely-used Eulerian averaging includes time averaging and volumetric averaging. The Eulerian time
Chapter 4 Generalized Governing Equations: Averaging Formulations 239
average is obtained by averaging the flow properties over a certain period of time, Δt, at a fixed point in the reference frame, i.e., 1 Φ = ³Δt Φ ( x, y, z , t )dt (4.1) Δt for this equation, the time period Δt is chosen so that it is larger than the largest time scale of the local properties’ fluctuation, yet small enough in comparison to the process macroscopic time scale of the process. During this time period, different phases can flow through the fixed point. Eulerian time averaging is particularly useful for a turbulent multiphase flow as well as for the dispersed phase systems summarized in Table 1.10. Eulerian volumetric averaging is usually performed over a volume element, ΔV, around a point ( x, y, z ) in the flow. For a multiphase system that includes Π different phases, the total volume equals the summation of the individual phase volumes, i.e.,
ΔV = ¦ ΔVk
k =1
Π
(4.2)
The volume fraction of the kth phase, ε k , is defined as the ratio of the elemental volume of the kth phase to the total elemental volume for all phases, i.e., ΔV εk = k (4.3) ΔV The volume fraction of all phases must sum to unity:
k =1
¦εk
Π
=1
(4.4)
Eulerian volume averaging is expressed as 1Π Φ= (4.5) ¦ ³ Φ k ( x, y, z, t )dV ΔV k =1 ΔVk where the volume element ΔV must be much smaller than the total volume of the multiphase system so that the average can provide a local value of Φ in the flow field. The volume element ΔV must also be large enough to yield a stationary average. Since the volume element includes different phases, information about the spatial variation of Φ for each individual phase is lost and Φ represents the average for all phases.
For any variable or property that is associated with a particular phase, Φ k , the phase-average value of any variable or property for that phase is obtained with the following equations
Intrinsic phase average:
Φk
k
=
1 ΔVk
³ΔV
k
Φ k dV
(4.6)
Extrinsic phase average:
240 Transport Phenomena in Multiphase Systems
1 (4.7) ³ Φ k dV ΔV ΔVk Intrinsic means that it forms to the inherent part of a phase and is independent of other phases in the volume element. In contrast, extrinsic means it is a property that depends on the phase’s relationship with other phases in the volume element. While the intrinsic phase average is taken over only the volume of the kth phase in eq. (4.6), the extrinsic phase average for a particular phase is taken over an entire elemental volume in eq. (4.7). These two phase-averages are related by k Φk = ε k Φk (4.8) The intrinsic and extrinsic phase averages defined in eqs. (4.6) and (4.7) are related to the volume average defined in eq. (4.5) by Φk =
Φ = ¦ Φk = ¦ ε k Φk
k =1 k =1
Π
Π
k
(4.9)
The deviation from a respective intrinsic phase-average value is k ˆ Φk ≡ Φk − Φk (4.10) When the products of two variables are phase-averaged, the following relations are needed: k k k k ˆˆ Φk Ψ k = Φk Ψ k + Φk Ψ k (4.11)
Φk Ψ k = ε k Φk
k
Ψk
k
ˆˆ + Φk Ψk
(4.12)
In order to obtain the volume-averaged governing equations, the volume average of the partial derivative with respect to time and gradient must be obtained. For a control volume ΔV shown in Fig. 4.1, the volume averaging of
ΔV
nk
dAk
Figure 4.1 Control volume for volume averaging.
Chapter 4 Generalized Governing Equations: Averaging Formulations 241
the partial derivative with respect to time is obtained by the following general transport theorem: ∂ Ωk ∂Ω k 1 = − (4.13) ³ Ωk VI ⋅ n k dAk ∂t ∂t ΔV Ak where Ak is the is the interfacial area surrounding the kth phase within control volume ΔV , ΔVk is the volume occupied by the kth phase in the control volume and ΔV , VI is the interfacial velocity, and nk is the unit normal vector at the interface directed outward from phase k (see Fig. 4.1). The volume average of the gradient is 1 ∇Ω k = ∇ Ω k + (4.14) ³ Ω k n k dAk ΔV Ak and the volume average of a divergence is 1 ∇ ⋅ Ωk = ∇ ⋅ Ωk + (4.15) ³ Ω k ⋅ n k dAk ΔV Ak The general quantity Ω k in eqs. (4.13) and (4.14) can be a scalar, vector, or tensor of the second order. It can be a vector or tensor of the second order in eq. (4.15). The formulation of macroscopic equations for multiphase systems can be classified into two groups: (1) the multi-fluid model (Section 4.3), and (2) the homogeneous model (Section 4.4), also known as the mixture or diffuse model. If the averaging is performed for each individual phase within a multiphase control volume, as shown in eqs. (4.6) and (4.7), one obtains the multifluid model, in which Π sets of averaged conservation equations – each set includes continuity, momentum and energy equations – describe the flow of a Π − phase system. The equations will also include source terms that account for the transfer of momentum, energy, and mass between phases. If only two phases are present, the multifluid model is referred to as the two-fluid model. However, if spatial averaging is performed over both phases simultaneously within a multiphase control volume, as indicated in eq. (4.5), the homogeneous model is obtained; in this case the mixture of a two-phase fluid would be considered a whole. The governing equations for the homogeneous model comprise a single set of equations including continuity, momentum, and energy equations, with one additional diffusion equation to account for the concentration change due to interphase mass transfer by phase change. Continuity, momentum, and energy equations for the mixture model can be obtained by adding together the governing equations for the multifluid models; a diffusion model must be developed to account for mass transfer between phases. In this chapter, it is assumed for the sake of simplicity that the reference frame is stationary.
242 Transport Phenomena in Multiphase Systems
4.2.2 Lagrangian Averaging
Lagrangian averaging is directly related to the Lagrangian description of a system, which requires tracking the motion of each individual fluid particle. Therefore, Lagrangian averaging is a very useful tool when the dynamics of individual particles are of interest. To obtain Lagrangian time averaging, it is necessary to follow a specific particle and observe its behavior for a certain time interval. Then, the behavior of this particle is averaged over the time interval. For a generalized function Φ = Φ ( X , Y , Z , t ) , X, Y, and Z are material coordinates moving with the particle, and X, Y, Z are functions of the spatial coordinates x, y, z, and time t, i.e., X = X ( x, y, z , t ), Y = Y ( x, y, z , t ), Z = Z ( x, y, z , t ) The most widely used Lagrangian averaging is time averaging, where the time average of the function Φ in time interval of Δt is 1 Φ = ³Δt Φ ( X , Y , Z , t )dt (4.16) Δt Lagrangian time averaging is performed for a distinct particle moving in the field; therefore, X, Y, and Z in the time interval Δt are not fixed in space. This focus on specific particles moving in space and time distinguishes Lagrangian averaging from Eulerian time averaging, which treats a fixed point in space relative to the reference frame. An example from daily experience will serve to illustrate this difference. In order to monitor traffic on the highway, the speed of all cars passing a point can be measured and averaged over a certain time interval – a case of Eulerian averaging. To catch an individual speeder, the police must follow the vehicle of interest to measure its speed as it moves in space over a certain time interval – a case of Lagrangian averaging.
4.2.3 Molecular Statistical Averaging
When the collective mechanics of a large number of particles is of interest, molecular statistical averaging may be employed. This relies on the concept of particle number density, which is the number of particles per unit volume. For a system with a large number of particles, the behavior of each individual particle is random because random collisions occur. To describe the behavior of each particle, it is necessary to track the motion resulting from each collision – an impractical and often unnecessary task. Although the behavior of each particle is random, the collection of particles may demonstrate some statistical behaviors that are different from those of the individual particles. When the number of molecules involved in the averaging process is large enough, the statistical average value becomes independent of the number of molecules involved. The statistical average value of the microscopic properties for a large number of molecules is related to the macroscopic properties of the system. For example, temperature is a statistical measure of the kinetic energy of the individual molecules, and the pressure of a gas in a container is the result of many molecules’ collisions with the wall. For some engineering problems, the
Chapter 4 Generalized Governing Equations: Averaging Formulations 243
macroscopic properties of the fluid are of interest as well as the microscopic properties are required for design or analysis. For this reason, a detailed presentation of Boltzmann statistical averaging including the discussion of Lattice Boltzmann model for both single and multiphase systems will be presented in Section 4.7.
4.3 Volume-Averaged Multi-Fluid Models
If spatial averaging is performed for each individual phase within a multiphase control volume, the multifluid model is obtained. Additional source terms are needed in these equations to account for the interaction between phases.
4.3.1 Continuity Equation
The volume average of the continuity equation for the kth phase is obtained by taking extrinsic phase averaging from eq. (3.38) ∂ρ k (4.17) + ∇ ⋅ ρ k Vk = 0 ∂t where the two terms on the left-hand side can be obtained by using eqs. (4.13) and (4.15), i.e., ∂ ρk ∂ρ k 1 = − ³ ρk VI ⋅ n k dAk ∂t ∂t ΔV Ak 1 ∇ ⋅ ρ k Vk = ∇ ⋅ ρ k Vk + ³ ρ k Vk ⋅ n k dAk ΔV Ak Substituting the above expressions into eq. (4.17), the volume-averaged continuity equation becomes ∂ ρk 1 (4.18) + ∇ ⋅ ρ k Vk = − ³ ρk (Vk − VI ) ⋅ n k dAk ∂t ΔV Ak The right-hand side of eq. (4.18) represents mass transfer per unit volume from all other phases to the kth phase due to phase change; it can be rewritten as Π 1 − (4.19) jk ³Ak ρ k (Vk − VI ) ⋅ n k dAk = ¦ m′′′ ΔV j =1( j ≠ k ) where m′′′ represents mass transfer per unit volume from the jth to the kth phase jk due to phase change. The value of m′′′ depends on the phase change process jk that takes place in the multiphase system, and the conservation of mass requires ′′′ that m′′′ = − mkj . jk The extrinsic phase-averaged density, averaged density,
ρ k , is related to the intrinsic phase-
ρk
k
, by
244 Transport Phenomena in Multiphase Systems
ρk = ε k ρk
k
(4.20)
Furthermore, the intrinsic phase-averaged density is equal to the density ρ k . Substituting eqs. (4.19) and (4.20) into eq. (4.18), and considering eq. (4.12), the continuity equation for the kth phase becomes Π ∂ k k k ˆˆ ε k ρ k + ∇ ⋅ ε k ρ k Vk + ρ k Vk = ¦ m′′′ (4.21) jk ∂t j =1( j ≠ k ) ˆˆ The dispersive term in eq. (4.21), ρ V , is generally small compared
(
)
(
)
k
k
with ε k ρ k Vk ; it is assumed that it can be neglected. The continuity equation for the kth phase becomes Π ∂ k k k ε k ρ k + ∇ ⋅ ε k ρ k Vk = ¦ m′′′ (4.22) jk ∂t j =1( j ≠ k )
k
k
(
)
(
)
4.3.2 Momentum Equation
The extrinsic phase-averaged momentum equation for the kth phase can be obtained by performing extrinsic phase-averaging on the momentum equation (3.52): ∂ ( ρ k Vk ) (4.23) + ∇ ⋅ ( ρ k Vk Vk ) = ∇ ⋅ ′ + ρ k X k k ∂t where the body force per unit mass is assumed to be the same for different species for sake of simplicity. Each term in eq. (4.23) can be obtained using eqs. (4.13) – (4.15), i.e., ∂ ρ k Vk ∂ρ k Vk 1 = − ³ ρ k Vk (VI ⋅ n k )dAk ∂t ∂t ΔV Ak 1 ∇ ⋅ ρ k Vk Vk = ∇ ⋅ ρ k Vk Vk + ³ ρ k Vk Vk ⋅ n k dAk ΔV Ak 1 ∇⋅ ′ =∇⋅ ′ + k k ³ ′k ⋅ n k dAk ΔV Ak ρk Xk = ρk Xk Substituting the above expressions into eq. (4.23) yields ∂ ρ k Vk 1 + ∇ ⋅ ρ k Vk Vk = ∇ ⋅ ′ − k ³ ρ k Vk (Vk − VI ) ⋅ n k dAk ∂t ΔV Ak (4.24) 1 + ³ ′k ⋅ n k dAk + ρk Xk ΔV Ak Considering eq. (4.12) and ignoring the product of the deviations, eq. (4.24) becomes
Chapter 4 Generalized Governing Equations: Averaging Formulations 245
∂ k k k k k = ∇ ⋅ εk ′ ε k ρ k Vk + ∇ ⋅ ε k ρ k Vk Vk k ∂t 1 1 k − ³Ak ρ k Vk (Vk − VI ) ⋅ n k dAk + ΔV ³Ak ′k ⋅ n k dAk + ε k ρk Xk ΔV k ′ where is the phase-averaged stress tensor for the kth phase: k ′ k
k
(
)
(
)
(
)
(4.25)
2 k k Tº k ª I + μk «∇ Vk + ∇ Vk » − 3 μ k (∇ ⋅ Vk )I (4.26) ¬ ¼ The second and third terms on the right-hand side of eq. (4.25) represent the momentum exchanges and interactive forces between all other phases and the kth phase. Π k 1 − ρ k Vk (Vk − VI ) ⋅ n k dAk = ¦ m′′′ Vk , I (4.27) jk ³Ak ΔV j =1( j ≠ k ) Π 1 ′ ⋅ n k dAk = ¦ F jk (4.28) k ³ ΔV Ak j =1( j ≠ k ) = − pk
k
(
)
where
Vk , I
k
is intrinsic phase-averaged velocity of the kth phase at the
interface. The difference between two adjacent phases results solely from the density difference between the two phases. F jk is an interactive force between the jth and the kth phase, and depends on the friction, pressure, and cohesion between different phases. Newton’s third law requires that the interactive forces satisfy F jk = − Fkj (4.29) The interactive force can be determined by F jk = K jk
(V
j j
− Vk
k
)
(4.30)
where Kjk is the momentum exchange coefficient between phases j and k. Determining the momentum exchange coefficient is a very challenging task because interphase momentum exchange depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k, as is the case with the dispersed phase system summarized in Table 1.10, one can assume that the secondary phase is spherical in shape and an appropriate empirical correlation can be used to obtain the momentum exchange coefficient. Since liquid-vapor flow is widely used in various applications, we will use liquid-vapor flow as an example to explain the determination of the momentum exchange coefficient. If liquid is the primary phase and vapor is the secondary phase, the vapor phase is dispersed in the liquid as vapor bubbles. If vapor is the primary phase and liquid is the secondary phase, the liquid phase is dispersed in the vapor as liquid droplets. Boysan (1990) suggested that the momentum exchange coefficient could be estimated by
246 Transport Phenomena in Multiphase Systems
ε j ρk 3 K jk = CD dj 4
Vj
j
− Vk
k
(4.31)
where phase k is the primary phase, and phase j is the secondary phase, and dj is the diameter of vapor bubbles or liquid droplets of the secondary phase j. CD is the drag coefficient based on the relative Reynolds number, which obtained by the following empirical correlations: 24 0.687 Re ≤ 1000 ) ° (1 + 0.15Re CD = ® Re (4.32) Re > 1000 °0.44 ¯ where
Re =
ρk
Vj
j
− Vk
k
dj (4.33)
μk
Substituting eqs. (4.27) and (4.28) into eq. (4.25), the multifluid volumeaveraged momentum equation becomes ∂ k k k k ε k ρ k Vk + ∇ ⋅ ε k ρ k Vk Vk ∂t (4.34) Π k k k = ∇ ⋅ εk ′ + ε k ρk Xk + ¦ F jk + m′′′ Vk , I k jk
(
)
(
(
)
j =1( j ≠ k )
(
)
)
Example 4.1 A mixture of water and its vapor at 1 atm flows in a 0.1 m ID tube with a mass flow rate of 0.1 kg/s. The liquid water is dispersed in the vapor phase in the form of 0.1- mm -diameter droplets and the quality of the mixture is x=0.8. While the vapor is saturated, the liquid droplets, are subcooled, both at 95 °C. The volume fraction of the liquid phase is ε A = 0.01 . Find the interactive force between the liquid and vapor phases. Solution: Since the liquid droplets are dispersed in the vapor phase, vapor ( v ) is the primary phase and liquid ( A ) is the secondary phase. At 1 atm, the properties of the saturated vapor are μv = 16.698 × 10−6 N-s/m 2 ρv = 0.596 kg/m3 The density of the subcooled water at 95 °C is ρA = 965.35kg/m3 The mass flow rates of the liquid and vapor phases are, respectively, mv = xm = 0.8 × 0.1 = 0.08kg/s mA = (1 − x)m = (1 − 0.8) × 0.1 = 0.02kg/s The total cross-sectional area of the tube is 1 1 A = π D 2 = × π × 0.12 = 7.85 × 10−3 m 2 4 4 The mass flow rates of the liquid and vapor phases are, respectively,
Chapter 4 Generalized Governing Equations: Averaging Formulations 247
mA = ρA wA AA = ρA wA ε A A
mv = ρ v wv Av = ρ v wv (1 − ε A ) A The average velocities of the liquid and vapor phases can be obtained by mA 0.02 A = = 0.264m/s wA = ρAε A A 965.35 × 0.01 × 7.85 × 10−3 mv 0.08 v wv = = = 17.27m/s ρv (1 − ε A ) A 0.596 × (1 − 0.01) × 7.85 × 10−3 The relative Reynolds number of the liquid droplets is
v
A
Re = =
ρv wA
A
− wv
v
dA
μv
= 60.70 16.698 × 10−6 The drag coefficient is then 24 24 CD = × (1 + 0.15 × 60.700.687 ) = 1.39 (1 + 0.15Re0.687 ) = Re 60.70 The momentum exchange coefficient can be obtained from eq. (4.31), i.e.,
0.596 × 0.264 − 17.27 × 0.1 × 10−3
ερ 3 K Av = C D A v 4 dA
v
wA
A
− wv
v
3 0.01 × 0.596 = × 1.39 × × 0.264 − 17.27 = 1056.63kg/(m3 -s) −3 4 0.1 × 10 The interactive force between liquid droplets and vapor is
FAv = K Av = 1056.63 × ( 0.264 − 17.27 ) = −1.797 × 104 N/m3 which is negative because the liquid droplet moves slower than the vapor phase.
(w
A
A
− wv
v
)
4.3.3 Energy Equation
The extrinsic phase-average of the energy equation, (3.79), is ∂ ( ρ k hk ) + ∇ ⋅ ρ k Vk hk ∂t
∂pk ′′′ = − ∇ ⋅ q′′ + qk + + Vk ⋅ ∇pk + ∇Vk : k k ∂t in which each term can be obtained by using eqs. (4.13) – (4.15):
(4.35)
248 Transport Phenomena in Multiphase Systems
∂ ρ k hk ∂ ( ρ k hk ) 1 = − ∂t ∂t ΔV ∇ ⋅ ρ k Vk hk = ∇ ⋅ ρ k Vk hk +
∇ ⋅ q′′ = ∇ ⋅ q′′ + k k
³A ρ k hk VI ⋅ n k dAk
k
1 ΔV
k
³A ρ k hk Vk ⋅ n k dAk
k
1 ΔV
k ³A q′′ ⋅ n k dAk
∂ pk ∂pk 1 = − ∂t ∂t ΔV
³A
k
pk VI ⋅ n k dAk
1 § · Vk ⋅ ∇pk Vk ⋅ ∇pk = Vk ⋅ ¨ ∇ pk + ³Ak pk n k dAk ¸ ΔV © ¹ 1 ∇Vk : k ∇Vk : k = ∇ Vk : k + ³ Vk n k dAk : k ΔV Ak where the products of deviations are neglected in the last two equations. Substituting these expressions into eq. (4.35), one obtains ∂ ρ k hk + ∇ ⋅ ρ k Vk hk ∂t ∂ pk ′′′ = −∇ ⋅ q′′ + qk + + Vk ⋅ ∇ pk + ∇ Vk : k k ∂t (4.36) 1 1 − k ³ ρ k hk (Vk − VI ) ⋅ n k dAk − ΔV ³Ak q′′ ⋅ n k dAk ΔV Ak 1ª + − p V ⋅ n dA + Vk ⋅ ³A pk n k dAk + ³A Vk n k dAk : k º » k k ¬³ kI k k ¼ ΔV « Ak where the terms in the square bracket on the right-hand side represent the work per unit volume done by the pressure and shear stress at the interface. These terms reflect conversion of mechanical energy to the thermal energy at the interface, and are usually negligible compared to the other terms. Considering eq. (4.12) and ignoring the product of the deviations, the volume-averaged energy eq. (4.36) becomes ∂ k k k k ε k ρ k hk + ∇ ⋅ ε k ρ k Vk hk ∂t
(
)
(
)
(
)
(
)
′′′ = −∇ ⋅ q′′ + qk + ε k k −
D pk Dt
k
+ ∇ Vk :
k
(4.37)
1 1 k ³Ak ρ k hk (Vk − VI ) ⋅ n k dAk − ΔV ³Ak q′′ ⋅ n k dAk ΔV The fifth term on the right-hand side of eq. (4.37) is interphase enthalpy exchange between all other phases and the kth phase due to phase change, i.e., Π k 1 − ρ k hk (Vk − VI ) ⋅ n k dAk = ¦ m′′′ hk , I (4.38) jk ³Ak ΔV j =1( j ≠ k )
Chapter 4 Generalized Governing Equations: Averaging Formulations 249
where
hk , I
k
is the intrinsic phase-averaged enthalpy of the kth phase at the
interface. The sixth term on the right-hand side of eq. (4.37) is heat transfer from all other phases and the kth phase, Π 1 − (4.39) k jk ³Ak q′′ ⋅ n k dAk = ¦ q′′′ ΔV j =1( j ≠ k ) where
q′′′ jk
is the intensity of heat exchange between phase j and k. It can be
obtained by using Newton’s law of cooling:
q′′′ = jk
hc ΔA j T j
(
j
− Tk
k
ΔV j
)
(4.40)
where hc is the convective heat transfer coefficient, Aj is the area of the interface between phases j and k, and Vj is the volume of the secondary phase in the volume element. Like the momentum exchange coefficient, the interphase heat transfer also depends on the structure of the interfaces. If a secondary phase j is dispersed in the primary phase k – as in the dispersed phase summarized in Table 1.10 – the following empirical correlation recommended is widely used:
§ μk · (4.41) Nu = 2 + (0.4 Re + 0.06 Re ¨ ¨μ ¸ ¸ k ,s ¹ © where the Reynolds number, Re, is obtained by eq. (4.33), the Nusselt number is defined as hc d j Nu = (4.42) kk and all thermal properties of the primary phase are evaluated at Tk except
1/ 2 2/3
1/ 4
) Prk0.4
μk , s , which is evaluated at T j . Equation (4.41) is valid for 3.5 < Re < 7.6 × 104
and 0.71 < Prk < 380 , which covers a wide variety of problems. If the secondary phase is liquid and the primary phase is vapor (gas), eq. (4.41) can be simplified to 1/ Nu = 2 + 0.6 Re1/ 2 Prk 3 (4.43) Substituting eqs. (4.38) and (4.39) into eq. (4.37), the volume-averaged energy equation becomes ∂ k k k k ′′′ ε k ρ k hk + ∇ ⋅ ε k ρ k Vk hk = −∇ ⋅ q′′ + qk k ∂t (4.44) k Π D pk ª q′′′ + m′′′ h k º +ε k + ∇ Vk : k + ¦ jk k ,I « jk » ¼ Dt j =1( j ≠ k ) ¬
(
)
(
)
Example 4.2 Estimate the interphase heat transfer rate between the liquid and vapor phases in Example 4.1.
250 Transport Phenomena in Multiphase Systems
Solution: The thermal conductivity and the Prandtl number of the vapor phase are kv = 0.0248W/m-K and Prv = 0.984 , respectively. Since the primary phase is vapor ( v ) and the secondary phase is liquid ( A ), eq. (4.43) can be used to estimate the heat transfer coefficient, i.e., 1/ Nu = 2 + 0.6 Re1/ 2 Prv 3 = 2 + 0.6 × 60.701/ 2 × 0.9841/ 3 = 6.65 The heat transfer coefficient is then Nukv 6.65 × 0.0248 hc = = = 1649.2W/m 2 -K dA 0.1 × 10−3 Therefore, the interphase heat transfer between the liquid and vapor phases can be obtained by eq.(4.40):
′′′ qAv = =
hc AA TA
(
A
− Tv
v
v
6hc TA
(
VA
A
) = h πd ( T
c 2 A
A
A
− Tv
v
− Tv
dA
) = 6 ×1649.2 × (95 − 100) = −4.95 ×10 W/m
8
π d A3 / 6
)
3
0.1 × 10−3
4.3.4 The Second Law of Thermodynamics
The extrinsic phase-averaged second law of thermodynamics can be obtained by applying volume average to eq. (3.99), i.e., ∂ ( ρ k sk ) § q′′ · q′′′ (4.45) + ∇ ⋅ ( ρ k Vk sk ) + ∇ ⋅ ¨ k ¸ − k ≥ 0 ∂t Tk © Tk ¹
Each term in eq. (4.45) can be evaluated by using eqs. (4.13) – (4.15) , i.e., ∂ ρ k sk ∂ ( ρ k sk ) 1 = − ³ ρ k sk VI ⋅ n k dAk ∂t ∂t ΔV Ak 1 ∇ ⋅ ρ k Vk sk = ∇ ⋅ ρ k Vk sk + ³ ρ k sk Vk ⋅ n k dAk ΔV Ak § q′′ · q′′ q′′ 1 k ∇⋅¨ k ¸ = ∇⋅ k + ³Ak T ⋅ n k dAk Tk ΔV k © Tk ¹
Substituting the above expression into eq. (4.45), one obtains ∂ ρ k sk q′′ q′′′ + ∇ ⋅ ρ k Vk sk + ∇ ⋅ k − k ∂t Tk Tk
1 + ΔV 1 ³Ak ρ k sk (Vk − VI ) ⋅ n k dAk + ΔV q′′ k ³Ak T ⋅ n k dAk > 0 k
(4.46)
Chapter 4 Generalized Governing Equations: Averaging Formulations 251
Considering eq. (4.12) and ignoring the product of the deviations, eq. (4.46) becomes q′′ q′′′ ∂ k k k k ε k ρ k sk + ∇ ⋅ ε k ρ k Vk sk + ∇ ⋅ k − k Tk Tk ∂t (4.47) q′′ 1 1 k + ³ ρk sk (Vk − VI ) ⋅ n k dAk + ΔV ³Ak T ⋅ n k dAk > 0 ΔV Ak k The last two terms on the right-hand side of eq. (4.47) represent interphase entropy exchange between all other phases and the kth phase due to phase change, as well as entropy generation due to heat transfer from all other phases and the kth phase, i.e., Π k 1 ρ k sk (Vk − VI ) ⋅ n k dAk = ¦ m′′′ sk , I (4.48) − jk ³Ak ΔV j =1( j ≠ k ) Π q′′ 1 k − ⋅ n k dAk = ¦ s′′′ (4.49) jk ³ ΔV Ak Tk j =1( j ≠ k )
(
)
(
)
where
sk , I
k
is the intrinsic phase-averaged entropy of the kth phase at the
interface. Substituting eqs. (4.48) and (4.49) into eq. (4.47), the volume-averaged second law of thermodynamics can be expressed as q′′ q′′′ ∂ k k k k ε k ρ k sk + ∇ ⋅ ε k ρ k Vk sk + ∇ ⋅ k − k ∂t Tk Tk (4.50) Π ª m′′′ s k + s′′′ º > 0 −¦ jk » « jk k , I ¼ j =1( j ≠ k ) ¬
(
)
(
)
4.3.5 Species
If the fluid undergoing phase change involves multiple components, it is also necessary to write the equation for conservation of the species mass in the kth phase. The extrinsic phase-average of conservation of species mass can be obtained by ∂ρ k ,i ′′′ (4.51) + ∇ ⋅ ρ k ,i Vk = − ∇ ⋅ J k ,i + mk ,i ∂t where each term can be obtained by ∂ ρ k ,i ∂ρ k ,i 1 = − ³ ρk ,i VI ⋅ n k dAk ∂t ∂t ΔV Ak ∇ ⋅ ρ k ,i Vk = ∇ ⋅ ρ k ,i Vk + 1 ΔV
³A ρ k ,i Vk ⋅ n k dAk
k
252 Transport Phenomena in Multiphase Systems
1 ³ J k ,i ⋅ n k dAk ΔV Ak Substituting the above expression into eq. (4.51), one obtains ∂ ρ k ,i ′′′ + ∇ ⋅ ρ k ,i Vk = −∇ ⋅ J k ,i + mk ,i ∂t (4.52) 1 1 − ³ ρ k ,i (Vk − VI ) ⋅ n k dAk + ΔV ³Ak J k ,i ⋅ n k dAk ΔV Ak where the third and fourth terms in the right-hand side of eq. (4.52) represent mass source (or sink) of the ith component in the kth phase due to phase change from other phases to the kth phase, as well as mass transfer at the interface due to diffusion, respectively, i.e., Π 1 1 − jk ³Ak ρ k ,i (Vk − VI ) ⋅ n k dAk + ΔV ³Ak J k ,i ⋅ n k dAk = ¦ m′′′ ,i (4.53) ΔV j =1( j ≠ k ) where m′′′ ,i represents the mass source (or sink) of the ith component in phase k jk ∇ ⋅ J k ,i = ∇ ⋅ J k ,i + due to phase change from phase j to phase k, as well as diffusive mass transfer at the interface between phases j and k. Substituting eq. (4.53) into eq. (4.52) and considering eq. (4.12), the volume-averaged species equation becomes ∂ ε k ρ k ,i ∂t
(
k
) +∇⋅ ε (
k
ρ k ,i
k
Vk
k
) = −∇ ⋅ J
k ,i
′′′ + mk ,i +
j =1( j ≠ k )
¦
Π
m′′′ ,i jk
(4.54) where the three terms on the right-hand side represent the effects of mass diffusion, species source/sink due to chemical reaction, and phase change at the interfaces. Example 4.3 A cooling tower is a multiphase system. Warm water is sprayed at the top, in the form of water droplets. Gravity drives these droplets downward. The water droplets are cooled by evaporated effects as they flow through the air. Air enters the cooling towers from the bottom on the sides. The system being investigated has a planar geometry. The schematic for this system is presented in Fig. 4.2. The temperature of the water collected at the bottom of the tank is desired, as well as the amount of water lost due to evaporation. Develop the governing equations by the volume-average method, as well as the boundary conditions. Solution: This problem can be solved using a volume-averaged multifluid formulation. The basic assumptions are that the flow field is steady. The planar geometry is long enough that the changes in the zdirected can be neglected; therefore a two-dimensional geometry is utilized. Also, since this is the case, a plane of symmetry exists half the
Chapter 4 Generalized Governing Equations: Averaging Formulations 253
y,
vk
k
x, z, Water in
uk
k
Air/water vapor
wk
k
g Insulated wall
Domain to model
Air in Water
Figure 4.2 Schematic of the cooling of droplets.
Air in
distance between the two adiabatic walls. The pressure effects on energy are considered negligible. The final assumption is that intrinsic phase average of the products of deviation is zero, k k k ˆˆ ˆˆ ˆˆ VV = Vh = Vω =0.
k k kk k k
The volume-averaged continuity equation, eq. (4.22), for two phases is: ∂ ∂ A A A A ′′′ ε A ρ A uA + ε A ρA vA = mvA (4.55) ∂x ∂y ∂ ∂ v v v v ′′′ ′′′ ε v ρ v uv + ε v ρ v vv = mAv = − mvA (4.56) ∂x ∂y The volume-averaged momentum equations for both phases in the x- and y- directions are obtained from eq. (4.34): Liquid x-momentum A A A A ∂ uA A A ∂ uA + ε A ρA vA = ε A ρ A uA ∂x ∂y
( (
) )
( (
) )
−
∂ε A pA ∂x
A
+
∂ uA ∂§ ¨ ε A μA ∂x ¨ ∂x ©
A
· ∂§ ∂ uA ¸ + ¨ ε A μA ¸ ∂y ¨ ∂y ¹ ©
A
· ′′′ ¸ + FvA , x + mvA uA , I ¸ ¹ (4.57)
A
254 Transport Phenomena in Multiphase Systems
Liquid y-momentum
ε A ρA
=−
A
uA
A
∂ vA ∂x
A
+ ε A ρA
A
vA
A
∂ vA ∂y
A
∂ε A pA ∂y
A
+
∂ vA ∂§ ¨ ε A μA ∂x ¨ ∂x ©
A
A
· ∂§ ∂ vA ¸ + ¨ ε A μA ¸ ∂y ¨ ∂y ¹ ©
A
A
· ¸ ¸ ¹
(4.58)
′′′ + FvA, y + mvA vA , I Vapor x-momentum v ∂ uv ε v ρ v v uv ∂x ∂ uv ∂§ + ¨ ε v μv ¨ ∂x ∂x © Vapor y-momentum
v v
− ε A ρA
v
g ∂ uv ∂y
v v
+ ε v ρv
vv
=−
∂ε v pv ∂x
v
· ∂§ ∂ uv ¸ + ¨ ε v μv ¸ ∂y ¨ ∂y ¹ ©
· v ¸ + FAv , x + mAv uv , I ′′′ ¸ ¹ (4.59)
v
ε v ρv
=−
v
uv
v
∂ vv ∂x
v
+ ε v ρv
v
vv
v
∂ vv ∂y
∂ε v pv ∂y
v
+
∂ vv ∂§ ¨ ε v μv ∂x ¨ ∂x ©
v
v
· ∂§ ∂ vv ¸ + ¨ ε v μv ¸ ∂y ¨ ∂y ¹ ©
v
v
· ¸ ¸ ¹
(4.60)
′′′ + FAv , y + mAv vv , I
− ε v ρv
g
The volume-averaged energy equations for the liquid and vapor phases are obtained from eq. (4.44): Liquid energy ∂ hA ∂x
A A
ε A ρA
A
uA
+ ε A ρA
A
vA
A
∂ hA ∂y
A
= (4.61)
∂ TA ∂§ ¨ ε A kA ¨ ∂x ∂x © Vapor energy v∂ ε v ρ v uv
· ∂§ ∂ TA ¸ + ¨ ε A kA ¸ ∂y ¨ ∂y ¹ © hv
v
A
· ′′′ ′′′ ¸ + qvA + mvA hAv ¸ ¹ ∂ hv ∂y
v
∂x
+ ε v ρv
v
vv
v
= (4.62)
v ∂ Tv · ∂ § ∂ Tv ∂§ ¨ ε v kv ¸ + ¨ ε v kv ∂x ¨ ∂x ¸ ∂y ¨ ∂y © ¹ © where the enthalpy is defined as: A hA = ³ c pA dT
v
· ¸ + qAv ′′′ ¸ ¹
(4.63)
Chapter 4 Generalized Governing Equations: Averaging Formulations 255
hv
v
= ³ c pv dT
N
(4.64)
There are multiple components in the vapor phase; therefore, the specific heat is a function of the mass fraction of those components.
c pv = ¦ ωi c pi
i =1
(4.65)
Air is considered to be one component, and water vapor is considered to be the other component. Therefore, only one species equation is needed. The species equation for water vapor, eq. (4.54), is: v v v v ∂ ωv v v ∂ ωv ε v ρ v uv + ε v ρv vv = ∂x ∂y (4.66) v v ∂ ωv · ∂ § ∂ ωv · ∂§ ¨ ε v Dv , air ¸ + ¨ ε v Dv , air ¸ + mAv ′′′ ∂x ¨ ∂x ¸ ∂y ¨ ∂y ¸ © ¹ © ¹ It is also known that εA + εv = 1 (4.67)
′′′ mvA = hm AAv (ωsat − ωv ,i ) ΔV (1 − ωsat )
(4.68)
−
ωsat =
p0 Rg Tv
v
′′′ qvA =
§h§1 exp ¨ − Av ¨ v ¨ Rg ¨ Tv ρv © © hc AAv (Tv − TA )
ΔV
v
1 ·· ¸¸ T0 ¸ ¸ ¹¹
(4.69)
(4.70)
A
FvA = − FAv = K vA Vv
(
v
− VA
)
(4.71) (4.72)
ρv
v
=
P Rg Tv
v
The surface to volume ratio can be calculated by assuming that the liquid is in the form of spherical droplets. Therefore, the surface to volume ratio is defined in terms of the liquid droplet diameter, d A , and the volume fraction of the liquid. AAv d = εA A (4.73) ΔV 6 There are 15 equations and 16 unknowns (9 equations and 6 constitutive relations). To solve this problem, it is assumed that the liquid is in droplets, and there is no pressure gradient within each droplet. Therefore the liquid pressure gradient is assumed to be equal to the vapor pressure gradient.
256 Transport Phenomena in Multiphase Systems
, (4.74) = = ∂x ∂x ∂y ∂y Now the liquid pressure is eliminated, therefore we have 13 equations and 13 unknowns, which is well posed. The boundary conditions are:
At walls v uv = vv ∂ TA ∂x At air inlet
A
∂ pA
A
∂ pv
v
∂ pA
A
∂ pv
v
v
= uA
v
A
= vA
v
A
=0 =0
(4.75) (4.76)
=
∂ Tv ∂x
v
=
∂ ωv ∂x
v
pv
= ρv
v
gh
(4.77) (4.78) (4.79) (4.80) (4.81) (4.82)
k
Tv
= Tin
εA = 0
At top TA pv
A
= Ttop
v
ε A = ε A ,in
=0
, ωv
v
At symmetry ∂φ k k = 0 ,φ = vk , Tk , ε A , pk ∂x k u k = 0 , k = A, v At bottom surface v ωv = ωv , sat −
(4.83) (4.84) (4.85)
ε v ρv
v
Dv ,i , air ∂ ωv
v v
v
1 − ωv
−ε v kv
∂y
= ε v ρv
v
v
vv
v
(4.86)
∂ Tv ∂y
uv
= ε v ρv
v
vv
v
hAv
(4.87) (4.88) (4.89)
= uA
v
A
≈0
A
Tv
= TA
Chapter 4 Generalized Governing Equations: Averaging Formulations 257
4.4 Volume-Averaged Homogeneous Model
The multifluid model presented above is obtained by performing phase averaging as defined in eqs. (4.6) and (4.7). If spatial averaging is performed for all phases within a multiphase control volume, the homogeneous (or mixture) model can be obtained. The relationship between volume averaging and phase averaging is given in eq. (4.9), which indicates that the homogeneous model can be obtained by summing the individual phase equations of the multifluid model.
4.4.1 Continuity Equation
The continuity equation for phase k in the multifluid model is expressed by eq. (4.22). Summing the continuity equations for all Π phases together, one obtains Π Π Π ∂§Π k· k k (4.90) ¦ ε k ρ k ¸ + ∇ ⋅ ¦ ε k ρ k Vk = ¦ ¦ m′′′ jk ∂t ¨ k =1 k =1 k =1 j =1( j ≠ k ) © ¹ The right-hand side of equation (4.90) must be zero because the total mass of all phases produced by phase change must equal the total mass of all phases consumed by phase change. Considering this fact and eq. (4.20), the continuity equation becomes Π ∂ρ k k + ∇ ⋅ ¦ ε k ρ k Vk = 0 (4.91) ∂t k =1 The bulk velocity of the multiphase mixture is the mass-averaged velocity of all the individual phases: 1Π k k V= (4.92) ¦ ε k ρ k Vk
ρ
k =1
Substituting eq. (4.92) into eq. (4.91), the final form of the continuity equation for a multiphase mixture is ∂ρ (4.93) +∇⋅ ρ V = 0 ∂t It can be seen that eq. (4.93) has the same form as the local continuity equation (3.38), except that the volume-averaged density and velocity are used in eq. (4.93), where
ρ = ¦ ε k ρk .
k k =1
Π
4.4.2 Momentum Equation
The momentum equation for phase k in the multifluid model is expressed in eq. (4.34). By adding together the momentum equations for all Π phases, one obtains
258 Transport Phenomena in Multiphase Systems
∂§Π ¦ ε k ρk ∂t ¨ k =1 © = ∇ ⋅ ¦εk
k =1
k
Vk
k
k
Π
′ k
+ ¦ ε k ρk
k =1 k
Π
(
· §Π ¸ + ∇ ⋅ ¨ ¦ ε k ρk ¹ © k =1
k
k
Vk Vk
k
)X + ¦
Π
k =1 j =1( j ≠ k )
¦
Π
(F
· ¸ ¹
+ m′′′ jk Vk , I
k
jk
)
(4.94)
The stress tensor of the multiphase mixture is
′ = ¦εk
k =1
Π
k
2 = − p I + μ ª∇V + ∇V T º − μ (∇ ⋅ V )I ¬ ¼3
Π Π
(4.95)
The summation of all interphase forces must be zero since
k =1 j =1( j ≠ k )
F jk = − Fkj , i.e.,
¦¦
F jk = 0
(4.96)
Considering eqs. (4.92), (4.95) and (4.96), the momentum equation becomes ∂ k k· §Π ρ V + ∇ ⋅ ¨ ¦ ε k ρ k Vk Vk ¸ = ∇ ⋅ ′ + ρ X + M ′′′ (4.97) I ∂t © k =1 ¹ where
(
)
M ′′′ = ¦ I
Π
k =1 j =1( j ≠ k )
¦
Π
m′′′ jk
Vk , I
k
(4.98)
Equation (4.98) represents the momentum production rate due to interaction between different phases along their separating interfaces. It must be specified according to the combination of phases in the multiphase system that is under consideration.
4.4.3 Energy Equation
By summing the energy equations for all Π phases in the multifluid model, eq. (4.44), one obtains ∂§Π k k· k k· §Π §Π ·Π ′′′ k ¨ ¦ ε k ρ k hk ¸ + ∇ ⋅ ¨ ¦ ε k ρ k Vk hk ¸ = −∇ ⋅ ¨ ¦ q′′ ¸ + ¦ qk ∂t © k =1 ¹ © k =1 ¹ © k =1 ¹ k =1
+¦ εk
k =1
Π
D pk Dt
k
+ ¦ ∇ Vk
k =1
Π
k
:εk
k k
+¦
Π
k =1 j =1( j ≠ k )
¦
Π
q′′′ + ¦ jk
Π
k =1 j =1( j ≠ k )
¦
Π
m′′′ hk , I jk
k
(4.99) The mass average enthalpy of the multiphase mixture is 1Π k k h= ¦ ε k ρ k hk (4.100)
ρ
k =1
The fifth term on the right-hand side of eq. (4.99) is for summation of all interphase heat transfer and it must be zero. The last term on the right-hand side of eq. (4.99) accounts for contribution of interphase phase change energy flux due to phase change; it can be defined as
Chapter 4 Generalized Governing Equations: Averaging Formulations 259
k =1 j =1( j ≠ k )
¦¦
Π
Π
m′′′ hk , I jk
Π
k
′′′ = qI m′′′ = 0 .
(4.101)
It is usually not zero although ¦
k =1 j =1( j ≠ k )
¦
Π
Considering eqs. (4.100) and (4.101), the energy equation (4.99) becomes Dp ∂ k k· §Π ′′′ + q′′′ + ∇V : + qI ρ h + ∇ ⋅ ¨ ¦ ε k ρ k Vk hk ¸ = −∇ ⋅ q′′ + ∂t Dt © k =1 ¹ (4.102)
(
)
4.4.4 The Second Law of Thermodynamics
Summing the equations of the second law of thermodynamics expressed by eq. (4.50), one obtains ∂§Π k k· k k· §Π ¦ ε k ρ k sk ¸ + ∇ ⋅ ¨ ¦ ε k ρk Vk sk ¸ ¨ ∂t © k =1 ¹ © k =1 ¹ (4.103) Π q′′ Π Π Π ′′′ qk ª m′′′ s j − s k + s′′′ º > 0 +∇ ⋅ ¦ k − ¦ −¦ ¦ j k jk » « jk ¼ k =1 Tk k =1 Tk k =1 j =1( j ≠ k ) ¬ where Π Π j k (4.104) ¦ ¦ ª m′′′ s j − sk + s′′′ º = 0 jk » « jk ¼ k =1 j =1( j ≠ k ) ¬ The mass-averaged entropy of the multiphase mixture is 1n k k (4.105) s= ¦ ε k ρk sk
(
)
(
)
ρ
k =1
Substituting eqs. (4.104) and (4.105) into eq. (4.103), the second law of thermodynamics for a multiphase mixture becomes Π Π q′′ Π q′′′ ∂ ( ρ s ) + ∇ ⋅ § ¦ ε k ρk k Vk sk k · + ∇ ⋅ ¦ Tk − ¦ Tk > 0 (4.106) ¨ ¸ ∂t k =1 k =1 © k =1 ¹ k k
4.4.5 Species
Summing the equations for conservation of species mass, eq. (4.54), for all phases yields k· k ∂§Π k· §Π ¦ ε k ρ k ,i ¸ + ∇ ⋅ ¨ ¦ ε k ρ k ,i Vk ¸ ¨ ∂t © k =1 ¹ © k =1 ¹ (4.107) Π Π Π Π § · ′′′ = −∇ ⋅ ¨ ¦ J k ,i ¸ + ¦ mk ,i + ¦ ¦ m′′′ ,i jk k =1 j =1( j ≠ k ) © k =1 ¹ k =1
260 Transport Phenomena in Multiphase Systems
By applying eq. (4.20) to the mass density of the ith component, one obtains
ρ i = ¦ ε k ρ k ,i
k =1
Π
k
(4.108)
In accordance with the conservation of mass, the mass source (or sink) of the ith component due to phase change in all phases must add up to zero, i.e.,
k =1 j =1( j ≠ k )
¦¦
Π
Π
m′′′ ,i = 0 jk
(4.109)
Substituting eqs. (4.108) and (4.109) into eq. (4.107), and using the massaveraged velocity defined in eq. (4.92), the conservation of species mass becomes ∂ ρi + ∇ ⋅ ρi V = −∇ ⋅ J i + mi′′′ (4.110) ∂t
Example 4.4 Solve problem in Example 4.3 using the homogeneous model using the same assumptions made in Example 4.3. The water droplets are small; therefore, make further reasonable assumptions to find a relationship between the liquid velocity and vapor velocity so that the homogeneous model can be applied. The flow variable to be solved is the mass averaged velocities, u and v , as well as the average temperature, T . Solution: Use the same assumptions listed in Example 4.3. To make this problem a homogenous problem, more assumptions need to be made. The inertial terms in the liquid momentum equation can be assumed to be negligible because the droplets are very small. It can also be assumed that the evaporation rates and pressure gradients do not affect liquid momentum equation. Also, assume that the temperature of the liquid v A and vapor are in local equilibrium, TA = Tv . Based on the assumptions above, the liquid momentum equation in the x-direction reduces to a statement that the velocities of the liquid and vapor are equivalent.
u A = uv = u (4.111) Similarly, the liquid momentum equation in the y-direction reduces to a statement that drag force between the liquid and vapor is equivalent to the gravitational force on the liquid droplets. vA
A A
v
= vv
v
−
ε A ρA
K vA
A
g
(4.112)
Chapter 4 Generalized Governing Equations: Averaging Formulations 261
Therefore, the y component of the mass averaged velocity is
(ε v=
v
ρv
v
+ ε A ρA
A
(ε
)v
v
v
v
− ε A2 ρA
A
(
A2
v
ρv
+ ε A ρA
)
)
gK v−A1
= vv
v
−
ε A2 ρA
(
A2
ρ K vA
)
g
(4.113) Therefore, the continuity equations, eq. (4.93), for the liquid phase and the mixture are: Liquid A2 A ·· ε2 ρ g § ε A ρA ∂ ∂§ ¨ ε A ρA v + A A ′′′ ¨ − 1¸ ¸ = mvA (4.114) ( ε A ρA u ) + ¨ρ ¸¸ ∂x ∂y ¨ K vA © ¹¹ © Mixture ∂ ∂ (4.115) ( ρ u ) + ∂y ( ρ v ) = 0 ∂x The homogeneous momentum equations, eq. (4.97), are: A2 § · ε A2 ρA g ¸ ∂u v ∂u v¨ ε v ρv u + ε v ρv ¨ v + ρ K vA ¸ ∂y ∂x ¨ ¸ (4.116) © ¹ ∂p ∂§ ∂u · ∂ § ∂u · =− + ¨ ε v μv ¸ ¸ + ¨ ε v μv ∂x ∂x © ∂x ¹ ∂y © ∂y ¹
(
)
§ ε A2 ∂¨ v ε v ρv u ¨ v + ∂x ¨ ©
(ρ )
A
A2
K vA
· g¸ ¸ + ε v ρv ¸ ¹
A2 § · ε A2 ρA g¸ v¨ ¨v + ρ K vA ¸ ¨ ¸ © ¹
(
)
§ A2 A2 § · § ·· ε A2 ρA ε l2 ρA g¸ g ¸¸ ¨ ∂p ∂¨ ∂ ∂¨ ⋅ ¨v + =− + ¨ ε v μv ¨ v + ¸ ρ K vA ¸ ρ K vA ¸ ¸ ∂y ¨ ∂y ∂x ¨ ∂x ¨ ¸ ¸ © ¹ © ¹¹ © § A2 § ·· ε l2 ρA g ¸¸ ∂¨ ∂¨ + ¨ ε v μv ¨ v + ¸−ε ρ g ρ K vA ¸ ¸ v v ∂y ¨ ∂y ¨ ¸ © ¹¹ ©
(
)
(
)
(
)
(4.117) The homogenous energy equation, eq. (4.102), is: A2 § · v A ε A2 ρA g ¸∂ h § · v ∂ hv v¨ A ∂ hA v ¸ + ε v ρv ¨ v + + ε A ρA u ¨ ε v ρv ¸ ∂y ¨ ∂x ∂x ¸ ρ K vA ¸ ¨ © ¹ © ¹
(
)
v
262 Transport Phenomena in Multiphase Systems
+ε A ρA +
A
A § ε A ρA g § ε A ρA ¨v + ¨ ¨ρ ¨ K vA © ©
A
· · ∂h ∂T · ∂§ − 1 ¸ ¸ A = ¨ ( ε v kv + ε A k A ) ¸ ¸ ¸ ∂y ∂x © ∂x ¹ ¹¹
∂T · ∂§ ′′′ ¨ ( ε v kv + ε A k A ) ¸ + mvA hAv ∂y © ∂y ¹ (4.118)
A
The enthalpies of the liquid and the vapor are: hA = ³ c pA dT (4.119) (4.120)
v
hv
= ³ c pv dT
N
There are multiple components in the vapor phase; therefore the specific heat is a function of the mass fraction of those components. c pv = ¦ ωi c pi
i =1
(4.121)
The homogenous species equation, eq. (4.110), is: A2 § · v ε A2 ρA g ¸∂ ω v ∂ ωv v¨ v + ε v ρv ¨ v + ε v ρv u ∂x ρ K vA ¸ ∂y ¨ ¸ © ¹
(
)
v
(4.122)
= Also
∂ ωv ∂§ ¨ ε v Dv ,air ∂x ¨ ∂x ©
v
· ∂§ ∂ ωv ¸ + ¨ ε v Dv ,air ¸ ∂y ¨ ∂y ¹ ©
v
· ¸ + mAv ′′′ ¸ ¹ (4.123) (4.124)
εA + εv = 1
′′′ mvA = hm AAv ωsat − ωv ΔV (1 − ωsat )
v
(
v
)
v
ωsat =
p0 Rg Tv
ρv
v
v
ρv
=
§h§1 exp ¨ − Av ¨ ¨ Rg ¨ Tv © © P
v
−
1 ·· ¸¸ T0 ¸ ¸ ¹¹
(4.125) (4.126)
Rg Tv
The surface to volume ratio can be calculated by assuming that the liquid is in the form of spherical droplets. Therefore, the surface to volume ratio is defined in terms of the liquid droplet diameter, d A , and the volume fraction of the liquid. AAv d = εA A (4.127) ΔV 6
Chapter 4 Generalized Governing Equations: Averaging Formulations 263
The problem was reduced to 6 equations with 3 constitutive relations from 9 equations. In total, the present formulation has 9 equations and 9 unknowns and is well posed. At walls u=v=0 (4.128) ∂T ∂x At top = ∂ ωv ∂x
v
=0
(4.129) (4.130) (4.131) (4.132)
T = Ttop
p =0
ε A = ε A ,in
At air inlet p = ρv
v
gh
(4.133) (4.134) (4.135)
v
T = Tin
εA = 0
At symmetry ∂φ = 0,φ = v , T , ε A , ωv ∂x u=0 At bottom surface v ωv = ωv , sat
(4.136)
(4.137) (4.138)
∂ ωv
v v
ρ v = − ε v ρv
(
v
+ ε A ρA
A
)1− ω
Dv ,air
v
∂y
A2
− ε A ρA
A
ε A ρA
K vA
A
g
(4.139)
§ · ε A2 ρA g¸ ¨ = ε v ρv ¨ v + hAv (4.140) ( ε v kv + ε A kA ) ∂y ρ K vA ¸ ¨ ¸ © ¹ From the above problem formulation, the homogenous model can still be used for counter-flowing phases. However, more assumptions were made to solve the problem as a homogeneous model. So the trade off for reducing the number of equations is the accuracy of the solution. ∂T
v
(
)
4.5 Area-Averaged Models for Channel Flows
The primary application of area averaging is in channel flows. Area-averaged governing equations are obtained by performing area averaging of the governing equations over the entire cross-section of the channel; this is distinct from
264 Transport Phenomena in Multiphase Systems
volume averaging, wherein the average is performed over a small-volume ΔV. The resulting area-averaged governing equations will be one-dimensional in nature, rather than three-dimensional like volume-averaged equations. Therefore, information about the distribution of properties in the cross-section is lost in the area-averaged governing equations. The jump conditions at the interfaces between two phases, as well as the exchange of momentum and energy between the channel wall and the fluid, must be expressed in the form of empirical correlations. For a multiphase system that includes Π different phases flowing in a channel, the total cross-sectional area of the channel equals the summation of the areas occupied by individual phase, i.e.,
Ac = ¦ Ak
k =1 th
Π
(4.141)
The volume fraction of the k phase, ε k , for channel flow is defined as the ratio of the cross-sectional area of the kth phase to the total area, i.e., A εk = k (4.142) Ac It follows from eq. (4.141) that
k =1
¦εk
Π
=1
(4.143)
Eulerian area averaging is performed over the cross-section of a channel and the governing equations are reduced to one-dimensional 1Π Φ= (4.144) ¦ ³ Φ k dA Ac k =1 Ak If the average is performed for the kth phase only, the phase averages are obtained Intrinsic phase-average
Φk
k
=
1 Ak
³A Φ k dA
k
(4.145)
Extrinsic phase-average
1 (4.146) ³ Φ k dA Ac Ak The area averaging and phase averaging defined in eq. (4.144) and eqs. (4.145) and (4.146) are related by Φk =
Φ = ¦ Φk = ¦ ε k Φk
k =1 k =1
Π
Π
k
(4.147)
The area-averaged homogeneous model treats the flow as a homogeneous mixture with all phases flowing at the same area-averaged velocity. Since averaging is not performed in the axial direction, the governing equations are still local instantaneous in the axial direction. The basic equations for onedimensional homogeneous equilibrium flow in a channel will be introduced here.
Chapter 4 Generalized Governing Equations: Averaging Formulations 265
4.5.1 Homogeneous Model
Continuity Equation
The area-averaged continuity equation is ∂ρ ∂ (4.148) + ( ρ w )=0 ∂t ∂z where ρ and w are, respectively, the area-averaged density and axial velocities of the multiphase mixture. The product of the deviations of density and velocity is neglected. For a mixture of liquid and vapor, it can be expressed as v A ρ = α ρv + (1 − α ) ρA (4.149) where α is the vapor or holdup void fraction, α = ε v , which represents the time-averaged volumetric fraction of vapor in a two-phase mixture.
Momentum Equation
The area-averaged momentum equation is §∂ w ∂w· ∂p P ρ¨ (4.150) +w − ρ g cos θ − τ w ¸=− ∂t ∂z ¹ ∂z Ac © where gravity is the only body force considered, τw is the average shear stress at the wall of the channel, P is the perimeter of the duct, θ is the inclination of the axial direction of the duct to the vertical, and Ac is the cross-sectional area of the channel.
Energy Equation
The area-averaged energy equation is 2 2 § § w ·º ∂ ª w ·º ∂ª «ρ¨e+ ¸» + « ρ w ¨ h + ¸» ¨ ¨ ∂t « 2 ¸ » ∂z « 2 ¸» © ¹¼ © ¹¼ ¬ ¬ P ′′ qw − ρ w g cos θ = Ac where
e
(4.151)
h are the area-averaged internal energy and enthalpy of the ′′ multiphase flow, and qw is the heat flux at the wall of the channel. Equation (4.151) is the total energy (including internal energy, kinetic energy, and potential energy) balance for the channel. The mechanical energy balance can be obtained by multiplying eq. (4.150) by w , i.e.,
and
266 Transport Phenomena in Multiphase Systems
2 ª∂ § w 2· ∂ § w ·º ¸+ w ¨ ¸» ρ« ¨ ¨ ¸ ∂z ¨ 2 ¸ » « ∂t © 2 ¹ © ¹¼ ¬ (4.152) ∂p P =− w − ρ g w cos θ − τ w w ∂z Ac Subtracting eq. (4.152) from eq. (4.151) and replacing internal energy with enthalpy ( h = e + p / ρ ), one obtains the final form of the energy equation:
∂p· ′′ Pqw P τw 1 §∂ p (4.153) +w + ¨ ¸+ w ∂t ∂z ∂z ¹ ρ © ∂t Ac ρ Ac ρ where the three terms on the right-hand side represent the effects of compressibility, viscous dissipation, and heat transfer from the wall. +w =
The Second Law of Thermodynamics
∂h
∂h
The second law of thermodynamics for multiphase flow in a channel is ′′ Pqw · Pτ ∂ ∂ 1§ (4.154) ( ρ s ) + ∂z ( ρ w s ) − T ¨ w A ρw + A ρ ¸ ≥ 0 ¨ ¸ ∂t c c © ¹
Species
The area-averaged conservation of species mass is ∂ ρi ∂ ρi ∂ Ji +w =− + mi′′′ (4.155) ∂t ∂z ∂z where the two terms on the right-hand side represent the effect of species mass flux and species mass source/sink. The above governing equations are valid for a channel with constant crosssectional area. The following example shows the derivation of the governing equations for a channel with variable cross-sectional area.
Example 4.5 Derive the continuity, momentum, and energy equations for the channel with an inclination angle of as shown in Fig. 4.3, using the area-averaged homogeneous model by making mass, momentum and energy balances on a channel element. The cross-sectional area Ac is changing along the length. Solution: The governing equation for a channel with variable crosssectional area can be obtained by applying various conservation laws to the channel element dz (control volume). For a channel element without mass generation, the mass balance in the channel can be written as
Chapter 4 Generalized Governing Equations: Averaging Formulations 267
z+dz
z
Figure 4.3 Homogeneous model for a channel with variable cross-section.
ª mass º ª mass º ª mass º «inflow » − «outflow » = «storage » (4.156) « »« »« » « rate » « rate » « rate » ¬ ¼¬ ¼¬ ¼ The mass inflow rate is the mass flow rate at z, and it is ρ w Ac . The mass outflow rate is the mass flow rate at z+dz, written as ∂ρ ∂ ρ w Ac + ( ρ w Ac ) dz . The mass storage rate is Ac dz . ∂t ∂z Therefore, eq. (4.156) becomes ∂ρ ∂ ρ w Ac − ρ w Ac − ( ρ w Ac ) dz = Ac dz (4.157) ∂z ∂t which can be rearranged to obtain the continuity equation: ∂ρ ∂ + ( ρ w Ac ) = 0 Ac (4.158) ∂t ∂z The momentum balance for the channel element dz can be written as ª momentum º ª momentum º ª momentum º ªsum of forces acting º «storage rate » + « outflow rate » − «inflow rate » = « on control volume » ¬ ¼¬ ¼¬ ¼¬ ¼ (4.159)
268 Transport Phenomena in Multiphase Systems
The momentum storage rate is ∂ ( ρ Ac dz w ) / ∂t . The momentum inflow rate is the momentum flow rate at z, and it is ρ w Ac w . The momentum outflow rate is the momentum flow rate at z+dz, and it can be ∂ obtained by ρ w Ac w + ( ρ w Ac w ) dz . The forces acting ∂z on the control volume include the pressure acting on surfaces z and z+dz, gravity, and the shear force acting on the wall of the channel, i.e., ªsum of forces acting º «on control volume » ¬ ¼
ª ∂ ( p Ac ) º = p Ac − « p Ac + dz » − ρ g cosθ Ac dz − τ w dzP ∂z « » ¬ ¼ Therefore, eq. (4.159) becomes ∂ ( ρ Ac dz w ) ª ∂ º + « ρ w Ac w + ( ρ w Ac w ) dz » − ρ w Ac w ∂t ∂z ¬ ¼ ª º ∂ ( p Ac ) dz » − ρ g cos θ Ac dz − τ w dzP (4.160) = p Ac − « p Ac + ∂z « » ¬ ¼ The momentum equation is obtained by rearranging eq. (4.160), i.e., ∂( ρ w ) ∂ ∂ ( p Ac ) Ac + ( ρ w Ac w ) = − − ρ g cos θ Ac − τ w P ∂t ∂z ∂z (4.161) Rewriting the left-hand side of eq. (4.161) yields ª ∂ρ º ª∂ w ∂wº ∂ w « Ac + ( ρ Ac w ) » + Ac ρ « +w » ∂t ∂z ∂z ¼ ¬ ¼ ¬ ∂t (4.162) ∂ ( p Ac ) =− − ρ g cos θ Ac − τ w P ∂z According to eq. (4.158),the first bracket at the left-hand side of eq. (4.162) is zero. Therefore, the momentum equation is ª∂ w ∂wº 1 ∂ ( p Ac ) P +w − ρ g cos θ − τ w ρ« »=− ∂t ∂z ¼ ∂z Ac Ac ¬ (4.163) The energy equation can be obtained by applying the conservation of energy to the control volume, i.e., ª energy º ª energy º ªenergy º «inflow » − « outflow » = «storage » (4.164) « »« »« » « rate » « rate » « rate » ¬ ¼¬ ¼¬ ¼
Chapter 4 Generalized Governing Equations: Averaging Formulations 269
The energy inflow rate includes the energy flow into the control volume at surface z,
ρ Ac w
′′ transferred into the control volume through the wall, qw Pdz . The energy outflow rate is the energy flow out from the control volume at surface z+dz
+
(h+w
)
2
/ 2 + gz cos θ , and the heat
)
and
it
is
2
ρ Ac w
∂ª 2 º « ρ Ac dz e + w / 2 » . Thus, the energy balance for the channel ¬ ¼ ∂t element becomes 2 § · w " ρ Ac w ¨ h + + gz cos θ ¸ + qw Pdz ¨ ¸ 2 © ¹ 2 § · w ° − ® ρ Ac w ¨ h + + gz cos θ ¸ ¨ ¸ 2 ° © ¹ ¯ 2 2 ½ § ·º ° ∂ ª § w w ·º ∂ª ¸» + « ρ Ac w ¨ h + + gz cos θ ¸ » dz ¾ = « ρ Ac dz ¨ e + ¨ ¸ » ° ∂t « ¨ ∂z « 2 2 ¸» © ¹¼ ¿ © ¹¼ ¬ ¬ (4.165) Rearranging eq. (4.165) yields the energy equation: 2 2 § § w ·º 1 ∂ ª w ·º ∂ª «ρ¨e+ « ρ Ac w ¨ h + ¸» + ¸» ¨ ¨ ∂t « 2 ¸ » Ac ∂z « 2 ¸» © ¹¼ © ¹ ¼ (4.166) ¬ ¬ P" = qw − ρ w g cos θ Ac When the channel’s cross-sectional area is constant, eqs. (4.158), (4.163), and (4.166) reduce to eqs. (4.148), (4.150), and (4.151).
∂ª ρ Ac w ¬ ∂z «
(h+w
(h+ w
2
/ 2 + gz cos θ
)
(
)
/ 2 + gz cos θ º dz. The energy storage rate is » ¼
4.5.2 Separated Flow Model
In the area-averaged homogeneous model, all flow parameters – both velocities and fluid properties – are assumed to be the same across all phases. A slightly more complex approach, the separated-flow model (Hewitt, 1998) is the intermediate between the homogeneous and multifluid models. It accommodates selected parameters that differ between phases, but maintains the single set of conservation equations that is characteristic of the homogeneous model. Recall that the two-fluid model, the simplest case of the multifluid model, recognizes that the two phases can have different properties and velocities.
270 Transport Phenomena in Multiphase Systems
Eulerian phase averaging is performed for each individual phase, yielding separate equations of continuity, momentum, energy, entropy, and species mass for each phase. These equations are solved simultaneously and in conjunction with rate equations that describe how the phases interact with each other and the system boundaries. A simple version of the two-fluid model is the separated-flow model, in which only one parameter, such as velocity, is allowed to differ between the two phases. That parameter is obtained for each phase via phase-averaging, as defined in eq. (4.145). A single set of conservation equations is written for the combined flow, but it uses both velocities. Therefore, the separated flow model is neither multifluid nor homogeneous, but has features of both. The multifluid feature is the recognition of different velocities between phases. The homogeneous feature is the definition of only one governing equation for each of the quantities – mass, momentum, energy, entropy, and species mass – which applies to all phases and across the entire flow cross-section. In addition, the pressure across any given cross-section of a channel carrying a multiphase flow is assumed to be the same for both phases. The following example shows the separated flow model derivation of the governing equations for steady flow in a channel with variable cross-sectional area by using mass, momentum, and energy balance for an element.
Example 4.6. Derive the steady-state continuity, momentum and energy equations for liquid-gas two-phase flow in the channel shown in Fig. 4.4 by using the area-averaged separated-flow model as discussed above. The channel has a cross-sectional area of Ac(z) and an inclination angle of θ. Solution: For a channel element without mass generation, the mass balance in the channel satisfies eq. (4.156). The mass inflow rate is the mass flow rate at z, and it is v v A A ρA wA AA + ρv Av Av . The mass outflow rate is the mass flow
rate at z+dz, and it can be obtained by +
ρA
A
wA
A
AA + ρv
v
wv
v
Av
d v v A A ρA wA AA + ρv wv Av dz . The mass storage rate for dz steady-state flow is zero. Therefore, eq. (4.156) becomes v v v v A A A A ρA wA AA + ρv wv Av − ª ρA wA AA + ρv wv Av ¬ d A A v v º (4.167) + ρA wA AA + ρv wv Av dz » = 0 dz ¼ The above equation can be rearranged to obtain the continuity equation: d v v A A ρA wA AA + ρv wv Av = 0 (4.168) dz
(
)
(
)
(
)
Chapter 4 Generalized Governing Equations: Averaging Formulations 271
z+dz
z
Figure 4.4 Separated flow model for a channel with variable cross-section.
The momentum balance for the channel element dz satisfies eq. (4.159), except that the momentum storage rate is zero for a steady-state flow. The momentum inflow rate is the momentum flow rate at z, and it is v v v A A A ρA wA AA wA + ρv wv Av wv . The momentum outflow rate is the momentum flow rate at z+dz, and it can be obtained by
ρA
A
wA
v
A
AA wA
v
A
+ ρv
v
v
wv
v
Av wv
v
+
d dz
+ ρv
wv
Av wv
) dz . The forces acting on the control volume
(ρ
A A
wA
A
AA wA
A
include pressures acting on surfaces z and z+dz, gravity, and the shear force on the wall of the channel, i.e., ªsum of forces acting º « on control volume » = p Ac ¬ ¼
ª º d p Ac v A dz » − [α ρA + (1 − α ) ρv ]g cos θ Ac dz − τ w dzP − « p Ac + dz « » ¬ ¼ Therefore, eq. (4.159) becomes d A A A v v v ρA wA AA wA + ρv wv Av wv dz = p Ac dz ª d ( p Ac ) º v A dz » − ªα ρA + (1 − α ) ρv º g cos θ Ac dz − τ w dzP − « p Ac + ¬ ¼ dz « » ¬ ¼ (4.169)
(
)
272 Transport Phenomena in Multiphase Systems
The momentum equation can be obtained by rearranging eq. (4.169), 1d A A A v v v ρA wA AA wA + ρv wv Av wv Ac dz (4.170) 1 d ( p Ac ) P v A =− − τ w − ªα ρA + (1 − α ) ρv º g cos θ ¬ ¼ Ac dz Ac The energy equation balance in the control volume satisfies eq. (4.164), except that the energy storage rate is zero. The energy inflow rate A A includes energy flow into the control volume at surface z, ρ A AA wA
(
)
⋅ hA
′′ and heat transferred into the control volume through the wall, qw Pdz . The energy outflow rate is the energy flow out from the control volume
at surface z+dz and it is + ρv ⋅ hA
v
(
A
+ wA
A2
/ 2 + gz cos θ + ρ v
)
v
Av wv
v
(h
v
v
+ wv
v2
/ 2 + gz cos θ ,
)
Av wv + wA
v
+ gz cos θ ) º dz . Thus, the energy balance for the channel element ¼
becomes
(
A
A2
( h + w / 2 + gz cosθ ) d ( h + w / 2 + gz cosθ ) + dz ª¬ ρ A w / 2 + gz cos θ ) + ρ A w ( h + w /2
ρA
A
AA wA
A
A
A2
A
A
v
v2
A
A
v
v
A
A
A
v
v
v
v2
v
v
v
v
v
(4.171) v2 § ·º wv v v v + ρv Av wv ¨ hv + + gz cos θ ¸ » = 0 ¨ ¸» 2 © ¹¼ Rearranging eq. (4.171) yields the following energy equation: v2 º A2 § · § · w w dª « ρA A AA wA A ¨ hA A + A ¸ + ρv v Av wv v ¨ hv v + v ¸ » ¨ ¨ 2¸ 2 ¸» dz « © ¹ © ¹¼ ¬ ′′ = qw P − ρ A
dª ′′ qw P − « ρ A dz « ¬
A
AA wA
A
§ ¨ hA ¨ ©
A
+
wA 2
A2
· + gz cos θ ¸ ¸ ¹
(
A
AA wA
A
+ ρv
v
Av wv
v
) g cosθ
(4.172)
4.6 An Extension: Porous Media
Transport in porous media is applicable to a wide range of fields, including mechanical, chemical, environmental, and petroleum engineering, as well as geology. Porous media can be found naturally in rocks and sand beds, and also can be fabricated as in wicks and catalytic pellets. They are an essential
Chapter 4 Generalized Governing Equations: Averaging Formulations 273
component in high-technology devices such as fuel cells and heat pipes. A fundamental formulation of the governing equations within a porous media will be presented here. The variation of some of the assumptions made has been quantified in recent studies. For these analyses of the variants within a particular porous media model, see Alazmi and Vafai (2000, 2001), and Nield and Bejan (1999). A porous medium is a solid matrix with a several voids, or pores, which are continuously connected. The voids are filled with one or more fluids that can pass through the medium because the voids are interconnected. The void fraction in the solid matrix is most frequently referred to as the porosity, ε . Conversely, 1 − ε is the volume fraction of the solid matrix, ε sm , that makes the porous medium. If g ( x, y, z, t ) is a geometric function that is identically equal to 1 in a void, and 0 in the solid matrix, the porosity is defined as: 1 ε= (4.173) ³ gdV . ΔV ΔV V is an elemental volume of the porous zone. When there is more than one phase present in a porous medium, there is a porosity of a single phase, ε k , which is simply the volume fraction of a particular phase in an elemental volume. 1 εk = (4.174) ³ gdV ΔV ΔVk The subscript, k, pertains only to phase k. The volume fraction of the solid matrix in an elemental volume is represented by ε sm . Therefore, 1 = ε + ε sm =
Π−1
k =1
¦ ε k + ε sm
(4.175)
To model a porous zone, it is imperative that one be familiar with the different length scales of the zone. The first dimension is the particle or void length scale, d, and the second dimension, L, is the system or porous zone length scale. If d is on the order of L, such as in a very thin porous layer on a heat transfer surface, where the bulk flow is normal to the surface, the porous zone can be directly modeled with minimal assumptions. However, this situation is usually not the case, and more often d L . When d L , direct simulation of transport characteristics in an individual pore is not practical, therefore the local mean velocity is used. Because of this difference in scale, the volume averages defined in eqs. (4.7) and (4.8) are often used to describe a porous system. Figure 4.5 (a) shows how the flow in an individual pore can be very difficult to model when a larger scale is needed, and therefore a volume-averaged velocity, seen in Fig. 4.5 (b), is often more useful. For analyses in a porous system, one of two velocities are often used: the intrinsic average velocity, k Vk , which is the volume-averaged velocity over the volume filled with a particular phase, or the extrinsic average velocity, Vk , which is the extrinsic averaged velocity over an entire volume including a solid. These two velocities are related by the Dupuit-Forchheimer relationship:
274 Transport Phenomena in Multiphase Systems
(a)
(b)
Figure 4.5 Comparison of velocity distributions in an elemental volume in a porous medium: (a) actual velocity, and (b) volume-averaged velocity.
Vk = ε k Vk (4.176) The governing equations for the fluid phase in a porous medium can now be developed. It is assumed that the porous medium is saturated with the fluid phase and that the relative velocity, Vf,rel, only refers to the relative motion of a fluid phase f with respect to the solid matrix. In the following sections, the solid matrix is assumed to be stationary, therefore V f = V f ,rel . In Sections 4.6.1 –
k
4.6.5, the conservation laws are derived for a porous zone by using a volumeaveraging technique for a single fluid occupying the pores. The momentum equations are derived in Section 4.6.2 in a manner in which the viscous interactions between a fluid and the solid matrix are modeled with a property of the porous wick called the permeability, K. The energy equation is derived in Section 4.6.3 in such a way that local thermal equilibrium between the solid matrix and the fluid can be assumed, or local thermal nonequilibrium conditions can be assumed. The second law is presented in Section 4.6.4, and the species equation is presented in Section 4.6.5 for a porous zone. Finally, multiphase transport phenomena in porous media are considered in Section 4.6.6.
4.6.1 Conservation of Mass
The fluid-saturated porous media can be considered as a two-phase system that contains fluid (f) and solid matrix (sm). According to eq. (4.22), the continuity equation for the fluid phase (f) in the porous media is f f f ∂ ε ρf Vf +∇⋅ ε ρf =0 (4.177) ∂t
(
)(
)
Chapter 4 Generalized Governing Equations: Averaging Formulations 275
where
Vf
f
is the intrinsic phase-averaged velocity. The right-hand side is
zero, because there is no phase change between the porous matrix and the fluid. If the fluid phase in the porous media can be assumed as incompressible, the volume-averaged density is equivalent to the density of a single phase,
ρf
f
= ρ f and eq. (4.177) can be simplified as ∂ ερ f + ∇ ⋅ ρ f V f ∂t
(
)
(
)=0
(4.178)
where
Vf = ε Vf
f
is the extrinsic phase-averaged velocity.
For an incompressible flow in a porous media of uniform porosity, the continuity equation is ∇ ⋅ Vf = ∇ ⋅ ε f Vf
(
f
) = 0.
(4.179)
4.6.2 Conservation of Momentum
Since a macroscopic approach is the most feasible way to model transport, there needs to be a way to model the bulk resistance to the flow caused by the porous zone. In 1856, Henry Darcy experimentally measured the resistance to a steady, one-dimensional, gravitationally-driven flow through an unconsolidated, uniform, rigid, and isotropic solid matrix. He came up with a relationship for the pressure gradients/resistance as a function of the dynamic viscosity, f, the extrinsic phase-averaged velocity, and the permeability, K, now known as Darcy’s law. μ 1 (4.180) ρ g − ∇(ε p ) = f V f K ε where the permeability, K, is the resistance to the flow, which is analogous to the thermal conductivity in the Fourier’s law of heat conduction. It has units of meters squared. Darcy’s law can be expanded into vectorial form for an anisotropic solid matrix. μ 1 (4.181) ρ g − ∇(ε p ) = f V f K ε where the permeability, K, becomes a second-order symmetric tensor. When the porous media is isotropic (having equal resistance in all directions), the permeability reduces to a scalar (single value), K. Darcy’s law is valid for Re<1, or the creeping flow regime, where the viscous forces dominate. The Reynolds number of the porous media is defined by the mean volumetric velocity, and the average characteristic length scale of the voids, d .
276 Transport Phenomena in Multiphase Systems
Re =
ρ f Vf d μf
(4.182)
In this flow regime, the wall effects are confined to within one or two particle diameters of the wall. Also, in the region where the flow enters the porous zone, the entrance region is confined to approximately three particle diameters. Therefore flow becomes very closely uniform in one main direction, and the walls bounding the porous zone have minimal effect. As the Reynolds number increases from unity to about 10, the drag smoothly transitions from linear to nonlinear. This linear to nonlinear transition of drag in a porous zone represents the flow moving from the creeping flow to the laminar flow regime. To account for this nonlinearity, when Re>10, Ward (1964) came up with a quadratic drag term dependent of /K1/2. μf Cf 1 ρ g − ∇(ε p ) = Vf + 1 2 ρ f Vf Vf (4.183) K ε K where Cf is a dimensionless drag constant, initially thought by Ward (1964) to be a universal constant, 0.55. Since the drag smoothly transitions from linear to nonlinear when 10>Re>1, the significance of the quadratic drag term can be formed as a linear function of Re. μ C ( Re− 1) 1 (4.184) ρ g − ∇(ε p ) = f V f + 1f 2 ρ f V f V f 9 K ε K It was later found to vary based on the nature of the porous medium. From the differential momentum equation, eq. (3.52), the advective term is ρ f V f ⋅ ∇ V f . From a scaling analysis, the velocity of the fluid is on the order
(
)
of the volume weighted average velocity over the porosity, for a single phase,
Vf ~ Vf / ε .
Therefore, the advective term is on the order of ρ f V f
2
/ε 2 .
A scaling analysis between the advective term and the quadratic drag term gives a ratio proportional to K1 2 /(C f ε 2 L) , where L is the macroscopic characteristic length scale. The ratio of these two terms is usually very small, therefore the advective term often can be assumed to have a negligible effect on the momentum equations. An extension of the Stokes drag force on a sphere in an infinite domain is to include the effects of the neighboring spheres, as with the case of packed beds made of spheres. The Stokes drag force is the drag exerted by a single sphere in an infinite domain on a flow with a Re 1 (negligible inertia). This was accomplished by superimposing the Stokes flow on the Darcy flow, i.e., μ 1 ρ g − ∇(ε p) = f V f − μ ′f ∇ 2 V f (4.185) ε K where the effective viscosity denoted by μ ′f = μ ′f μ f , ε , ℑ is a function of the
(
)
dynamic viscosity, the porosity, and the tortuosity ℑ of the porous media. The
Chapter 4 Generalized Governing Equations: Averaging Formulations 277
tortuosity is a measure of the connectivity of void space in a porous zone. Equation (4.185) is called the Brinkman equation, and it is usually noted that it is valid for high porosity, ε > 0.8 . When the porosity is low, the stresses felt on a fluid in one pore are communicated to the fluid in another pore mainly by pressure because the solid matrix prevents direct viscous interaction of the fluid in separate pores. Even though this equation is only valid for higher porosities, the Laplacian operator of V f is needed when the boundaries of the porous media affect the flow field. An empirical momentum equation can be heuristically obtained if Darcy’s law, eq. (4.180), is combined with the inertial drag component, eq. (4.183), and Brinkman’s equation, eq. (4.185). μf Cf 1 ρ g − ∇(ε p ) = − μ ′f ∇ 2 V f + Vf + 1 2 ρ f Vf Vf (4.186) K ε K However, for a more complete understanding, this equation is compared to a volume-averaged momentum equation derived in Section 4.3.2. If the fluid phase is incompressible, i.e.,
ρf
f
= ρ f , the volume-averaged momentum
equation for the fluid phase, eq. (4.34), becomes f f f ∂ ερ f V f + ∇ ⋅ ερ f V f Vf ∂t f § ˆˆ = ∇ ⋅ ε ′f + ερ f g + Fsm , f − ∇ ⋅ ¨ ερ f V f V f ©
(
(
)( )
f
)
f
· ¸ ¹ Vf
f
(4.187)
where
Vf Vf
f
in eq. (4.34) is replaced by
Vf
f
ˆˆ + Vf Vf
and Xf
is replaced by g, i.e., gravity is the only body force. Since the flow in a porous media is usually laminar because of the characteristic diameter, the product of the ˆˆ velocity deviations, V V , is also small and therefore neglected henceforth.
f
The third term on the right hand side of eq. (4.187) represents interaction between the solid matrix and the fluid phase. For an incompressible fluid, the volume average of the stress tensor, eq. (4.26), becomes f f f f Tº ª ′f = − p f I + μ f «∇ V f + ∇ V f (4.188) » ¬ ¼ Since the local static pressure is needed in the equation of state for an ideal gas, or for the temperature at a liquid vapor interface, the local pressure is more accurately described as the intrinsic phase average pressure, not the extrinsic phase averaged pressure. If the deviation of pressure within a fluid volume ˆ element is small, pk ≈ 0 , then the volume-averaged pressure is
(
)
pf = ε pf
f
=ε p.
(4.189)
278 Transport Phenomena in Multiphase Systems
Substituting eq. (4.188) into eq. (4.187) and considering eq. (4.189), the momentum equation becomes f f f ∂ + ∇ ⋅ ερ f V f = −∇ ( ε p ) + ερ f V f Vf ∂t (4.190) f 2 + ερ f g + Fsm , f μ f ∇ ε Vf
(
(
)( )
)
Using the continuity equation, eq. (4.178), and assuming constant porosity, the momentum equation becomes ª 1 ∂ Vf º Vf + 2 ∇ ⋅ Vf » ρf « ε « ε ∂t » ¬ ¼ (4.191) μf 2 1 = −∇p + ∇ V f + ρ f g + Fsm , f
(
)
ε
ε
Table 4.1 Properties of porous media materials (Nield and Bejan, 1999; Faghri 1995)a Solid Matrix Black Slate Powder Brick Cigarette Cigarette filters Coal Concrete (ordinary mixes) Copper powder (hotcompacted) Cork board Fiberglass Hair (on mammals) Hair felt Leather Limestone (dolomite) Porosity 0.57 – 0.66 0.12 – 0.34 0.17 – 0.49 0.02 – 0.12 0.02 – 0.07 0.09 – 0.34 0.88 – 0.93 0.95 – 0.99 0.56 – 0.59 0.04 – 0.10 8.3×10-10 – 1.2×10-9 9.5×10-14 – 1.2×10-13 2×10-15 – 4.5×10-14 2.0×10-11 – 1.8×10-10 5.2×10-11 3.3×10-10 – 1.5×10-9 2.4×10-11 – 5.1×10-11 (5.6 – 7.7)×104 Permeability K (m2) 4.9×10-14 – 1.2×10-13 4.8×10-15 – 2.2×10-13 1.1×10-9 Effective Pore Radius reff (m)
ΔAsm ΔV
( m −1 )
7×105 - 8.9×105
(1.2 – 1.6)×106
Sand 0.37 – 0.50 (1.5– 2.2)×104 -6 Screen, SST, 200 mesh 0.733 58×10 305×10-6 Screen, Nickel, 50 mesh, -10 0.625 6.63×10 sintered Felt, Sintered, SST 0.916 5.46×10-10 94×10-6 -11 Felt, Nickel, A30 0.815 3.06×10 120×10-6 Powder, Sintered, Nickel 0.597 58×10-6 -11 Beads, Monel, 70-80 Mesh 0.4 7.75×10 96.9×10-6 -16 -12 Sandstone (“oil sand”) 0.08 – 0.38 5×10 – 3×10 Silica grains 0.65 Silica powder 0.37 – 0.49 1.3×10-14 – 5.1×10-14 (6.8 – 8.9) ×105 -13 -11 Soil 0.43 – 0.54 2.9×10 – 1.4×10 Wire crimps 0.68 – 0.76 3.8×10-9– 1.0×10-8 (2.9– 4.0)×103 Reproduced with kind permissions of Springer Science and Business Media, and Routledge/Taylor & Francis Group, LLC.
Chapter 4 Generalized Governing Equations: Averaging Formulations 279
The empirical momentum relationships are heuristically related to the volumeaveraged momentum equation through reasonable observations. In the second term on the right hand side of eq. (4.191), the value of μ f / ε represents μ ′f in eq. (4.186). The fourth term on the right-hand side of eq. (4.190) is replaced by the Darcian relationship and the quadratic drag term. The Darcian drag term comes from the deformation stress caused by the solid/fluid interaction, and the quadratic drag term comes from the pressure coefficient associated with the fluid flowing around the solid matrix. ª 1 ∂ Vf º Vf + 2 ∇ ⋅ V f » = −∇p + ρf « ε « ε ∂t » ¬ ¼ (4.192) μf 2 μf ρ f Cf ∇ Vf + ρ f g − Vf − 1 2 Vf Vf ε K K This equation reduces to eq. (4.186) if the inertia terms are negligible. In eq. (4.192), the viscous interaction between the fluid and the solid matrix is replaced by Darcy loss and inertial loss terms. The shear losses between the fluid and the solid matrix are accounted for without having any information on the characteristics of the flow in each pore. The permeability, K, can be obtained by empirical relations, which are well tabulated for various types of porous media. For packed beds of spherical particles of diameter d, the permeability is d 2ε 3 K= (4.193) 150(1 − ε ) 2 Further details on permeability and porosity of various porous materials, including porous wicks in heat pipes, can be found in Kaviany (1995) and Faghri (1995). For convenience, the porous zone properties for a variety of applications are presented in Table 4.1, including K, ε , ΔAsm / ΔV and reff. Example 4.7 A cylindrical tube has a porous zone inside it. The porous zone consists of N smaller cylindrical tubes within it as shown in Fig. 4.6. The large cylinder has a radius of R0, while the smaller internal cylinders all have a radius of r0. Derive the permeability from the momentum equation derived by the volume averaging technique. Assume that the flow is incompressible, a fully-developed flow within each tube, and that the volume-averaged velocity, Vk , is constant throughout the larger cylinder. Solution: There are two parts in solving this problem, (1) knowing the frictional losses in each pore, and (2) then substituting these losses into the simplified porous media momentum equation, eq. (4.180). Since the flow is steady, incompressible and one-dimensional, continuity is already specified by having a uniform flow through the cylinder. Also
280 Transport Phenomena in Multiphase Systems
the volume averaged momentum equation, eq. (4.25), can be simplified by: f f ∂ 1 f +∇⋅ ερ f Vf Vf = ∇⋅ ε τ − ερ f Vf ³ ρ f Vf ( Vf − VI ) ⋅ ndA ∂t ΔV ΔAsm 0 0 0 1 + ³ τ ⋅ ndAk + ερ f X f ΔV Ak 0 With constant properties, and the assumed fully developed profile, the above equation becomes
(
)(
)(
)
(4.194) ³ D f ⋅ n dA ∂x ΔV ΔAsm The fully developed velocity profile in a single tube is: r02 − r 2 § dp f · Vf ⋅ i = (4.195) ¨− ¸ 4μ © dx ¹ where i is the unit vector in the x-direction. From the velocity profile, we can obtain the rate of deformation, D f ⋅ n , with respect to the smaller
ε
∂ pf
f
=
2μ f
(
)
(
)
cylindrical tubes surface. 1 ∂V f ⋅ i r0 dp f (4.196) Df ⋅n = = 2 ∂r 4 μ dx Integrating this value over all of the pore wall areas to give the total frictional loss.
§ r0 · dp f (4.197) ³ΔAsm D f ⋅ n dA = N ¨ R ¸ dx ΔV © 0¹ The extrinsic average velocity is equal to the area average velocity. Nr04 § dp f · 2N r V f ⋅ i = 2 ³00 V f ⋅ irdr = (4.198) − ¸ 2¨ 8μ R0 © dx ¹ R0 2μ f N
(
)
2
x, i R0
2r0
Figure 4.6 Cross-sectional and side view of a cylinder with cylindrical pores.
Chapter 4 Generalized Governing Equations: Averaging Formulations 281
Now the frictional loss can be converted to a function of substituted into the original momentum equation, eq. (4.194).
d pf dx
f
V f , and
ε
=
−8μ f V f ⋅ i r02
(4.199)
The porosity is:
Nr02 2 R0 To maintain the form of Darcy’s law,
ε=
(4.200)
dx K Therefore, the permeability is Nr 4 K = 02 8 R0
−
d pf
f
=
μ
Vf ⋅ i =
2 8R0 μ Vf ⋅ i Nr04
(4.201)
(4.202)
Example 4.8 When using a Darcy’s law in a porous media, the effect of a no-slip boundary condition cannot be accounted for because of the order of the equation. The shear stress caused by a wall (Fig. 4.7) can be captured using Brinkman’s equation. Use the analysis of a semi-infinite, fully developed velocity in a porous zone in no gravitational field to determine the criteria when the wall effects should be incorporated in the governing equations. Develop a solution using a no-slip wall condition and use this solution to determine the significance of the wall for different porous properties. Also, find the volume average velocity profile of a flow through a porous media bounded by parallel plates under a given pressure gradient. Use no-slip boundary conditions, and compare the mean flow rate to the mean flow rate using Darcy’s law. Solution: Brinkman’s equation in a one-dimensional porous zone is: 2 § uf dp 1 ∂ uf · ¨ ¸ − = μf − (4.203) ¨K dx ε ∂y 2 ¸ © ¹ where u f The extrinsic phase-averaged velocity in the x-direction
and the effective viscosity is assumed to be μ ′ = μ / ε . conditions at the wall and at the free stream of: y = 0, u f = 0 , u f = 0
y = ∞, ∂ Vf ∂y =0
With boundary (4.204) (4.205)
282 Transport Phenomena in Multiphase Systems
Elemental Control Volume
x, uf y, vf
Figure 4.7 Porous region near a wall
Assume that the volume average velocity at the wall is zero (knowing it may not be correct). The solution to the semi-infinite Brinkman’s equation is: μf uf § § ε ·· (4.206) − = ¨1 − exp ¨ − y ¸¸ ¨ ¨ dp K ¸¸ © ¹¹ © K dx which is plotted in Fig. 4.8. It can be seen in Fig. 4.8 that the wall effects are confined to a distance of approximately 6 K / ε . The average pore length in the ydirection is about K / ε ; therefore, the wall effects will only be noticeable in cells that are within approximately six pore lengths from the wall. For computational considerations, if the grid is not resolved within six average pore lengths, then there is no advantage to using Brinkman’s equation over Darcy’s law. To quantify the overall size of the solution to determine what scale is large enough such that wall effects can be neglected, we must solve Brinkman’s equation for a bounded flow between two parallel plates ( 0 ≤ y ≤ h ).
2 § uf dp 1 ∂ uf ¨ − = μf − ¨K ε ∂y 2 dx ©
· ¸ ¸ ¹
(4.207)
Chapter 4 Generalized Governing Equations: Averaging Formulations 283
−
μf uf
K dp dx
y ε /K Figure 4.8 Solution to Brinkman’s equation in a semi-infinite domain.
The boundary conditions are: y = 0, u = 0
(4.208)
y=h, u =0 (4.209) The general solution to this equation after applying boundary conditions in eq. (4.208) and (4.209) is: §ε· § K dp § ε ·· (4.210) uf = − ¨ 1 − A exp ¨ ¨ K y ¸ − B exp ¨ − K y ¸ ¸ ¸ ¨ ¸¸ μ f dx ¨ © ¹ © ¹¹ © where the constants A and B are:
ª § § ε ·º ε · ºª § ε · A = « exp ¨ − h ¸ − 1» « exp ¨ − h ¸ − exp ¨ ¨ ¨ ¨ K h ¸» ¸ K ¸ »« K¸ « © ¹ ¼¬ © ¹ © ¹» ¬ ¼
−1
(4.211)
ª § ε ·º ª § § ε ·º ε· B = «1 − exp ¨ h ¸ » «exp ¨ − h ¸ − exp ¨ ¨K¸ ¨ ¨ K h ¸» ¸ K¸ « © ¹» « © ¹ © ¹» ¬ ¼¬ ¼ The average velocity found by integrating u f between the plates is 1h ³ u f dy h0 Therefore, the nondimensional average velocity is: μu 2 ( exp ( −η ) − 1) ( exp (η ) − 1) − =1− dp η exp ( −η ) − exp (η ) K dx uf =
−1
(4.212) over the distance
(4.213)
(4.214)
284 Transport Phenomena in Multiphase Systems
Figure 4.9 Non-dimensional average velocity vs. non-dimensional distance.
where η is equal to the non-dimensional distance, h ε / K . Equation (4.214) is presented in Fig. 4.9. By analyzing Fig. 4.9, it can be noted that the wall effects are insignificant when the overall length of the system η is greater than 103, and therefore Darcy’s law is a good approximation. However, when η is less than 103, Brinkman’s equation should be applied. If η is on the order of 1, the volume averaging approach is not good, and the local velocities should be modeled.
4.6.3 Energy Equation
The derivation of the energy equation in a porous medium is similar to that of continuity and momentum equations. The volume-averaged energy equation (4.44) is valid for both fluid phase and the solid matrix. Assuming the fluid phase is incompressible, and neglecting the viscosity dissipation, the energy equation for the fluid phase becomes f f f ∂ ′′′ + ∇ ⋅ ερ f V f = −∇ ⋅ q′′ + q′′′ + qsm, f (4.215) ερ f h f hf f f ∂t
(
)(
)
where
Vf hf
f
in eq. (4.44) is replaced by
Vf
f
hf
f
in eq. (4.215), i.e.,
the product of the deviations were neglected. matrix to the fluid phase.
′′′ qs m , f
is heat transfer from solid
Chapter 4 Generalized Governing Equations: Averaging Formulations 285
By applying the continuity equation (4.178), eq. (4.215) becomes ′′′ (4.216) = −∇ ⋅ q′′ + q′′′ + qsm, f f f ∂t The volume-averaged energy equation in the solid matrix is: sm sm ∂ hsm ′′′ ′′′ = −∇ ⋅ q′′ + qsm − qsm , f ε sm ρ sm (4.217) sm ∂t where the convection term is dropped since the solid matrix is in the solid state and has no velocity. The thermal boundary conditions at the fluid-solid interface, Asm, are: T f = Tsm (4.218)
ερ f
∂ hf
f
+ ερ f V f
f
⋅ ∇ hf
f
q′′ ⋅ n f = q′′ ⋅ n sm . f sm
(4.219)
These boundary conditions between the fluid and the solid matrix are what make ′′′ up the heat transfer term between the fluid and solid, qsm, f . Equations (4.218) and (4.219) are very general, and no assumptions about the thermal equilibrium between the fluid and the solid matrix have been made. If the fluid and the solid matrix are considered to be in local thermal equilibrium, Tf
f
= Tsm
∂ hf
sm
= T , the energy equations for the solid and the fluid can be
∂ hsm
sm
added to find the total porous media energy equation:
f
ερ f
∂t ∂t ′′′ = −∇ ⋅ q′′ − ∇ ⋅ q′′ + q′′′ + qsm f sm f
+ ε sm ρ sm
sm
+ ερ f V f
f
⋅ ∇ hf
f
(4.220)
The energy equation can be further simplified by assuming that the enthalpy is only a function of temperature and constant specific heat: ∂T ′′′ + ( ρ c p ) f V f ⋅ ∇T = ∇ ⋅ keff ∇T + qeff ( ρ c p )eff (4.221) ∂t The effective heat capacity, thermal conductivity and heat generation rates per unit volume, denoted by the subscript “eff,” are defined by (4.222) ( ρ c p ) = ε f ( ρ c p ) f + ε sm ( ρ c p )sm ,
eff
keff = ε keff , f + ε sm keff , sm , ′′′ qeff = ε q′′′ f
f
(4.223) . (4.224)
′′′ + ε sm qsm
sm
where keff , f and keff , sm are, respectively, the effective thermal conductivity of the fluid and solid-matrix in the porous media, and both of them depend on porosity and pore structure in the porous media. When the thermal conductivities of the fluid and the solid matrix are close to each other, k f k sm , one can assume that keff , f ≈ k f and keff , sm ≈ k sm , and use eq. (4.223) to evaluate the effective thermal conductivity. In a case where thermal conductivities of the
286 Transport Phenomena in Multiphase Systems
fluid and the solid matrix differ significantly, the following equation should be used (Hadley, 1986): keff ε f 0 + (1 − ε f 0 ) ksm / k f = (1 − α 0 ) 1 − ε (1 − f 0 ) + ε (1 − f 0 ) kf (4.225) 2(1 − ε )(k sm / k f ) 2 + (1 + 2ε )k sm / k f + α0 (2 + ε )k sm / k f + 1 − ε where
f 0 = 0.8 + 0.1ε −4.898ε ° log α 0 = ®−0.405 − 3.154(ε − 0.0827) ° −1.084 − 6.778(ε − 0.298) ¯
(4.226) 0 ≤ ε ≤ 0.0827 0.0827 ≤ ε ≤ 0.298 0.298 ≤ ε ≤ 0.580 (4.227)
The above energy equation, eq. (4.221), is only valid when the fluid and the solid matrix are assumed to be in local thermal equilibrium. While local thermal equilibrium is an often used hypothesis for heat transfer in porous media, it is well recognized that local equilibrium between different phases in a system of solid and fluid cannot be achieved when thermal properties for different phases differ widely, or during rapid heating or cooling. The fluid and solid-matrix energy equations need to be kept separate, and this situation is referred to as local thermal non-equilibrium. In this situation, heat transfer from the solid matrix to ′′′ the fluid phase, qsm , f , in eqs. (4.216) and (4.217) can be modeled as a convective boundary condition, with a heat transfer coefficient of hc . When Fourier’s law of heat conduction is applied, the energy equations for the fluid and solid matrix are § ∂T f · (ρcp ) f ¨ ε + V f ⋅ ∇T f ¸ © ∂t ¹
ΔAsm (4.228) Tsm − T f ΔV ∂T ΔA ′′′ ε sm ( ρ c p ) sm sm = ∇ ⋅ (ε sm keff , sm ∇Tsm ) + ε sm qsm − hc sm Tsm − T f (4.229) ∂t ΔV Equations (4.228) and (4.229) are more accurate representations of the thermal field in a porous zone. However, they require an empirical relationship for the heat transfer coefficient between the solid matrix and the fluid phase. = ∇ ⋅ (ε keff , f ∇T f ) + ε q′′′ f
f
+ h f , sm
(
)
(
)
4.6.4 The Second Law of Thermodynamics
The second law of thermodynamics can be derived in the same manner as the first law. The volume-averaged second law of thermodynamics, eq. (4.50), is valid for both the fluid phase and the solid matrix. If the fluid phase is incompressible, the second law of thermodynamics becomes
Chapter 4 Generalized Governing Equations: Averaging Formulations 287
∂ ερ f s f ∂t
(
f
) + ∇ ⋅ (ερ
f
f
Vf
f
sf
f
)+∇⋅
Vf
f
q′′ f
Tf sf
q′′′ f Tf
f
−
q′′′ f Tf
′′′ − ssm , f > 0 (4.230)
where
Vf s f
is approximated as
. Applying the continuity
equation (4.178), the second law of thermodynamics for the fluid phase becomes
ερ f
∂ sf ∂t
f
+ ρ f v f ∇⋅ sf
f
+∇⋅
q′′ f Tf
−
′′′ − ssm , f > 0
(4.231)
For the solid matrix, we have sm ∂ ssm q′′ q′′′ ′′′ + ∇ ⋅ sm − sm + ssm , f > 0 ε sm ρ sm ∂t Tsm Tsm
(4.232)
The second law of thermodynamics in the fluid and solid matrix can be summed to obtain:
§ q′′ · § q′′′ · + ∇ ⋅¨ ¸ − ¨ ¸ > 0 ∂t ∂t © T ¹eff © T ¹eff (4.233) where the effective entropy generation due to heat conduction and internal heat generation are defined as the sum of the entropy generated in the solid matrix and the fluid.
ερ f
∂ sf
f
+ ε sm ρ sm
∂ ssm
sm
+ ρ f Vf ∇ ⋅ s f
f
q′′ § q′′ · f ¨ ¸ =ε Tf © T ¹eff
q′′′ § q′′′ · f ¨ ¸ =ε Tf © T ¹eff
f
+ ε sm
f
q′′ sm Tsm
′′′ qsm Tsm
sm
(4.234)
sm
+ ε sm
(4.235)
4.6.5 Species
The volume-averaged species eq. (4.54) can be readily applied to the fluid phase in the porous media: (4.236) = −∇ ⋅ J f ,i + m′′′,i f ∂t Applying the continuity equation (4.178) and the definition of the mass fraction, eq. (1.77), the averaged-species equation becomes
f ,i
∂ ε ρ f ,i
(
f
) +∇⋅ ε ρ (
f
f
Vf
f
)
ερ f
∂ ω f ,i ∂t
+ ρ f v f ⋅ ∇ ω f ,i
f
= −∇ ⋅ J f ,i + m′′′,i f
(4.237)
288 Transport Phenomena in Multiphase Systems
4.6.6 Multiphase Transport in Porous Media
The framework for multiphase transport in a porous zone (see Fig. 4.10) is already laid out through the single phase transport in porous media. When the porous medium is saturated with both liquid and vapor, the governing equations for liquid and vapor as well as solid matrix must be specified. In addition, the solid-liquid-vapor interactions play an important role in heat and mass transfer applications. The governing equation for the liquid and vapor phase must also consider possible phase change between liquid and vapor. Chapter 5 is entirely dedicated to problems of involving solid-liquid-vapor interfacial phenomena. Problems that directly solve the solid-liquid-vapor phenomena are typically applied to a single pore. The premise for both the multifluid model (MFM) and the multiphase mixture model (MMM) for solving multiphase porous media are the volumeaveraged Navier-Stokes equations. The major assumption for these models used in a porous medium is that the flow is considered noninertial. One of the benefits of a noninertial flow is that it behaves very well and is in the laminar flow regime. Therefore, the deviatoric values from the volume-averaged equation are neglected without compromising the accuracy of the solution. This approximation is reasonable as the pore size is very small, which makes viscous affects dominate. The MFM solves both phases separately, along with the details of each phase, such as temperature, and species mass fraction. The MMM adds the governing equations for each phase, and converts the flow variables into their mass-averaged values counterparts. Solid matrix
Liquid
Vapor
Figure 4.10 Typical elemental volume with periodic dispersion of phases.
Chapter 4 Generalized Governing Equations: Averaging Formulations 289
Multi-Fluid Model (MFM)
In this model, there are numerous liquid/vapor interfacial regions within each elemental pore; therefore, the average change of these interfacial characteristics is considered as a model for this type of problem. The transport equations for a single-component substance coexisting in the vapor phase, v, and the liquid phase, A , are given. In these equations, the liquid and vapor are considered to be in local thermal equilibrium. However, the liquid and vapor phases are in local thermal non-equilibrium with the solid. Whitaker (1977) developed governing equations for the periodic dispersion of phases within a porous zone, as seen in Fig. 4.10, The continuity equation, (4.178), for the vapor and liquid phases become ∂ (4.238) (ε v ρv ) + ∇ ⋅ ( ρv Vv ) = mA′′′v ∂t m′′′ ∂ε A (4.239) + ∇ ⋅ VA = − Av ∂t ρA where the density in the liquid phase is assumed constant. The volumetric mass ′′′ transfer term, mAv , with units of kg/m3-s, represents the evaporation rate when it has a positive value and the condensation rate when it has a negative value. When gravity is the only body force, the momentum equation (4.190), for the vapor and liquid phases becomes ª 1 ∂ Vv § V ·º V + v ⋅ ∇ ¨ v ¸» ρv « εv « ε v ∂t © ε v ¹» ¬ ¼ (4.240) Cf μv 2 μv v = −∇ pv + ∇ Vv − Vv − 1 2 ρv Vv Vv + ρv g Kv εv K
ρA «
ª 1 ∂ VA «εA ¬ ∂t
+
VA
εA
+ ρA g » ¼ The first two terms in the bracket on the right-hand side of eq. (4.241) represent the pressure gradient due to the change in capillary pressure, pc, as a function of volume fraction of liquid and temperature. The two coefficients in eq. (4.241) are defined as ∂p Cε = c (4.242) ∂ε A ∂pc CT = (4.243) Av ∂ TAv + ∇ pv
μ =− VA + A ∇ 2 VA − ªCε ∇ε A + CT ∇ T f « ¬ KA εA
μA
§ V ·º ⋅ ∇ ¨ A ¸» © ε A ¹» ¼
(4.241)
f vº
290 Transport Phenomena in Multiphase Systems
where the capillary pressure is defined as 2σ v pc = Av = pv − pA reff
A
(4.244)
and reff is the effective radius of curvature of the interface between the vapor and the liquid. For more information on interfacial phenomena, the reader is referred to Chapter 5. The effective radius is a function of the volume fraction of the liquid and the pore geometry. The energy equation for the fluid (liquid and vapor) phase, where TA
A
= Tv
A
v
= Tf
f
, is
v
(ε ρ c
A pA
+ ε v ρv
c pv
)
∂ Tf ∂t
f
+ ρA c p A v A + ρv
(
v
c pv v v ⋅ ∇ T f
f ΔAsm ,A sm Av = ∇ ⋅ ª ε A keff ,A + ε v keff ,v ∇ TAv º + hsm,A Tsm − Tf ¬ ¼ ΔV f ΔAsm ,v sm + hsm,v Tsm − Tf (4.245) ΔV where hc , sm ,A and hc , sm ,v are convective heat transfer coefficients from solid-
(
)
(
)
(
)
f
)
matrix to liquid and vapor phases, respectively. The energy equation for the solid matrix is
(ε sm ρ smc p,sm )
−hsm,A Tsm
∂ Tsm ∂t − Tf
sm
= ∇ ⋅ ε sm keff , sm ∇ Tsm
f
(
sm
)
f
− hsm ,v Tsm − Tf ΔV ΔV The volume fractions of all the phases must satisfy the following constraint: ε A + ε v + ε sm = 1 (4.247) The equation of state for the vapor phase, assuming it is an ideal gas, is pv
v
(
sm
)
ΔAsm ,A
(
sm
)
ΔAsm ,v
(4.246)
= ρv
v
Rg T f
f
(4.248)
The thermodynamic relation for saturation pressure is ª § ·º ½ 2σ Av h 1 1° v ° (4.249) pv = p0 exp ®− « + lv ¨ − ¸» ¾ « r ρ R T f Rg ¨ T f T0 ¸ » ° « eff A g f ©f ¹» ° ¼¿ ¯¬ where T0 is saturation temperature corresponding to a reference pressure p0 . The above governing equations for multiphase transport in porous media are based upon the assumption that the liquid and vapor phases are in thermal equilibrium. Gray (2000) addressed macroscale equilibrium conditions for liquidgas two-phase flow in porous media. Duval et al. (2004) used the method of volume averaging to derive a three-temperature macroscopic model that
Chapter 4 Generalized Governing Equations: Averaging Formulations 291
considered local thermal nonequilibrium among the three phases. A closed form of the evaporation rate at the macroscopic level can be obtained depending on the macroscopic temperatures and the effective properties. Many multiphase problems have a clear liquid/vapor interface within the porous zone. These problems can be dealt with by solving the phases as separate solutions with boundary conditions between the phases. Example 4.9 demonstrates this class of problem.
Example 4.9 A heat pipe uses a thermal cycle in which a liquid is evaporated from one end, and the vapor is condensed at the other end. The liquid motion is caused by capillary action through a wick, which is a porous zone (see Fig. 4.11). The maximum pressure drop that the wick can handle is based on the effective radius, reff, by Δpmax = 2σ / reff . If the
effective radius is 2.5×10-5 m, calculate whether the wick has enough pumping capability to compensate for 10 kW/m2 of heat entering the evaporator and leaving the condenser. The adiabatic section, La, has a length of 0.5 m, while the evaporator and condenser sections have a length of 0.1 m (Le=Lc). The permeability, K, of the wick is 5×10-9 m2, while the wick thickness, wick, is 0.005 m. The porosity of the wick is 0.3. The working fluid is water, and the only effect the vapor phase has on the wick is that the pressure applied to the interface is 1 atmosphere. Assume that the inertial losses in the wick are negligible, and that the fluid flow in the wick is one-dimensional. Also assume that the wick is thin enough that conduction through the wick can be considered onedimensional, and that the flow of the liquid has a negligible effect on the temperature distribution in the wick.
Solution: Because the saturation temperature of water at 1 atm is 100 °C, the density of the liquid is ρv =958.77 kg/m3, the dynamic viscosity is 0.000279 N-s/m2, and the latent heat of vaporization is hAv = 2251.2 kJ/kg. The surface tension is σ = 0.05891 N/m.
x, u
Condenser Lc
Adiabatic/Transport Section La
Evaporator Le
Porous Zone
′ qout qi′n
wick
Figure 4.11 Evaporator, condenser and adiabatic sections of a wick in a heat pipe.
292 Transport Phenomena in Multiphase Systems
Since the volume-averaged velocity is uniform across the wick, the volume-averaged continuity equation, eq. (4.22), can be integrated over the thickness of the wick and reduced to: du ′′ = mvA ρAδ wick (4.250) dx ′′ where mvA is the condensation mass flux (negative for evaporation). The mass flux is related to the heat flux by the latent heat: q′′ ′′ mvA = (4.251) hAv where q′′ is the heat flux. Therefore, du q′′ ρAδ wick (4.252) = dx hAv Multiplying both sides by dx, and then integrating, the volume-averaged velocity become: q′′x u= , x < Le (4.253) ρδ wick hlv q′′Le , Le < x < Le + La (4.254) u= ρδ wick hlv q′′ u= (4.255) ( 2 Le + La − x ) , x > Le + La ρδ wick hlv Darcy’s law can now be applied: dpA μ =− u (4.256) dx K The overall pressure drop in the liquid needed can be obtained by integrating Darcy’s law with respect to x, μ q′′x 2 ΔpA = , x < Le (4.257) K 2 ρδ wick hAv
§ L2 · e (4.258) ¨ Le x − ¸ , Le < x < Le + La K ρδ wick hAv © 2¹ L · x2 · q′′ § μ § ΔpA = Le ( 2 x − Le − La ) + La ¨ x − a ¸ − ¸ , x > Le + La ¨ K ρδ wick hAv © 2¹ 2¹ © (4.259) The pressure vs. distance along the wick is plotted in Fig. 4.12. The maximum difference in pressure with a value of 3102.3 Pa is at the end of the evaporator. The maximum pressure drop the wick can handle is 2σ 2 × ( 0.05891N/m ) Δpmax = Δpv + ΔpA = = = 4712.8Pa reff 2.5 × 10−5 m
ΔpA =
μ
q′′
Chapter 4 Generalized Governing Equations: Averaging Formulations 293
0
-500
-1000
-1500
Δp
-2000 -2500 -3000
Condensor
Adiabatic Section
Evaporator 0.5 0.6 0.7
-3500
0
0.1
0.2
0.3 x (m)
0.4
Figure 4.12 Liquid pressure drop (Pa) in the wick section of a heat pipe.
Since the maximum pressure drop the wick can handle is greater than the calculated pressure drop in the liquid, assuming the pressure drop in the vapor is zero, the wick has enough pumping capacity. The basic equations for the MFM model can also be laid out in terms of phase saturation. The volume averaged continuity equation for phase k in terms of phase saturation is: ∂ (4.260) (ε sk ρ k ) + ∇ ⋅ ε sk ρk Vk k = mk′′′ ∂t The porosity of the porous medium is ε , the phase saturation is sk , the
(
)
and the volumetric mass production ′′′ due to phase change from all other phases to phase k is mk . The summation of the phase saturation is unity, ¦ sk = 1 .
k =1 Π
intrinsic phase average velocity is
Vk
k
The phase saturation is the fraction of a
pore occupied by phase k; therefore, the product of the phase saturation and the porosity is the total volume fraction of that phase, ε k . ε k = ε sk (4.261) The summation of the total volume fraction of each phase and the solid matrix is unity,
k =1
¦ ε k + ε sm = 1 .
Π
The momentum equation in a porous medium
294 Transport Phenomena in Multiphase Systems
for phase k is Darcian if inertia and macroscopic shear effects (shear stresses between pores) are neglected. KK rk k ε sk Vk = − (4.262) ( ∇pk − ρ k g )
μk Equation (4.262) can be directly inserted into the continuity equation to get a Laplace-type equation for the pressure of phase k. The volume-averaged
k
species equation for the species mass fraction, ωk ,i
, is:
k k ∂ k ′′′ = − ∇ ⋅ J k ,i + mk ,i (4.263) ε sk ρ k ωk ,i + ∇ ⋅ ε sk ρ k Vk ωk ,i ∂t The diffusive mass flux of species i is J k ,i , and the species generation rate
(
)(
)
′′′ due to phase change or chemical reaction is mk ,i .
The products of chemical
reactions in phase k are also in phase k; therefore the sum of the species generation term is simply the mass production due to phase change: ′′′ ′′′ ¦ mk ,i = mk
i =1 N
(4.264)
Also, the summation of the species generation for one component overall the phases is simply the species generation rate due to chemical reactions, mi′′′ , only because phase change processes do not change the species.
k =1
′′′ ¦ mk ,i = mi′′′
Π
(4.265)
The final equation is the energy equation for phase k, and it is: ∂ k k k ε sk ρ k hk + ∇ ⋅ ε sk ρ k Vk hk ∂t
(
)
(
)
(4.266)
′′′ ′′′ = − ∇ ⋅ q′′,i + qk + mk hk k
k
′′′ + qk , E
The second term on the right hand side represents the total volumetric heat transfer rate from all other phases to phase k, and the third term on the right hand side is the heat added through the mass production of phase k through phase change. The last term is the heat generated in phase k to an external heat source such as radiation or electrical current. The summation of the second and third terms over all the phases is zero, because of the interfacial energy balance.
k =1
′′′ ′′′ ¦ ( qk + mk
Π
hk
k
)=0
(4.267)
Example 4.10 Transport phenomena in the gas diffusion layer of a proton exchange membrane fuel cell (PEMFC; see Section 1.6.1) can be simplified into a one-dimensional steady-state problem, as shown in Fig. 4.13 (Nam and Kaviany, 2003). Gaseous oxygen is supplied from the channel and diffused into the gas diffusion layer. Water vapor enters at the boundary between the catalyst layer and the diffusion layer, and
Chapter 4 Generalized Governing Equations: Averaging Formulations 295
Channel
Gas Diffusion Layer
Cathode Catalyst Layer
m′′,O2 g
T∞ cO2 ,∞ cH2 O,∞
′′ qk
′′ ′′ mv , H 2O , mA , H 2O y y=0
′′ mv , H 2O y=L
Figure 4.13 One-dimensional model of a gas diffusion layer
condenses into liquid water in the gas diffusion layer. The condensate is then depleted to the channel. The gas diffusion layer also receives chemical reaction heat from its interface with the catalyst layer. Specify the governing equations and corresponding boundary conditions for transport phenomena in the gas diffusion layer.
Solution: Assuming that different phases in the gas diffusion layer are in thermal equilibrium, and that the transport phenomena in the gas diffusion layer are steady-state, the energy equation, eq. (4.102), is dT· d§ ′′′ (4.268) ¨ keff ¸ = hAv mvA ,H2O dy © dy ¹ is effective thermal conductivity, T is volumewhere keff
′′′ averaged temperature, mvA ,H 2O is the volumetric condensation rate of water (m3/kg-s). The boundary condition at the interface between the channel and the gas diffusion layer is at a prescribed temperature, i.e., T = 0, y = 0 (4.269) At the interface between the catalyst and gas diffusion layers, we have ∂T ′′ k = qk , y = L (4.270) ∂y
296 Transport Phenomena in Multiphase Systems
The species conservation equation, eq. (4.110) for oxygen uses a Fickian diffusion term and an advection transport term. Fickian diffusion assumes that the diffusive transport of species is driven by the concentration gradient of that species. dxO2 J· d§ (4.271) + xO2 M O2 ¦ i ¸ = 0 ¨ −cM O2 Dm,O2 dy © dy i Mi ¹ where c = p( Rg T ) is total molar concentration, xi is molar fraction of species i, Mi is molecular mass for species i, and effective mass diffusivity of oxygen. The boundary conditions for eq. (4.271) are cO2 = cO2 ,∞ , y = 0 −cM O2 Dm,O2 where cO2 = cxO2 = J O2 , y = L dy is the mole concentration of oxygen.
dxO2 Dm,O2
is the
(4.272) (4.273)
The species conservation equation for water vapor also has a diffusive and advection transport term, plus a term representing the condensation rate. dxH 2O J· d§ ′′′ + xH 2O M H2O ¦ i ¸ = −mvA ,H 2O (4.274) ¨ −cM H2O Dm,H2O dy © dy i Mi ¹ subject to the following boundary conditions cH 2O = cH2O,∞ , y = 0
cM H2O Dm,H 2O dxH 2O
(4.275)
where cH 2O = cxH 2O
= J H2O , y = L (4.276) dy is the mole concentration of water vapor.
pv ,H 2O − psat ,H 2O (T ) Rg T
′′′ The volumetric condensation rate, mvA ,H 2O , is obtained by ′′′ mvA ,H 2O = M H 2Oγ (4.277)
is the volumetric condensation coefficient, and where γ pv ,H2O and psat ,H2O (T ) are partial pressure of the water vapor and saturation pressure of water at temperature T, respectively. The liquid flow in the gas diffusion layer is driven by capillary force, assuming that the gas phase pressure is constant, and the Reynolds number is less than unity, i.e., KK rA § dpc · dS º dª ′′′ (4.278) « − ρA » = mvA ,H 2O dy ¬ μA ¨ dS ¸ dy ¼ © ¹
Chapter 4 Generalized Governing Equations: Averaging Formulations 297
where S = ( s − sir ) /(1 − sir ) is reduced water saturation (s is saturation defined as the fraction of pore space occupied by liquid water, sir is irreducible saturation), K is permeability, and K rA = S 3 is relative permeability. The relative permeability accounts for the increased resistance to the flow in the liquid when the liquid saturation decreases.
Multiphase Mixture Model (MMM)
The general procedure to relate the MFM to the MMM is to sum each equation over all of the phases. The introductions of the mass-averaged density and velocity, as well as the relative mobility, are given where appropriate. The MMM was initially developed by Wang and Cheng (1996) and was applied to model multiphase flow in fuel cells. The continuity equation summed over two phases, k and j, is: jº ∂ª k ª º (4.279) ¬ε sk ρ k + s j ρ j ¼ + ∇ ⋅ «ε sk ρ k Vk + s j ρ j V j » = 0 ¬ ¼ ∂t Defining mixture density and the mixture velocity ρ = sk ρ k + s j ρ j (4.280)
(
)
(
)
ρ V = ε sk ρ k Vk
k
+ ε s j ρ j Vj
j
(4.281)
the continuity equation can be rewritten as ∂ (4.282) (ερ ) + ∇ ⋅ ( ρ V ) = 0 ∂t The momentum equation can be summed over both phases. ªK º K rj j k ε sk ρ k Vk + ε s j ρ j V j = − K « rk ( ∇pk − ρ k g ) + ∇p j − ρ j g » νj « νk » ¬ ¼ (4.283) A capillary pressure that relates the pressure in phase k to the pressure in phase j is introduced: pc = p j − pk (4.284)
(
)
The capillary pressure is often expressed as the Leverette function (Leverette, 1940). This function relates the capillary pressure to the wetting phase saturation, sw . A phase is said to wet the porous material if the contact angle that phase makes with it is less than 90° . Therefore, if θ k < 90° , the wetting phase is phase k, sw = sk ; if θ j < 90° , then phase j is the wetting phase,
sw = s j .
§ε · pc = σ cos θ ¨ ¸ ©K¹
1/ 2
ª1.417 (1 − s ) − 2.120 (1 − s )2 + 1.263 (1 − s )3 º (4.285) w w w ¬ ¼
298 Transport Phenomena in Multiphase Systems
This function was developed to describe the capillary pressure in soils engineering; however, its use has been extended to other technology such as fuel cells. The reason for the over usage of the Leverette function is the lack of functions to describe the capillary pressure for other types of porous media. With the definition of the mixture velocity and the capillary pressure, the mixture momentum equation can be rewritten as: ª§ K §ρ K K· ρK · º K ρ V = − K «¨ rk + rj ¸ ∇p j − rk ∇pc − ¨ k rk + j rj ¸ g » (4.286) ¨ νk νj ¸ νk νj ¸ ¼ «¨ ¹ © ¹» ¬© ν k Note that the capillary pressure gradient can be reduced into its relative components. k ∂p ∂pc N −1 ∂σ jk ∂pc ∂σ jk ∇pc = c ∇sk + ∇ ωk ,i + ∇T (4.287) ¦ k ∂sk ∂σ jk i =1 ∂ ω ∂σ jk ∂T
k ,i
Now the mixture kinematic viscosity, ν , and relative mobility, λk , are introduced. §K K rj · ν = ¨ rk + ¸ ¨ νk νj ¸ © ¹ K rk λk = ν
−1
(4.288) (4.289)
νk The momentum equation, eq. (4.286), can be rewritten as:
Kª ∇p j − λk ∇pc − λk ρ k + λ j ρ j g º ¼ ν¬ A definition of the mixture pressure is introduced so that: ∇p = ∇p j − λk ∇pc
ρV = −
(
)
(4.290) (4.291)
A density correction factor, γ ρ , is also introduced.
γρ =
1
ρ
(ρ λ
kk
+ ρ jλ j
)
(4.292)
The momentum equation (4.290) can be written in the form K ρ V = − ª ∇p − γ ρ ρ g º ¼ ν¬
(4.293)
In order to simplify the MMM derivation for the species and energy equations, a diffusive phase-mass flux is introduced. This term is analogous to the diffusion mass flux in multicomponent mixtures, but refers to each phase, rather than each component in a phase. This value relates the actual mass flux of phase k to the mixture mass flux. k jk = ε sk ρ k Vk − λk ρ V (4.294) From this relation, it can be shown that the diffusive phase-mass flux is:
jk = − j j =
λk λ j ª K ∇pc + ( ρ k − ρ j ) g º ¼ ν¬
(4.295)
Chapter 4 Generalized Governing Equations: Averaging Formulations 299
This relation will be useful, and will be used henceforth. The mixture species equation for species i is obtained by adding the species i equation for phase k and phase j. j j j k k ∂ k ε sk ρ k ωk ,i + ε s j ρ j ω j ,i + ∇ ⋅ ε sk ρ k Vk ωk ,i + ε s j ρ j V j ω j ,i ∂t ′′′ = −∇ ⋅ J k ,i + J j ,i + mk ,i + m′′′,i j
(
)(
)
(4.296) From eq. (4.265), the summation of the species production over all the phases is the species production due to chemical reaction. The mixture mass fraction is:
ρωi = sk ρ k ωk ,i
k
+ s j ρ j ω j ,i
j
(4.297)
The correction factor for species advection, γ i , is introduced. j k 1 γi = λk ωk ,i + λ j ω j ,i
ωi
(
)
(4.298)
The mixture species equation is: ∂ (ερωi ) + ∇ ⋅ (γ i ρ Vωi ) ∂t (4.299) § j ω k − ω j · + m′′′ = −∇ ⋅ J k ,i + J j ,i − ∇ ⋅ ¨ k ¸ k ,i j ,i i © ¹ It should be noted that eqs. (4.282), (4.293), and (4.299) are also applicable to multiphase and not necessary to two-phase as developed here. It is important to point out that the mixture species equation still contains the species mass fraction of each phase. Therefore, the species mass fraction in phase k and phase j must be related to the mixture mass fraction through thermodynamic equilibrium. Expanding all the terms in eq. (4.297) yields:
(
)
(s ρ
k
k
+ s j ρ j ωi = sk ρ k ωk ,i
)
k
+ s j ρ j ω j ,i
j
(4.300)
Since there are two phases presented, the phase saturations add to unity. Therefore, the phase saturation of phase k can be calculated by:
sk =
ρ j ωi − ω j ,i ρ k ωi − ωk ,i
(
(
j
k
)(
)
+ ρ j ω j ,i
k
− ωi
)
(4.301)
It is important to note that when the saturation of phase k is calculated in this manner, one phase continuity equation is not solved; instead, all N species equations are solved. For more discussion on this approach to the calculation of the phase saturation, refer to the comparison of MFM and MMM models below. The last transport equation that must be examined is the energy equation. The energy equation of both phases are added together plus an unsteady and heat conduction term represents the influence of the solid matrix on the energy equation.
300 Transport Phenomena in Multiphase Systems
∂§ ¨ (1 − ε ) ρ sm hsm ∂t ©
sm
+ ε ª sk ρ k hk « ¬
k
k
+ s j ρ j hj hj
j
j
º· + »¸ ¼¹
º · = −∇ ⋅ q′′ + q′′ + q′′ + q′′′ k ,i j ,i sm ,i E »¸ ¼¹ (4.302) ′′′ The external heat generation is represented by qE , and it is the summation of the external heat generation over all of the phases. The mixture enthalpy is
k
§ ∇ ⋅ ¨ ε ª sk ρ k Vk ¬ ©«
hk
+ s j ρ j Vj
j
ρ h = sk ρ k hk
k
+ s j ρk h j
j
(4.303)
The correction factor for energy advection, γ h , is also introduced: j 1 k γ h = λk hk + λ j h j (4.304) h An assumption that corresponds with the mixture enthalpy is that each phase
(
)
is in thermodynamic equilibrium, Tk
k
= Tj
j
.
If Fourier’s law governs the
heat conduction and effective thermal conductivity can be used, the energy equation can be rewritten using an effective thermal conductivity and the diffusive phase-mass flux as: ∂ sm (1 − ε ) ρ sm hsm + ερ h + ∇ ⋅ (γ h ρ Vh ) = ∂t (4.305) § j h k − h j · + q′′′ ∇ ⋅ ( keff ∇T ) − ∇ ⋅ ¨ k k ¸E j © ¹ where the effective thermal conductivity is related to the phase saturation, porosity, geometry of the porous medium (tortuosity, ℑ ) and the conductivity of the porous media and both phases, keff = f ( k sm , kk , k j , ε , sk , ℑ) .
(
(
)
)
It is important to note that the heat generation and consumption of chemical reactions and latent heat are all embedded in this governing equation.
Comparison of MFM and MMM models
The major advantage of MFM over MMM is that MFM has the potential to capture the details of each phase, including interface. Furthermore, the nonequilibrium thermodynamics for temperature and species concentration in each phase can be obtained by MFM. Another benefit is the potential to account for the pore size distribution in a porous media and model how the saturation of each pore is not equivalent. The pore size distribution is not well researched in the literature. This is why the Leverette function of capillary pressure is commonly used, even though this function applies to soil engineering and probably is not the best choice for other applications. In the MFM model, the momentum equation can be directly inserted into the continuity equation, therefore eliminating the intrinsic velocity of each phase. It
Chapter 4 Generalized Governing Equations: Averaging Formulations 301
also can be directly inserted into the energy and species equations. Also, the phase saturation of phase j is known if the saturation of phase k is known. Similarly, only N − 1 species equations are needed in each phase because the species balance is maintained through continuity. Therefore, the main flow variables in the MFM model are sk , p j , ωk ,i
k
, ω j ,i
j
, Tk
k
, Tj
j
.
The number of equations and the number of unknowns using the MFM model for two phases when the energy equation is needed is 5 + 2 × ( N − 1) . These equations are: two continuity equations, three energy equations, and N-1 species equations for each phase. Applying the appropriate boundary conditions, this is a well posed problem. In the MMM model, applying the momentum equation directly into the continuity equation to eliminate the velocity, the main flow variables are: ρ , p, ωi , h . The variables are the mixture density, the mixture pressure, N − 1 species mass fraction, and the temperature, i.e., the number of variables is 3 + ( N − 1) . The main equations are one mixture continuity equation, N − 1 mixture mass fraction equations, and one energy equation, i.e., the number of equations = 2 + ( N − 1) . One more equation needs to be solved in order for this approach to be well posed, since there is one more dependent variable than equations. The first approach to this problem is to solve the continuity equation for one of the phases. Using the MMM model variables, the continuity equation for phase k is: ∂ (4.306) (ε sk ρ k ) + ∇ ⋅ ( λk ρ V ) = −∇ ⋅ J k + mk′′′ ∂t The other approach, which is suggested by Wang and Cheng (1996) and used in Wang and Wang (2003) to model fuel cells, is to solve species equations for all N species, and not solve the continuity equation for phase k. It is very important to note that if this approach is taken, the mixture species mass fraction as well the species mass fraction in each phase must add up to unity.
¦ ωi = ¦ ωk ,i
i =1 i =1
N
N
k
= ¦ ω j ,i
i =1
N
j
=1
(4.307)
If these criteria are met, the saturation calculated by equation (4.301) is correct. Assuming the relations in eq. (4.307) are upheld, than the species equation for Nth component can be written in terms of all the other components. N −1 ∂ N −1 ∂ ( ερ ) − ¦ (ερωi ) + ∇ ⋅ ( ρ V ) + ¦ ∇ ⋅ ( γ i ρ Vωi ) = ∂t i =1 ∂t i =1 (4.308) N −1 N −1 N −1 § j ω k − ω j · − m′′′ ¸¦i ¦ ∇ ⋅ J k ,i + J j ,i + ¦ ∇ ⋅ ¨ k k ,i j ,i © ¹ i =1 i =1 i =1 This equation is simply the continuity equation subtracted by N–1 species equations. Mathematically, this equation is not a new (independent) equation, and therefore solving for the Nth species equation still leaves the system of
(
)
302 Transport Phenomena in Multiphase Systems
equations ill-posed. This leads us to using the continuity equation for phase k to complete the system of equations. The formulation of MMM reduces the total number of equations, which, at first glance, would seem to decrease the computational time. However, when looking at the number of terms that these equations incorporate, it is not very different from the MFM. For example, examining the species equation in the MFM model for Π phases and Dx dimensions, and comparing it to the MMM model for the same number of phases and dimensions can show how many terms are solved (Table 4.2; Rice 2006a). For two-phase, two-dimensional systems, the MFM model solves a total of 10 terms, while the MMM model solves a total of 11 terms. So the computational time difference between the MMM and the MFM models are probably not as dramatic as it first appears. Since the total number of terms solved is very close, the computational time is probably limited by the convergence rate more than anything else.
Table 4.2 Comparison of MFM and MMM, Total Terms Computed in Species Equation Model MFM 1st Order Time 1st Order space (Advection) ΠDx 2nd Order Space (Diffusion) ΠDx
Π Dx
Total
Π
1
Π ( 2 Dx + 1)
MMM
(1 + Π ) Dx
Dx ( 2Π + 1) + 1
Unsaturated flow theory comes from MFM, with additional assumptions. In unsaturated flow theory (UFT), only the liquid phase is considered, and the gas pressure is assumed to be equal to the hydrostatic pressure of the gas. If the capillary pressure is assumed to be only a function of liquid saturation, and there is no mass transfer considered, then Richard’s equation is valid. Richard’s equation is: ªK º dp ∂ (4.309) (ε sk ρ k ) + ∇ ⋅ « K rk § c ∇sk + ρ k − ρ j g · » = 0 ¨ ds ¸ ∂t ¹¼ ¬ νk ©
(
)
4.7 Boltzmann Statistical Averaging
4.7.1 Boltzmann Equation
Another approach is to look at the molecular level of fluids from a statistical standpoint (Rice 2006b). From this standpoint the independent variables are space, velocity and time, while the dependent variable is a molecular distribution function for species i, fi ( x, c, t ) . The Boltzmann equation relates the
distribution function at ( x, c, t ) to the distribution function at (x + Δx, c + Δc, t + Δt ). The location in space is x, and the particle velocity is c. It is important
Chapter 4 Generalized Governing Equations: Averaging Formulations 303
to note that the particle velocity is directly related to the mass average velocity, V, that is used throughout this book. This distribution function can be related to the Navier-Stokes equations as well as other transport equations; these relationships give insight to the origin of transport coefficients such as viscosity. The species distribution function is related to the number of particles, Ni, at a location less than x0, and with a velocity less than c0 at time t. N i = particles ( x < x0 , c < c0 , t )
(4.310) ∂ 6 Ni ∂6 N = 3 3i ∂x1∂x 2 ∂x3∂c1∂c 2 ∂c3 ∂ x∂ c For simplicity, the partial derivatives with respect to all spatial directions and all velocity components are written as ∂ 3 x = ∂x1∂x 2 ∂x3 and ∂ 3c = ∂c1∂c 2 ∂c3 , respectively. The number density of particles, ni, is the average number of particles of species i per unit volume. ni ( x, t ) = ³ f i ( x, c, t ) d 3c (4.311)
fi ( x, c, t ) =
The limits of the integral are not shown because it incorporates all the velocities from −∞ to ∞ . From the distribution function, macroscopic flow variables such as density ( ρ ) , the mean species velocity (Vi), the mass mean velocity (V), the average kinetic energy (e ) and the average random kinetic energy or internal energy (e) can be defined. ρi = ³ mi fi ( x, c, t ) d 3c (4.312)
ρ V = ¦ ³ mi cf i ( x, c, t ) d c
3
i
ρi Vi = ³ mi cfi ( x, c, t ) d 3c
(4.313) (4.314)
1 2 3 (4.315) ³ mi c fi ( x, c, t )d c i2 1 ρ e = ¦ ³ mi w i2 fi ( x, c, t )d 3c (4.316) i2 The mass of a single particle is mi. The random velocity ( w i ) for species i is defined as the difference between a molecule’s actual velocity and the mass mean velocity. wi = c − V (4.317) The temperature is a function of the average random kinetic energy. 3 ρ k BT = e (4.318) n 2 It is important to note that the average kinetic energy and average random kinetic energy neglect both rotational and vibrational kinetic energy and are only accurate in cases where these energies are small compared to translational energy, as is the case for a monatomic ideal gas. Since this description is the most accurate for a monatomic gas, the terms molecule and particle are
ρe = ¦
304 Transport Phenomena in Multiphase Systems
interchangeable. Also, for the following discussion, all the gas particles are assumed to have a constant mass. As is discussed above, the Boltzmann equation relates the distribution function at (x, c, t) to the distribution function at a projected position (x+ x, c+ c, t+ t). The change in spatial position and velocity are related to the particle velocity and particle acceleration (a), which is the external body force exerted on a particle such as gravity or a Lorentz force. Δx = cΔt (4.319) (4.320) Δc = aΔt The Boltzmann equation for species i is (Chapman and Cowling, 1970): ( fi ( x + cΔt , c + aΔt , t + Δt ) − fi ( x, c, t ) ) d 3xd 3c = Ωi ( fi ) d 3xd 3cdt (4.321) The collision term, the rate of change of fi due to collisions, is denoted by i(fi). This function is discussed in more detail below. Expanding the first term by Taylor series yields: § ∂f · fi ( x + cΔt , c + aΔt , t + Δt ) = f i ( x, c, t ) + ¨ i + ∇fi ⋅ c + ∇c fi ⋅ a ¸ Δt (4.322) © ∂t ¹ The gradient operators without a subscript and with a subscript c are ª∂ ª∂ ∂ ∂º ∂ ∂º ∇=« » and ∇c = « » , respectively. Therefore, ¬ ∂x1 ∂x 2 ∂x3 ¼ ¬ ∂c1 ∂c 2 ∂c3 ¼ the Boltzmann equation can be rewritten as: Df i ∂fi = + c ⋅ ∇f i + a ⋅ ∇ c f i = Ω i ( f ) (4.323) Dt ∂t If no collisions occurs, the Boltzmann equation states that the distribution function does not change as it moves along the trajectory of a particle, Dfi / Dt = 0 . In most real systems, however, collisions do occur. Particles exhibit long-range and short-range forces. When the long-range forces are negligible, as is the case for a neutral gas, the particles are assumed to interact only when the short-range forces are prevalent. The short-range forces are prevalent when the distance between the particles is on the order of the particle diameter. Since the particles are this close to one another only for a very short time, the interaction is treated as a collision. Also, it is assumed that only binary collisions exist, therefore the higher order collisions, Ωi ( fi fi ), Ω( f i f i f i ), Ω( fi fi fi fi )… are not considered. This assumption is not good for dense gases or fluids. The collision between two particles with velocities c and c A before collision results in the particles having velocities c′ and c′A after collision. The particles are assumed to retain their mass before and after collision, therefore the conservation of momentum and energy require: mc + m Ac A = mc′ + m Ac′A (4.324) 121 1 1 2 mc + mAc 2 = mc′2 + mAc′A (4.325) A 2 2 2 2
Chapter 4 Generalized Governing Equations: Averaging Formulations 305
Line of Impact
c′A
b
cA
c
φ
θ
c′
c − cA
c′ − c′A
Figure 4.14 Elastic Collision of Two Particles.
Collision Plane
Also, the collisions are considered elastic, therefore: c − c A = c′ − c′A
(4.326) Squaring eqs, (4.323) and (4.325) and adding them together yields the energy equation, (4.324). The impact of two particles is presented in Fig. 4.14. Since the collisions are frictionless, momentum and energy are only transferred in the normal direction (along the line of impact). Also, the impact occurs on a single plane; therefore, the relative velocity of the two particles can be deflected by an angle θ while the magnitude of the relative velocity is maintained. The deflection angle θ , is only a function of the impact parameter b, when the intermolecular forces are given. The impact parameter is the distance of closest approach of the center of mass of the two particles, if the particle’s trajectory was not deflected by the collision. Since we are considering many particles, and want to know how they collide with each other, the differential cross section ϑ and the solid angle d Ω are defined to account for all possible angles the plane of interaction may be as well as the deflection angle. d Ω = sin (θ ) dθ dφ (4.327) ϑ d Ω = bdbdφ (4.328) The azimuthal angle is φ , which defines the plane in which a collision occurs. The collision function, Ωi ( f i ) , for only binary collisions is (Chapman and Cowling, 1970): ′ Ωi ( f i ) = ³³ ( f i′f A − fi f A ) c − c A ϑ d Ωd 3c A (4.329)
306 Transport Phenomena in Multiphase Systems
′ The distribution functions after a collision are f ′ and f A , and those before the collisions are f and f A . Since the collision function comes from the impact of particles that conserve mass, momentum and energy, the integration of the elementary collision invariants and the collision function are zero (Cercignani, 1975). 3 (4.330) ³ψ Ωi ( f ) d c = 0
The elementary collision invariants are ª º « mi » « » (4.331) ψ = « mi c » «1 » 2 « mi c » ¬2 ¼ The velocity, c, has three components. To gain a general transport equation, the Boltzmann equation can be multiplied by a transport variable, ψ i , and
integrating over d 3c yields: ª ∂fi º3 3 (4.332) ³ψ « ∂t + c ⋅ ∇fi + a ⋅ ∇c fi » d c = ³ψ Ωi ( f ) d c ¬ ¼ If the transport variable is an invariant, than the right hand side of eq. (4.332) is zero. The first term in the derivative can be rewritten using the chain rule and Liebnitz integral rule. ∂f i 3 ∂ψ i 3 ∂ 3 (4.333) ³ψ ∂t d c = ∂t ³ψ fi d c − ³ fi ∂t d c The second term can be rewritten as: ∂fi 3 ∂ ∂ψ 3 3 (4.334) ³ψ c j ∂x d c = ∂x ³ψ c j fi d c − ³ fi c j ∂x d c j j j Note that the subscript i refers to a specific species and the subscript j refers to a direction. The independent variable is time, t, and x and c are both functions of time [ x = x(t ) and c = c(t ) ]; therefore, ∂c / ∂x = 0 . Also, the acceleration is assumed to be independent of the particle velocity, which is valid when gravitational or electromagnetic forces are involved. ∂fi 3 ∂ ∂ψ 3 3 (4.335) ³ψ a j ∂c d c = ∂c ³ψ a j fi d c − ³ fi a j ∂c d c j j j For simplicity, any function K multiplied by fi and integrated is written as: 3 (4.336) ³ Kfi d c = K Therefore, the Boltzmann equation can be rewritten as a general transport equation. This equation relates the microscale physics to macroscale physics. ∂ ∂ψ (4.337) ψ− + ∇ ⋅ ψ c − c∇ψ + ∇c ⋅ ψ a − a∇cψ i = ψ Ωi ∂t ∂t
()
()
Chapter 4 Generalized Governing Equations: Averaging Formulations 307
To retain species continuity, ψ = mi , and the mass of a particle does not change ∂mi ∂mi ∂mi = = = 0 . Therefore, with time, location or velocity, ∂t ∂x ∂c ∂ρ ∂ mi + ∇ ⋅ mi c = i + ∇ ⋅ ( ρi Vi ) = 0 (4.338) ∂t ∂t To retain the momentum equation, the invariant may be set to the particle momentum, ψ = mi c j .
()
∂mi c j ∂ mi c j − + ∇ ⋅ mi c j c = 0 ∂t ∂t The fourth term drops out because c is not a function of x. of the general transport equation drop out, because ∇c ⋅ mi c j a − a∇c mi c j = mi ( a − a ) = 0
(
)
(4.339) The last two terms (4.340)
(
)
The body force term comes from the second term in eq. (4.339), ∂c j = mi a j . The third term in eq. (4.339) may be the term of because mi ∂t greatest interest, because it contains information about the microscopic and the bulk movements of particles.
mi c j c = V j + w i , j
(
)(V + w ) = ρ V V + m w
i i j i
i, j w i
(4.341) Now
The random velocity of species i in the j direction is represented by w i , j .
the momentum equation can be rewritten in terms of the average velocity. ª ∂ρi V j º (4.342) ¦ « ∂t + ∇ ⋅ ( ρi V j V ) = ρi X j − ∇ ⋅ ρi w i , j w i » i¬ ¼
(
)
Since this momentum equation is identical to the Navier-Stokes equation, the stress tensor is: (4.343) τ = − ¦ ρi w i w i
i
The scalar gas pressure is defined as one-third of the trace of the stress tensor. 1 p = −¦ ρi w i2,1 + w i2,2 + w i2,3 (4.344) i3 Pressure always has a positive value, as it does in this definition. To retain the remainder of the stress tensor, τ ′ , the pressure can be subtracted from the total stress tensor. 1 τ ′ = ¦ − ρ i w i w i + ρi w i2,1 + w i2,2 + w i2,3 I 3 i
(
)
(
)
·§ ·º p ¸ + ¨ ¦ ρi w i2,3 − p ¸ » I ¹ ©i ¹¼ (4.345) Reducing the stress tensor to be governed by Newton’s law of viscosity, ·§ p ¸ + ¨ ¦ ρi w i2,2 − ¹ ©i
ª º 1 ª§ = − « ¦ ρi w i w i − pI » + «¨ ¦ ρi w i2,1 − ¬i ¼ 3 ¬© i
308 Transport Phenomena in Multiphase Systems
¦ ρi w i w i − pI = − μ ( ∇V + ( ∇V )
i
T
)
(4.346)
Now the total stress tensor can be written as: (4.347) 3 Summing the momentum equation over all the species, applying continuity and substituting the stress tensor, the momentum equation can be written as: DV ρ = ∇ ⋅ (τ ) + ρ X (4.348) Dt For the energy equation, the invariant may be set to the kinetic energy of a 1 particle, ψ = mi c 2 . 2 ª º 1 ∂ mi c 2 «∂ 1 §1 ·» (4.349) ¦ « ∂t 2 mi c2 − 2∂t + ∇ ⋅ ¨ 2 mi c2c ¸ » = 0 i« © ¹» « » ¬ ¼ 2 The fourth term drops out because c is not a function of xj. The last two terms of the general transport equation drop out, because §1 · 1 ∇c ⋅ ¨ mi c 2a ¸ − a∇c mi c2 = mi c ⋅ ( a − a ) = 0 (4.350) 2 ©2 ¹
i
2 T τ = −¦ ρi w i w i = − pI + μ ª∇V + ( ∇V ) º − μ ( ∇ ⋅ V ) I
¬ ¼
The energy added due to the body force comes from the second term in eq. 1 ∂c 2 (4.349), because mi = mi c ⋅ a = ρi Vi ⋅ X . The first term can be rewritten in 2 ∂t terms of the average velocity and the deviant velocity. ∂ §1 ∂§ 1 ∂§ § V2 · · 2· 2 2· (4.351) ¨ mi c ¸ = ¨ ¦ ρi w i + V ¸ = ¨ ρ ¨ e + ¸¸ ¨ ¸ ∂t © 2 ∂t © i 2 2 ¹¹ ¹ ∂t © © ¹
(
)
The third term in equation (4.349) can also be expanded, resulting in several terms, such as advection, heat flux, and viscous dissipation terms. §1 · §1 · ¦ ∇ ⋅ ¨ 2 mi c2c ¸ = ¦ ∇ ⋅ ¨ 2 ρi w i2 + 2w i V + V 2 ( V + w i ) ¸ © ¹ i © ¹i
(
)
§ §1 §1 · ·· = ¦ ∇ ⋅ ¨ ρi w i2 + V 2 V ¸ + ¦ ∇ ⋅ ¨ ρi ¨ w i2 w i + w i ( w i ⋅ V ) ¸ ¸ ©2 ¹i ¹¹ i © ©2
(
)
(4.352)
§§ § §1 2 V2 · · ·· = ∇⋅¨ ρ ¨e + ¸ V ¸ + ¦ ∇ ⋅ ¨ ρi ¨ w i w i + w i ( w i ⋅ V ) ¸ ¸ ¨ ¸i 2¹ ¹ ¹¹ © ©2 ©© The last term includes the mechanical work of the stress tensor, as well as the heat flux components. 1 (4.353) ¦ 2 ρi w i2 w i = q′′ i
Chapter 4 Generalized Governing Equations: Averaging Formulations 309
¦ ρi w i ( w i ⋅ V ) = τ ⋅ V
i
(4.354)
Therefore, the energy equation can be rewritten in the standard form. § § ∂§ § V2 · · V2 · · ¨ ρ ¨e + ¸ + ∇ ⋅ ¨ ρV ¨ e + ¸¸ ¸ ¸ = −∇ ⋅ q′′ + ∇ ⋅ (τ ⋅ V ) + ρ V ⋅ X (4.355) ¨ ¨ ¸ ∂t © © 2 ¹¹ 2 ¹¹ © © The bulk transport effects are derived from a statistical molecular model, to get equations governing conservation of mass, momentum and energy. Since only binary collisions are considered, and the rotational and vibrational energies of a molecule are neglected, the model is valid for dilute monatomic gases. Another insight to the Boltzmann equation is looking at the equilibrium distribution function of a gas mixture at rest, known as the Maxwell-Boltzmann distribution. At equilibrium, the rate of change of the distribution function is zero: Dfi =0 (4.356) Dt In order for the rate of change of the distribution function to be zero, the collision function must also be zero, Ωi = 0 . Therefore: ′ fi′f A = fi f A (4.357) The solution to this function is the Maxwell-Boltzmann distribution function: 3/ 2 ª mi ( c − V )2 º § mi · » (4.358) fi ( x, c, t ) = ni ¨ ¸ exp « − 2k B T » « © 2π k BT ¹ ¬ ¼ This function describes how the particle velocities vary at each location. Modeling of the Boltzmann equation can be done using the Lattice Boltzmann model (McNamara and Zanetti, 1988). This model discretizes the
c3
c4 c5
(a)
c2 c0 c6
Δxα = cα Δt
(b)
c1
Figure 4.15 Lattice Boltzmann model discretization and lattice structure for the (a) velocity vectors corresponding to f i of a node, (b) Lattice structure connecting nodes.
310 Transport Phenomena in Multiphase Systems
Boltzmann equation into six velocity components (seven if a resting particle is included) at different locations on a mesh. A representation of the Lattice Boltzmann model is presented in Fig. 4.15. The mesh is constructed so that the grid-spacing of the nodes is equivalent to a velocity component multiplied by the time-step. The distribution function is broken into two parts: a streaming step and a collision step. The particles will stream along the grid “lattices” during the streaming step. During the collision step, a collision function is applied at the node. This collision function usually incorporates mass and momentum conservation, but generally the collisions should conserve mass, momentum, and energy.
4.7.2 Lattice Boltzmann Model (LBM)
Most numerical codes are based on the Navier-Stokes equations, which treats a fluid as a continuous field. It is well known that a fluid is made of a discrete number of particles or molecules. Since the number of molecules is extremely large (Avogadro’s number = 6.022×1023 atoms/mole) for almost all practically sized systems, it may never be computationally viable to track each particle and its interactions with other particles. The number of molecules in a given region and the molecular interaction are described through the fluid’s density and transport coefficients (i.e., viscosity) in the continuous model. Modeling the individual molecules for a small system over a small period of time has been achieved by molecular dynamic simulations (MDS). The computational requirements needed in these simulations can be greatly reduced if the degrees of freedom of the system are reduced. Also, instead of considering individual molecules, groups of molecules can be considered. The degrees of freedom can be reduced by restricting the movement of the molecules to a lattice. A lattice is simply a predefined direction in which a molecule can move. The first of such models is the lattice gas automata (LGA; Hardy et al., 1973). The evolution equation of the LGA is: nα ( x + Δxα , t + 1) − nα ( x, t ) = Ωα ( n ( x, t ) ) (4.359) (α = 0,1" , M ) Each node has a Boolean variable, nα ( x, t ) , for each connecting lattice; if there is a molecule present it is unity, and if none is present it is zero. The number of connecting lattices is M, and the collision operator, Ωα , is dictated by collision rules. On each node, only one particle can be found moving in a given direction at a given time. This condition is imposed for memory efficiency and leads to a Fermi-Dirac local equilibrium distribution (Frisch et al. 1987). The two subsequent steps of the LGA simulations are streaming and collision. In the first step, a particle can move in one lattice site corresponding to the link direction. In the second step, collision rules are applied on sites where two or more particles meet. The rules are: the particles change configuration when two or more particles meet. Because of the Boolean variables, the LGA method is quite noisy. The Lattice Boltzmann Method (LBM) is introduced by modeling
Chapter 4 Generalized Governing Equations: Averaging Formulations 311
the LGA equation by single particle distribution functions, fi ( x, t ) , instead of a Boolean variable (McNamara and Zanetti, 1988). The Lattice Boltzmann Equation (LBE) is: fα ( r + cα Δt , t + Δt ) − fα ( r , t ) = Ωα ( f ( x, t ) ) Δt (α = 0,1,", M ) (4.360) A triangular lattice structure should be used because it has sufficient symmetry to make a solution independent of how the lattice is laid out on the geometry (Frisch et al., 1986). The velocity vectors surrounding a node as well as the lattice structure are presented in Figure 4.15. The velocity vector, c0 , symbolizes a particle at rest. The discretization of the LGA and the LBE is the same when Δxα / Δt = cα , therefore the lattice lengths Δxα are equal to cα Δt . The LBE is a discrete velocity model of the Boltzmann equation. In the LBE, the macroscopic flow variables, density and velocity, are defined as:
ρ = ¦ fα
α
(4.361) (4.362)
ρ V = ¦ fα cα
α
The collision operator can also be discretized. The collision operator for the th component of the distribution, Ωα , contains a summation over all the velocity components, as well as all the collision cross sections, ϑβα , that the th component of the distribution function makes with the th component. For a more in-depth description of the collision operator, refer to the discussion on the Boltzmann equation. The scattering angle is predetermined by the discretization of f.
Ωα ( f ) = ¦ ( fα′ f β′ − fα f β )ϑβα cα − c β
β
(4.363)
The collision operator can be linearized, and the dimension of the matrix can be greatly reduced by expanding the distribution function into an equilibrium, fα , and non-equilibrium, fα* , component. fα = fα + fα* (4.364) The equilibrium distribution function depends on the local macroscopic variables and should satisfy the following constraints:
¦ fα = ρ
α
(4.365) (4.366)
¦ fα cα = ρ V
α
If fα* fα , then the collision operator can be linearized. Ωα ( f ) = Ωα fα + ¦ Ωαβ fα* + O fα*2
β
()
()
()
(4.367)
312 Transport Phenomena in Multiphase Systems
The linearized collision matrix is defined as: Ωαβ ≡
∂Ωα fα ∂f β
( ).
Also, the
collision operator of an equilibrium distribution is zero, Ωα fα = 0 , the LBE is: fα ( r + cα Δt , t + Δt ) − fα ( r, t ) = ¦ Ωαβ f β − f β Δt
β
()
(
)
(α = 0,1,", M )
(4.368) Since the linearized collision matrix must conserve mass and momentum, it must satisfy the following constraints:
¦ Ωαβ
β
b
=0 =0
(4.369) (4.370)
¦ Ωαβ c β
A further assumption of the linearized collision matrix is that the local particle distribution relaxes to an equilibrium state at a single rate, τ . 1 Ωαβ = − δαβ (4.371) The delta function is unity when α = β , and zero when α ≠ β . This approximation is called the BGK collision term (Bhatnagar et al., 1954). With a single rate relaxation period, the Lattice Boltzmann BGK (LBGK) equation is: 1 fα ( r + cα Δt , t + Δt ) − fα ( r, t ) = − fα − fα Δt (4.372)
τ
τ
(
)
Now, the relaxation period, as well as the equilibrium distribution function must be specified. The equilibrium distribution function can be found by taking a second-order expansion of the local velocity of the Maxwell-Boltzmann distribution function (Házi et al., 2002). 2 (4.373) fα = wα ª A + Bcα ⋅ V + C ( cα ⋅ V ) + DV 2 º ¬ ¼ For a two dimensional triangular grid as shown in Fig. 4.15, the constants are: 1 (4.374) wα = ρ , A = d 0 , B = C = 0, D = 2 for α = 0 c (1 − d0 ) 8 1 D ,B= 2, C= , D = − 2 for α = 1,...,6 (4.375) wα = ρ , A = 4 6 6c 12c 6c 2 where c is:
¦ ( cα ⋅ x1 )( cα ⋅ x1 ) = 3c 2
α
(4.376)
where x1 is Cartesian directions, and α is a single component of a lattice vector. From these expressions, the macroscopic momentum equations can be written. The pressure and kinematic viscosity are:
Chapter 4 Generalized Governing Equations: Averaging Formulations 313
(4.377) 2 c2 § 1· ν = ¨τ − ¸ (4.378) 4© 2¹ Therefore the relaxation time is related to the kinematic viscosity, and the constant d 0 is related to the pressure. The Lattice Boltzmann equations are solved for the distribution function fα at every node. It is important to note that the Navier-Stokes equations are matched with these criteria only if the flow has a nearly constant density. Also, note that the above derivations apply only for isothermal systems, because the conservation of energy was never considered. A brief review of models that have attempted to capture thermal effects is discussed in Házi et al. (2002). Even though these models can theoretically capture arbitrary Prandtl number fluids, they are strictly limited by stability constraints. There are two major forms of boundary conditions for hydrodynamic problems: wall boundary conditions and flow boundary conditions, in which flow enters or exits a computational domain. The no-slip wall boundary condition in LBM can be implemented by a bounce-back method. In this method, the distribution functions that arrive at a boundary are bounced back with the same distribution function in which they arrived. In a flow boundary condition, the distribution function needs to be specified during the streaming phase. This specification is driven by the macroscopic velocity or pressure. The closure to this open set of equations was achieved through the work of Zou and He (1997).
p=
ρ
(1 − d0 ) c 2
4.7.3 LBM for Multiphase Flows
One of the most important aspects of modeling multiphase flow is capturing interfacial dynamic effects such as surface tension. The first group of models is based on the color model of Gunstensen et al. (1991). In this model, red and blue particle distribution functions ( f R and f B , respectively) were introduced to represent two different fluids. The total distribution function is the sum of the partial distribution functions: fα = fαR + fαB (4.379) The LBE can be written for each phase, k. k fαk ( x + ci Δt , t + Δt ) − fαk ( x, t ) = Ωα ( x, t ) Δt (4.380) The collision operator can be split into two components: the collisions that lead to the local equilibrium similar to the LBGK model and the collisions that contribute to the dynamics of the interfaces between the two different particles. k k k Ωα = Ωα,1 + Ωα,2 (4.381) The first collision operator is the same as the BGK operator, but the relaxation parameter only pertains to phase k.
314 Transport Phenomena in Multiphase Systems
k Ωα,1Δt = −
1
τk
(f
α
k
− fαk
)
(4.382)
The density of each phase is:
ρ k = ¦ fαk
α
(4.383)
The total mixture density, ρ , is the summation of the phase densities,
ρ = ¦ ρ k . Similarly, the phase velocity is:
k
ρ k Vk = ¦ fαk cα
α
(4.384) is the summation of the phase
The
total
mixture
k
velocity,
V,
velocities, ρ V = ¦ ρ k Vk .
The second collision operator, which contributes to
the dynamics in the interfaces and generates surface tension, is ª ( e ⋅ F )2 1 º A k Ωα,2 Δt = k F « α 2 − » 2 2» «F ¬ ¼ where F is the local color gradient, defined as:
(4.385)
F ( x ) = ¦ eα ª ρ r ( x + eα ) − ρb ( x + eα ) º ¬ ¼
α
(4.386)
Note that F will vanish in a single phase fluid where the density is constant. The direction and length of each lattice link is represented by the vector eα . The parameter Ak determines the surface tension. Defining a second collision operator alone does not keep different phases separated. To keep the phases separated, the local color momentum j, must be aligned with the direction of the local color gradient after collision.
j = ¦ fαr − fαb eα
α
(
)
(4.387) (4.388)
jF = jF
One of the major drawbacks of the color model is that it does not have any thermodynamic background; therefore it is limited to hydrodynamic studies. Another two-phase flow model was originated by Shan and Chen (1993), in which the surface interaction can be maintained automatically. The second collision term is different from the color model. It is: k Ωα,2 Δt = eα ⋅ F k (4.389) th Here, the effective force on the k phase comes from a pairwise interaction between the different phases.
Chapter 4 Generalized Governing Equations: Averaging Formulations 315
F k ( x ) = −¦¦Vkk ′ ( x, x + eα ) eα
k′ α
(4.390)
The interaction potential between the phases is Vkk ′ , and it is defined as:
Vkk ′ = Gkk ′ ( x, x′ )ψ k ( x )ψ k ′ ( x′ )
(4.391)
The strength of the interaction is Gkk ′ and ψ k ( x ) is a density function for phase k at x. A discussion of the form of ψ is included in Chen and Doolen (1998) and gives a nonideal equation of state, which is needed for a two-phase flow model. The strength of the interaction only accounts for the nearest neighbor on a lattice link, x − x′ = ei , and is zero for any other combination. The surface tension, σ Av , is equal to: M Gkk ′ σ Av ≈ (4.392) 2 D ( D + 1) The constant D is the constant in front of the V2 term in the equilibrium distribution function. The order-parameter, M , is 1q2 2 M ≡ ª( ρ − ρ ) º (4.393) « » ρ¬ ¼ The rounded overbar represents the spatial average of a value over the entire k lattice. It should be noted that the collision operator, Ωα,2 , for this model does not satisfy local momentum conservation. Another approach to solving multiphase flow problems using the LBM method is a technique that uses time splitting (Premnath et al., 2005). In this method, the interface is directly tracked using a kinematic equation for the level set function, φ . The level set function varies between -1 and 1, and the interface corresponds to the location where the level set function equals 0. ∂φ + V ⋅ ∇φ = 0 (4.394) ∂t The first collision term is the BGK collision term. The second collision term is: 1· § (4.395) Ω 2 = ¨ 1 − ¸ S i ( x, t ) © 2τ ¹ The effect of force interactions is include in the source term, Sα . This term includes the mean-field force, F, as well as an external body force, Fext, representing gravity or a Lorentz force. ( c − u ) ⋅ ( F + Fext ) (4.396) Sα = α fα ρ cs2 The mean-field force can be broken into two components, the force from the nonideal part of the equation of state, ψ , and the surface tension force, Fs . F = −∇ψ + Fs (4.397)
1
316 Transport Phenomena in Multiphase Systems
The nonideal part of the equation of state was modeled by the Carnahan-Starlingvan der Waals equation of state. The surface tension force is related to the density and its gradient by: Fs = κρ∇∇ 2 ρ (4.398) where the surface tension parameter, κ , is related to the surface tension, σ Ag , by:
§ ∂ρ · (4.399) σ A g = κ ³ ¨ ¸ dn © ∂n ¹ where n is the normal direction of the interface. Finally, the macroscopic mass and momentum balance can be incorporated into the distribution function in nodes near the interface by writing the distribution function in terms of the gradient of the momentum flux. This is done by expanding the distribution function around its local equilibrium in terms of the Knudson number and applying the Chapman-Enskog analysis to the Taylor series expansion of the LBM equation. In this method, the interface is directly tracked; therefore, there are nodes that lie near the interface. If a node lies near an interface, the distribution function for the phase that exists on the node is known. However, the distribution on the nodes in the other phase is not known, but is calculated through the interfacial mass and momentum balances; these balances are related to the distribution function from the results of the Chapman-Enskog analysis. A final approach to multiphase flow modeling using the LBM is the free energy approach, in which the equilibrium distribution function is defined consistently with thermodynamics. This model is introduced by Swift et al. (1995, 1996). The models mentioned here do not take into account thermal transport. Other models that capture heat transfer of two-phase systems are reviewed by Házi et al. (2002) and include modeling the energy transfer with an added distribution function. Recently work has been done to incorporate solid/liquid/vapor interaction with the LBM method. Latva-Kokko and Rothman (2005) used the LBM to simulate the capillary rise between two horizontal plates as well as in a rectangular tube. Their model used the color-gradient LBM to distinguish the phases and showed that it gave appropriate static contact angles for both imbibition and drainage. The development of the LBM method has the potential to be very useful in computational flow modeling of multiphase systems, however, some issues first need to be resolved. The first issue is that the models that conserve energy are not stable across a wide variety of problems. Since many multiphase systems utilize the latent heat associated with phase change, the LBM method must have the capability of conserving energy to make it a useful model for practical applications.
2
Chapter 4 Generalized Governing Equations: Averaging Formulations 317
References
Alazmi, B. and Vafai, K., 2000, “Analysis of Variants within the Porous Media Transport Models,” ASME Journal of Heat Transfer, Vol. 122, pp. 303-326. Alazmi, B. and Vafai, K., 2001, “Analysis of Fluid Flow and Heat Transfer Interfacial Conditions between a Porous Medium and a Fluid Layer,” International Journal of Heat and Mass Transfer, Vol. 44, pp. 1735-1749. Basu, S., 2005, Personal Communication, University of Connecticut, Storrs,
CT.
Bhatnagar, P.L., Gross, E.P., and Krook, M., 1954, “A Model for Collision Processes in Gases. I Small Amplitude Processes in Charged and Neutral OneComponent Systems,” Physical Review, Vol. 94, pp. 511-525. Boysan, F., 1990, A Two-Fluid Model for Fluent, Flow Simulation Consultants, Ltd., Sheffield, England. Cercignani, C., 1975, Theory and Application of the Boltzmann Equation, Scottish Academic Press. Chapman, S. and Cowling, T.G., 1970, The Mathematical Theory of Non-uniform Gases, Cambridge University Press, Cambridge, UK. Chen, S. and Doolen, G.D., 1998, “Lattice Boltzmann Method for Fluid Flows,” Annual Review of Fluid Mechanics, Vol. 30, pp. 329-364. Duval, F., Fichot, F., and Quintard, M., 2004, “A Local Thermal NonEquilibrium Model for Two-Phase Flows with Phase-Change in Porous Media,” International Journal of Heat and Mass Transfer, Vol. 47, pp. 613-639. Faghri, A., 1995, Heat Pipe Science and Technology, Taylor & Francis, New York. Frisch, U, d’Humiéres, D., Hasslacher, B., Lallemand, P., Pomeau, Y. and Rivet, P.P.,1987, “Lattice Gas Hydrodynamics in Two and Three Dimensions,” Complex System, Vol. 1, pp. 649-707. Frisch, U., Hasslacher, B., and Pomeau, Y., 1986, “Lattice-Gas Automata for the Navier-Stokes Equations,” Physical Review Letters, Vol. 56, pp. 1505-1508. Gray, W.G., 2000, “Macroscale Equilibrium Conditions for Two-Phase Flow in Porous Media,” International Journal of Multiphase Flow, Vol. 26, pp. 467-501. Gunstensen, A.K., Rothman, D.H., Zaleski, S., and Zanetti, G., 1991, “Lattice Boltzmann Model of Immiscible Fluids,” Physical Review A., Vol. 43, pp. 43204327. Hadley, G.R., 1986, “Thermal Conductivity of Packed Metal Powders,” International Journal of Heat and Mass Transfer, Vol. 29, pp. 909-202.
318 Transport Phenomena in Multiphase Systems
Hardy, J., Pomeau, Y., and Pazzis, O., 1973, “Time Evolution of a TwoDimensional Model System. I. Invariant States and Time Correlation Functions,” Journal of Mathematical Physics, Vol. 14, pp. 1746-1759. Házi, G., Imre, A.R., Mayer, G., and Farkas, I., 2002, “Lattic Boltzmann methods for two-phase flow modeling,” Annals of Nuclear Energy, Vol., 29, pp. 14211453. Hewitt, G. F., 1998, “Multiphase Fluid Flow and Pressure Drop,” Heat Exchanger Design Handbook, Vol. 2, Begell House, New York. Ishii, M., 1975, Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles. Kaviany, M., 1995, Principles of Heat Transfer in Porous Media, 2nd ed., Springer Verlag, New York. Latva-Kokko, M. and Rothman, D.H. 2005, “Static Contact Angle in Lattice Boltzmann Models of Immiscible Fluids,” Physical Review E, Vol. 72, pp. 046701(1-7). Leverette, M.C., 1940, “Capillary Behavior in Porous Solids,” Transaction of AIME, Vol. 142, pp. 152-169. McNamara, G.R. and Zanetti, G., 1988, “Use of the Boltzmann Equation to Simulate Lattice-Gas Automata,” Physical Review Letters, Vol. 61, pp. 23322335. Nam, J.-H., and Kaviany, M., 2003, “Effective Mass Diffusivity and Water Saturation Distribution in Single- and Two-Layer PEMFC Diffusion Medium,” International Journal of Heat and Mass Transfer, Vol. 46, pp. 4595-4611. Nield, D.A., and Bejan, A., 1999, Convection in Porous Media, 2nd ed., SpringerVerlag, New York. Premnath, K.N., Nave, J-C., and Banerjee, S., 2005, “Computation of Multiphase Flows with Lattice Boltzmann Methods,” Proc. 2005 ASME IMECE, Nov. 5-11, Orlando, FL, USA. Rice, J. 2006a, Personal Communication, University of Connecticut, Storrs, CT. Rice, J. 2006b, Personal Communication, University of Connecticut, Storrs, CT. Shan, X. and Chen, H., 1993, “Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components,” Physical Review E., Vol. 47, pp. 18151819. Swift, M.R., Osborn, W.R. and Yeoman, J.M, 1995, “Lattice Boltzmann simulation of nonideal fluids,” Physical Review Letters, Vol. 75, pp. 830-833. Swift, M.R., Orlandini, S.E., Osborn, W.R. and Yeoman, J.M, 1996, “Lattice Boltzmann simulations of liquid-gas and binary-fluid systems,” Physical Review Letters E, Vol. 54, pp.5041-5052.
Chapter 4 Generalized Governing Equations: Averaging Formulations 319
Wang, C.Y. and Cheng, P., 1996, “A Multiphase Mixture Model for Multiphase Multicomponent Transport in Capillary Porous Media Part I: Model Development,” International Journal of Heat and Mass Transfer, Vol. 39, pp. 3607-3618. Wang, Z.H. and Wang, C.Y. 2003, “Mathematical Modeling of Liquid-Feed Direct Methanol Feed Cells,” Journal of The Electrochemical Society, Vol. 140, pp. A508-A519. Ward, J. C., 1964, “Turbulent Flow in Porous Media,” ASCE J. Hyd. Div., Vol. 90, HY5, pp. 1-12. Whitaker, S., 1977, “Simultaneous Heat, Mass and Momentum Transfer in Porous Media: a Theory of Drying,” Advances in Heat Transfer, Vol. 13, pp. 119-203. Zou, Q.S. and He, X.Y., 1997, “On pressure and velocity boundary conditions for the lattice Boltzmann BGK model,” Physics of Fluids, Vol. 9, pp. 1591-1598.
Problems
4.1. The volume fraction of phase k in a multiphase mixture containing Π phases is αk. Supposing the density for phase k is ρk, what is the mass fraction of phase k? What is the total density of the multiphase mixture? 4.2. A saturated mixture of water and its vapor at 180 °C flows in a 0.1 m ID tube with a mass flow rate of 0.5 kg/s. The liquid water is dispersed in the vapor phase in the form of 0.1- mm -diameter droplets, and the quality of the mixture is x=0.75. The average velocity of the vapor phase is 30m/s. Find: (a) the average velocity of the liquid droplets, and (b) the interactive force between the liquid and vapor phase. 4.3. If the liquid droplet is subcooled at a temperature of 175 °C, estimate the interphase heat transfer rate between liquid and vapor phases in the above problem. 4.4. The mass flux in the channel is defined as the mass flow rate per unit area. What is the total mass flow rate for the multiphase flow discussed in Example 4.1? 4.5. A capillary porous structure with partial heating and evaporation on its upper surface is shown in Fig. P4.1. The entire porous structure is saturated with liquid from the bottom surface (y=0), which is connected to a pool for ′′ liquid supply. A constant heat flux, q0 , is applied over part of the upper surface (0<x<Lxf), which is impermeable to the fluid. The rest of the upper surface (Lxf<x<Lx) exposes liquid in the pores of the capillary structure to the vapor space above. Heat is transferred from the left portion of the upper
320 Transport Phenomena in Multiphase Systems
surface to the right portion of the upper surface, where evaporation takes place. The left- and right-side walls are adiabatic and impermeable. At steady-state, the liquid is drawn from the bottom of the porous structure and flows to the right portion of the upper surface due to evaporation. If the liquid velocity at the bottom surface and the vapor velocity at the upper surface are uniform, determine the liquid velocity at the bottom, vA , and the vapor velocity at the upper surface, vv .
Figure P4.1
4.6. The area-averaged velocity components in the x- and y-directions in Problem 4.5 can be estimated by Darcy’s law: K ∂p K ∂p u=− v=− μ ∂x μ ∂y where K is permeability. The continuity equation is ∂u ∂v + =0 ∂x ∂y Determine the analytical solution for the pressure and velocity distributions in the porous structure. 4.7. Specify the energy equation and the corresponding boundary conditions for Problem 4.5; then nondimensionalize them. Write a computer program to obtain dimensionless temperature distribution in the porous media. 4.8. For a porous bed made of spherical balls packed in a body-centered-cubic (BCC) formation, use the continuity equation (3.41) and momentum equations (3.61) – (3.63), to develop the boundary conditions and numerically solve the flow field, and determine the permeability of the
Chapter 4 Generalized Governing Equations: Averaging Formulations 321
packed bed by plotting the product of v vs. dp/dx. Assume that the elemental cube’s sides all have a length of 10-3 m and that the spheres all have a radius of 4.6×10-4 m. (Note: the spheres are slightly compressed against each other, and therefore their radii are slightly larger than if they were only contacting one another at a single point). Model the top right quadrant of an elementary cubic cell, Fig. P4.2. Numerically solve the problem for different fluids and determine the permeability.
Quadrant used to model elemental volume y,v z,w
x,u
Figure P4.2
4.9. For most heat pipes, the wick structures are very thin, so the liquid flow in the wick can be simplified to one-dimensional axial flow. Since the liquid velocity and its gradient in the axial direction are very small, the liquid K dpA , where K is flow in the wick satisfies Darcy’s law, i.e., wA = − μA dz permeability. For an axially-grooved wick structure as shown in Fig. P4.3, determine the permeability analytically. 4.10. The designer of a coffee maker wants to make a filtration system in which the water is heated and filtered at the same time. The filter can be modeled as a porous zone within a tube, and the heat can be modeled as a uniform heat flux applied to the tube’s outer wall. Set up the governing equations and boundary conditions so that the designer can get a temperature and pressure field, and thus know where the liquid may start to boil. 4.11. A conventional heat pipe with multiple heat sources, as shown in Fig. P4.4, is divided into three regions: vapor-flow region, liquid-wick region, and solid wall. The effect of gravity can be neglected so that the fluid flow and heat transfer can be treated as two-dimensional. Specify governing equations and corresponding boundary conditions for all three regions.
322 Transport Phenomena in Multiphase Systems
Figure P4.3.
Figure P4.4.
4.12. The boundary of a porous zone is subjected to a constant heat flux. The porous zone is in local thermal nonequilibrium. Specify the thermal boundary conditions for the fluid phase and the thermal boundary condition for the solid phase. Also discuss the physical significance of your boundary conditions. 4.13. Consider flow entering from different sides of a tee-junction in a duct. The fluid entering is a mixture of air and liquid water, with known velocities and volume fraction at each inlet. The outlet is held at a constant pressure. A schematic of this problem is presented in Fig. P4.5. Model the problem using a two-dimensional Cartesian homogeneous model to find the mass averaged velocity. Assume that the flow is in the laminar regime and that gravitational effects are important. Also, assume that the flow field is isothermal. The volume fraction of air is small; therefore, the air exists in the form of bubbles. The diameter of these bubbles is known from experimental observation. 4.14. A parametric study on the effect of the inlet velocities and gas volume fraction from Problem 4.13 is needed. However, before performing this study using the homogeneous model, the inaccuracies of the model must first be quantified. The inaccuracies should be quantified by comparing the homogeneous results to a volume-averaged multifluid model. Therefore, set up the governing equation using a volume-averaged multifluid model with the corresponding boundary conditions with same assumptions used in Problem 4.13.
Chapter 4 Generalized Governing Equations: Averaging Formulations 323
Inlet 2: Liquid water/air mixture entering velocity of each component known
g
y, v x, u Outlet: Constant pressure
Inlet 1: Liquid water/air mixture entering velocity of each component known
Figure P4.5
4.15. There are many processes where large pressure drops in a liquid occur. If the pressure is low enough at one side, the liquid may locally vaporize. This process is called cavitation. Cavitation may occur in a liquid flowing past a sharp-edged circular orifice. The schematic of this problem is presented in Fig. P4.6. Liquid enters the orifice from the left. At the orifice there is a sharp and large pressure drop, causing vapor bubbles to form. The walls of the orifice are well insulated and can be considered adiabatic. The flow is inherently turbulent; therefore, use the k − ε model to model turbulent kinetic energy and dissipation. Develop the two-dimensional steady state governing equations using a homogeneous volume-averaged model. Assume that the vapor and liquid have the same velocities, therefore the mass-averaged velocity of the mixture is identically equal to each intrinsic phase average velocity.
324 Transport Phenomena in Multiphase Systems
Liquid In r, v x, u pin Axis
Adiabatic Wall
Liquid and Vapor Out
pout> pin
Figure P4.6
4.16. Consider the Direct Methanol Fuel Cell (DMFC) with a liquid feed presented in Fig. P4.7. Assume that the methanol solution is passing over the fuel cell at a very high velocity, and that solution has a constant concentration at the top of the anode gas diffusion layer. The oxygen supplied to the cathode comes from air. Assume that the catalyst layers are very thin, therefore their thickness can be neglected and the reaction rates can be considered to be a surface reaction rate. The membrane can be considered saturated at all times. Develop the governing equations and boundary conditions for the transport of liquid and gas through the fuel cell as well as the electric potentials of the membrane and gas diffusion layers for a one-dimensional steady state simulation.
Methanol Solution
x Catalyst Layers
Anode Gas Diffusion Layer Membrane Cathode Gas Diffusion Layer Air
Figure P4.7
Current Collectors
4.17. The multiphase transport model in porous media can also be used to model two-phase flow in microchannel heat exchangers. The transport phenomena in a microchannel heat exchanger can be characterized as unidirectional flow combined with three-dimensional conduction. Set up
Chapter 4 Generalized Governing Equations: Averaging Formulations 325
the governing equations for this system, while leaving the terms caused by phase interaction and solid/fluid interaction unclosed. The energy equation shall be written in terms of internal energy, and the solid, liquid and vapor equations shall be set up separately. 4.18. Consider a case of annular flow in which the liquid flows as a film on the wall of a pipe while the gas flows down the central cylindrical core. For a pipe of diameter D let the interfacial and wall shear stresses be τ δ and τ w as given inputs, respectively. Assuming vertical co-current flow configuration as shown in Fig. P4.8, develop the formulation for onedimensional, unsteady, laminar flow with multifluid modeling.
z
L
Gas
Liquid layer
r
Figure P4.8
4.19. In a fluidized bed particles are supported by an upward flow of fluid around them as shown in Fig. P4.9. Formulate the problem as multifluid model and deduce the fluid flux and pressure gradient needed to cause fluidization in a bed with void fraction ε f by neglecting inertia. Assume one-dimensional, incompressible, steady flow with negligible interparticle forces and no phase change. 4.20. Develop the general momentum and continuity equations for flow in a duct with phase change in terms of velocity and quality. The duct is inclined at an angle to the horizontal (Fig. P4.10). Assume laminar, incompressible, one-dimensional and steady flow. Velocities of the two phases are specified at the inlet of the duct.
326 Transport Phenomena in Multiphase Systems
z
z
θ
Fluid flow
Figure P4.9
Figure P4.10
4.21. Liquid fuel spray is coming out of an injector nozzle in an Ar-H2 plasma at a high temperature (see Fig P4.11). Develop the governing equations for mass, momentum, and energy based on a mixed Eulerian-Lagrangian technique, with Eulerian for gas phase and Lagrangian for the liquid. Your model should capture the essential features of primary atomization, using models that can obtain interactions between the fluid and droplets that can be directly applied to the standard atomization models used in practice. The problem can be assumed to be laminar, two-dimensional, unsteady, and with gravity neglected.
x
Figure P4.11 (Basu, 2005)
4.22. Deduce the general one-dimensional governing equations for mass and momentum for stratified turbulent flow of two phases (1) and (2) with evaporation over a flat plate as shown in Fig. P4.12.
Chapter 4 Generalized Governing Equations: Averaging Formulations 327
4.23. Formulate the analysis for heat transfer for Problem 4.22 using drift flux model. Assume steady, incompressible flow with negligible viscous dissipation. Use the results obtained in Problem 4.22 in carrying out the analysis.
z
Steam Flow (1) Liquid Film (2)
Figure P4.12
4.24. Develop the general expressions corresponding to pool boiling for maximum heat flux using drift flux method. Indicate the physical idea behind pool boiling. Boiling in this problem is a process when a pool of liquid is heated from below through a horizontal surface as shown in Fig. P4.13. A uniform bulk temperature is maintained far from the wall, meaning that most of the liquid is at a uniform temperature except near the wall, where there is a temperature difference in a thin layer. As shown in the figure, the dominant component of qw is simply the enthalpy difference between the upward and downward fluxes of vapor and liquid respectively. Also, find the maximum heat flux expression.
Vapor bubble
qw
Figure P4.13
Tw
4.25. Formulate steady flow in a one-dimensional nozzle under compressible inviscid conditions. The problem should be modeled as a homogeneous two-component flow with no phase change. Neglect the effects of surface tension. Obtain velocity and volume fraction distribution along the nozzle. 4.26. Derive the equations for mass, momentum, and energy for the flow of compressible gas through a channel containing dispersed particles. The flow occurs at high speed. Neglect the following effects: exchange of momentum, heat and mass with the channel walls, phase change, chemical reactions, viscous drag and viscous dissipation and gravity. Assume the
328 Transport Phenomena in Multiphase Systems
flow is one-dimensional, steady, and constant particle phase density. Do the analysis using the area-averaged multifluid models for channel flows. 4.27. Electro-osmosis is flow induced by an external electric field. When surfaces acquire a finite charge density and are in contact with a polar solution, they induce a distribution of electrical charges within the electrically charged solution. Electro-osmosis is used in biochip technology as a pumping mechanism in MEMS devices for chemical and biological analysis and diagnoses. Two-fluid electro-osmotic flow is used in many microfluids. The working fluid needs to be a conducting liquid with significant electrical conductivity to function properly. As shown in Fig. P4.14, a high electro-osmotic mobility liquid is at the bottom of the channel, and a low electro-osmotic mobility liquid is in the upper section of the microchannel. When an electric field is applied across the channel, pressure and electro-osmosis effects drive the liquids. The flow depends on the viscosity ratio of the two fluids, the external electric field, electroosmosis characteristics of the high mobility fluid, and the interfacial curvature between the two fluids. Develop the steady formulation for velocity and the electric potential due to an electrically applied field caused by charges on the wall.
y
Liquid 2 u2,e Interface Liquid 1 u1,e
e 0
H
x Ex
Figure P4.14
4.28. Develop the governing equations and boundary conditions for fluid mechanics, heat transfer, and volume of fluid (VOF) for a controlled liquid impinging jet on a rotating disk. Include the conjugate effects with and without evaporation as shown in Fig. P4.15 as a free surface flow. Fig. P4.16 shows a computational domain for the three regions: disk, liquid, and vapor, with an extended domain to include the edge effects. The basic assumptions are constant properties, incompressible and laminar
Chapter 4 Generalized Governing Equations: Averaging Formulations 329
flow. For evaporation, neglect the effect of mass transfer. constant heat flux along the heated section.
Assume
Fig. P4.15
Computational Domain
Extended Domain for Edge Effects
u x Water
in Collar
in
v
g
Vapor
Free surface 100 Liquid
in
r rin
Adiabatic
Disk
Heater (Constant Heat Flux) Adiabatic
dd
rhin rhout rout
axis
Fig. P4.16
330 Transport Phenomena in Multiphase Systems