2. One-Dimensional Two-Phase Flow Model

Mass Continuity

[eqn: {Rate~of~Creation} choose {of~Mass} ~ = ~{Mass~Outflow} choose Rate ~ - ~ {Mass~Inflow} choose Rate ~+~ {Mass Storage} choose Rate]

(ROC)        =       [eqn: dot O]        -       [eqn: dot I]        +       [eqn: dot S]

[eqn: A partialrho over partialt + partialw over partialz = 0] ,       where [eqn: w=rho u A]

For each phase,

[eqn: partial over partialt (A``rho_k`` alpha_k ) + partial over partialz (A``alpha_k `` rho_k ``u_k )~=~ GAMMA_k~,~~~~~~~~ SUM from { k }`` GAMMA_k ``=`` 0]

Momentum Equation

[eqn: {ROC ~~of} choose Momentum ~=~ {Momentum} choose {Outflow~ Rate} - {Momentum} choose {Inflow~ Rate} + {Momentum} choose { Storage~ Rate}]

[eqn: =~{S um ~of~ the~ Forces} choose {acting~ on~ the~ C.V}][eqn: ( Sum`` F)]

[eqn: therefore ~~ Sum`` F ~=~ TRIANGLE p_static ~+~ TRIANGLE p_fric ~+~ TRIANGLE p_accel ~+~ TRIANGLE p_grav ~+~ TRIANGLE p_local]

For vapor phase,

[eqn: -A_g dp -dF_g -S - A_g`` dz`` rho_g`` g ``sin theta ~=~ w_g`` du_g + dw_g`` u_g - dw_g``u_f + partialw_g over partialt A_g] (4.2.13)

For liquid phase,

[eqn: -A _{f} dp-dF _{f} +S-A _{f} ``dz`` rho _{f} ``g``sin theta `=~w _{f} ``du _{f} + {Partial w _{f}} over {Partial t} A _{f}]             (4.2.14)

Adding Eqs. (4.2.13) & (4.2.14) and [eqn: dw_g ~=~-dw_f]

[eqn: -A dp - dF_g - dF_f - g ``sin theta``dz`` (A_f rho_f + A_g rho_g )]

[eqn: ~=~ d(w_f`` u_f + w_g`` u_g ) + partial over partialt (w_g`` A_g + w_f`` A_f )]               (4.2.15)

Here, interfacial shear stress,  [eqn: S ~=~ tau_i``P_i`` TRIANGLE z]

where   [eqn: tau_i ~=~ (C_f )_i``` rho_g over 2g_c (u_g - u_l )^2]       with no entrainment

[eqn: (C_f )_i ~=~ 0.005[1+75(1-<alpha>)]]  ; interfacial friction factor

And  [eqn: dw_f ~=~ -dw_g ~=~ delta w']

[eqn: =~ {{dot q}'' P_H} over h_fg]   ; evaporation amount per unit axial length

★ See Lahey's Table 5-II

Energy Equation

[eqn: {Rate~of~Creation} choose {of~Energy} ~ = ~{Energy} choose {Outflow~ Rate} ~ - ~ {Energy} choose {Inflow~Rate} ~+~ {Energy} choose {Storage~Rate}]

[eqn: w``(``` delta Q - delta W````) ~=~ w``d h + d`` left[ {w_g`` u_g^2} over 2 + {w_f`` u_f^2} over 2 right] + w``g``sin theta``dz]             (4.2.16)

(1) Homogeneous Model

○ [eqn: u_g ~=~ u_l ~=~ u]

○ thermodynamic equilibrium between phases

○ use suitably-defined  f for wo-phase flow

★ valid for

① bubbly and wispy-annular flow pattern

② high velocity flow in large channel

[eqn: <rho_H > ~=~ 1/(v_f +<x>v_fg )]                                     (4.2.17)

[eqn: <h > ~=~ h_f +<x>h_fg]                                         (4.2.18)

[eqn: G ~=~ <rho_H ><j>``,]  and <j> = u for homogeneous flow         (4.2.19)

For constant flow area,

1) mass continuity

[eqn: partialrho over partialt ~+~ partialG over partialz ~=~0]                                                (4.2.20)

2) momentum equation

[eqn: {<rho_H >} over g_c Du over Dt ~=~ - partialp over partialz ~-~ g over g_c <rho_H > sin theta ~-~ {tau_w P_f} over A]                   (4.2.21)

3) energy equation

[eqn: <rho_H > D<e> over Dt ~=~ {dot q}'' left( P_H over A right)]   where [eqn: <e> ~=~ left( <h> + u^2 over {2g_c J} + g over g_c {z sin theta} over J right)] (4.2.22)

★ thermal energy equation  ;  Eq. (4.2.21) x u

[eqn: <rho_H > D<h> over Dt ~=~ {dot q}'' ``left( P_H over A right) ~+~ u ``left( {tau_w P_f} over JA right)]                         (4.2.23)

4) quality propagation equation

From Eqs. (4.2.18) and (4.2.19), for constant system pressure, if we neglect

△p

compared to the system pressure,

[eqn: D<x> over Dt ~-~ OMEGA`` <x> ~=~ Omega v_f over v_fg]                          (4.2.24)

where [eqn: Omega]  :  characteristic frequency of phase change

(즉 speed at which phase change take place)

[eqn: under TRIANGLE {dot q}'' `` P_H over A v_fg over h_fg]

(2) Separate Flow Model

○ [eqn: u_f ~neq~u_g]

○ thermodynamic equilibrium between phases

○ [eqn: phi^2 , ~alpha]

★ valid for annular flow

For constant flow area

1) continuity  equation

[eqn: {partial{< bar rho >}} over partialt ~+~ partialG over partialz ~=~0]                                    (4.2.25)

where   [eqn: < bar rho > ~=~ rho_f ``(1-<alpha>) ``+`` rho_g``<alpha>]

[eqn: G ~=~ rho_f ``(1-<alpha>)u_f ``+`` rho_g``<alpha>u_g]

2) momentum equation

○ liquid phase

[eqn: -(1-<alpha>) partialp over partialz - g over g_c rho_f``(1-<alpha>)sin theta - {tau_w P_f} over A - {tau_i P_i} over A]

[eqn: =~ 1 over g_c left[ partial over partialt [ rho_f ``(1-<alpha>) u_f ] + partial over partialz [ rho_f ``A (1-<alpha>) u_l^2 ] right] + deltaw' over gA u_i]    (4.2.26)

○ vapor phase

[eqn: -<alpha> partialp over partialz - g over g_c rho_g``<alpha>sin theta - {tau_i P_i} over A]

[eqn: =~ 1 over g_c left[ partial over partialt [ rho_g ``<alpha> u_g ] + partial over partialz [ rho_g ``A <alpha> u_g^2 ] right] - deltaw' over gA u_i]     (4.2.27)

○ two-phase mixture

[eqn: - partialp over partialz - g over g_c < bar rho >``sin theta - {tau_w P_f} over A ~=~ 1 over g_c left[ partialG over partialt + partial over partialz left( G^2 over <rho'> right) right]]       (4.2.28)

where  [eqn: 1 over <rho'> ~ under TRIANGLE ~ left[ (1-<x>)^2 over {rho_f`` (1-<alpha>)} + <x>^2 over {rho_g``<alpha>} right]]     ;  momentum density

3) Energy Equation

[eqn: rho_f``(1-<alpha>) De_f over Dt + rho_g ``<alpha> De_g over Dt ~=~ {{dot q}'' P_H} over A - deltaw' over A (e_g - e_f )]       (4.2.29)

Fig.4.2 Control volume for separated flow(each phase) [Lahey's Fig.5.3]

■  Interfacial exchange for two-fluid model

1) interfacial mass transfer

[eqn: GAMMA ~=~ deltaw over A ~=~ a_i`` m]                                         (4.2.30)

where  [eqn: a_i] ;  local interfacial area concentration  (IAC)

- the probability of the existence of interfacial area at local                             position in the two-phase flow

m ;  averaged mass transfer rate per unit area

2) interfacial energy transfer

[eqn: h_ki `` GAMMA_k ~+~ q_ki over L_s ~=~ a_i ``[m_ki ``h_ki ~+~ q_ki``]]                            (4.2.31)

(3) Drift-Flux Model

[Ref] G.B. Wallis, "One-dimensional Two-phase Flow" McGraw-Hill

★ Two-fluid model ; numerical complication and uncertainties in

specification of interfacial transport

In drift-flux model, concern the relative motion between phases rather than each phase motion

o mixture center of mass velocity

o drift velocity ; volumetric flux of each phase relative to a surface

Fig.4.2 Graphical solution

of jgf

moving at average velocity, j

The drift flux is defined in many ways as,

[eqn: V_gj ~(local~ drift ~flux)~=~ v_g - j],

[eqn: j_gf ~=~ alpha(v_g - j) ~=~ alpha V_gj]       in Wallis

[eqn: bar {(V_gj )} ~=~ {<alpha(v_g j)>} over <alpha>]    in Lahey & Zuber

[eqn: =~<<V_gj >>]            in Ishii

Some relations ;

[eqn: j_gf ~=~ u_r alpha (1-alpha)]              (4.2.32)

[eqn: j_g ~=~ alpha`` j + j_gf]

Thus, all properties ([eqn: alpha,~rho_m,], etc) are sum of     the homogeneous flow value and a correction    factor which is a function of drift flux or additional term.

[eqn: alpha ~=~ j_g over j (1- j_gf over j_g )] ,     [eqn: rho_m = {j_l`` rho_l + j_g`` rho_g } over j ~+~ (rho_l - rho_g ) j_gf over j]

[eqn: therefore]  If the flow is homogeneous,  [eqn: j_gf ~=~ 0]

☞   [eqn: j_g = alpha over {1-alpha} j_f + 1 over {1-alpha} j_gf]

Then, for a constant [eqn: alpha] , a plot of jg vs jl

will be a straight line of slope [eqn: alpha over {1-alpha}]

■ Modeling

o area averaging

o no information on change of variables normal to flow direction

o need empirical correlations or simplified models for the transfer of

momentum and energy between the wall and the fluid.

A.  Cross-Section Averaged Kinematic Equation

<F> = [eqn: 1 over A int from A F dA]

For constant system pressure, the vapor-phase continuity equation, will be

[eqn: {partial alpha} over {partial t} + partial over partialz [alpha(j+V_gj )] = GAMMA over rho_g]                                (4.2.33)

Taking [eqn: {1} over {A} int _{} ^{} {} [Eq.(4.2.33)]dA],

[eqn: {partial <alpha>} over {partial t} + partial over partialz [<alpha j>+<alpha V_gj >] = <GAMMA_g> over rho_g]

[eqn: therefore~~{partial <alpha>} over {partial t} + partial over partialz [C_o <alpha><j>+ {<alpha V_gj >} over <alpha> <alpha>] = <GAMMA_g> over rho_g]             (4.2.34)

where  [eqn: C_o ~ under TRIANGLE ~ <alphaj> over <alpha><j>]

[eqn: therefore~~{partial <alpha>} over {partial t} + C_o <j> partial<alpha> over partialz + C_o <alpha> partial<j> over partialz + bar{V_gj} partial<alpha> over partialz = <GAMMA_g> over rho_g]   (4.2.35)

Since   [eqn: partial<j> over partialz = OMEGA]  ,   where  [eqn: OMEGA = v_fg {delta w'} over A_x-s = {q_o '' P_H v_fg} over {A_x-s h_fg}]

[eqn: {partial <alpha>} over {partial t} + [C_o <j>+bar{V_gj} ] partial<alpha> over partialz = ( rho_f over TRIANGLE rho - C_o <alpha> ) OMEGA]          (4.2.36)

Now,  let's define   CK (kinematic velocity) = [eqn: C_o <j> + bar{V_gj}]

Then,

[eqn: {partial <alpha>} over {partial t} + C_K partial<alpha> over partialz = ( rho_f over TRIANGLE rho - C_o <alpha> ) OMEGA]

(4.2.37)

=> Void Propagation Equation

B. Drift-Flux parameters

① Co  ; concentration parameter

[eqn: C_o ~under TRIANGLE ~ 1 over <j><alpha>A`` int int from A (j alpha) dA]                                    (4.2.38)

[eqn: = ~ <alphaj> over <alpha><j>]     ; global slip due to cross-sectional averaging

From definition of average velocity of vapor phase,

[eqn: v_g ~=~ {<j_g >} over <alpha>]                                          (4.2.39)

Averaging Eq.(4.2.32),  [eqn: <j_gf ``> ~=~ <j_g > - <alpha j>]

[eqn: therefore~~{<j_g >} over <alpha> ~=~ C_o <j> + bar V_gj]                                 (4.2.40)

Or    [eqn: <alpha> ~=~ {<j_g >} over {[C_o <j> + bar{V_gj} ]}]

Note that, from Eq.(4.1.18),

[eqn: S`=~ under{{1-C _{o} < alpha >} over {1-< alpha >}} `+~ under{{V _{gj} (1-< alpha >)} over {(1-C _{o} < alpha >)<j _{f} >}}]                    (4.2.41)

integral slip         local slip

For homogeneous model ( [eqn: V_gj ~=~ 0] ), local slip will be vanished.

And,

[eqn: <alpha> ~=~ <x> over { C_o ``left[ <x> + rho_g over rho_f (1-<x>) right] + {rho_g V_gj } over G}]                   (4.2.42)

=> Zuber-Findley Void-Quality Correlation

Since [eqn: <j> ~=~ {Q_1 + Q_2 } over A ~,~~~<j_k > ~=~ Q_k over A]   and  [eqn: bar V_gj ~=~ {<j_21 >} over <alpha>] ,

Eq. (4.2.40) can be rewritten as,

[eqn: <alpha> ~=~ {Q_2 - A<j_21 > } over {C_o (Q_1 + Q_2 )}] [eqn: ～ ~ 1 over C_o Q_2 over {Q_1 + Q_2}]

1) Bankoff model

Neglecting the local slip with  [eqn: V_gj ~=~ 0] , Eq. (4.2.40) will be

[eqn: <alpha> over <beta> ~=~ 1 over C_o ~under TRIANGLE ~ kappa] ,     where κ = 0.71 + 0.0001 p (psia) for water            Here β, volumetric flow fraction, is defined as, [eqn: under TRIANGLE ~{<j_g >} over <j>]

☞ However, Vgj is important for the flow of low G.

☞ applicable to determine the two-phase level in steam generators

2) Dix model

; Co is function of <β> which is dependent on pressure and quality

[eqn: C_o ~=~ <beta> [1+(1/<beta> -1)^b ]]            where  [eqn: b ~ under TRIANGLE ~ (rho_g rho_f )^0.1]

If Co = 1,

[eqn: <beta> ~=~ <alpha> ~=~ <x> over {[<x> + rho_g over rho_f (1-<x> ) ]}]                 (4.2.43)

② Vgj

From fractional analysis for bubbly/slug flow, the important forces on a              bubble are follows;

net bouyancy force    ;  Fb  =  [eqn: 1 over 6 pi`` D_b^3 ``(rho_l -rho_g ``) g over g_c]

inertia force           ;  Fi  =  [eqn: rho_l over g_c U_t^2 ``( 1 over 4 pi`` D_b^2 ``)]

surface tension force   ;  Fσ  =  [eqn: pi``D_b ``sigma]

In the churn-turbulent flow, the most important force groups are

1) ratio Fi to Fb

2) ratio Fi to Fσ

Then,   [eqn: {rho_l`` U_t^2} over {D_b``g``(rho_l -rho_g )} ~ under TRIANGLE ~ k_1^2]        and      [eqn: {rho_l`` U_t^2 ``D_b } over {g_c``sigma} ~ under TRIANGLE ~ k_2^2]

The product of these dimensionless groups yields,

[eqn: {rho_l^2`` U_t^4} over {g``g_c``sigma``(rho_l -rho_g )} ~=~k_1^2``k_2^2 ~ under TRIANGLE ~ k_3^4]

Then, the terminal rise velocity will be,

[eqn: U_t ~=~ k_3 ``left[ {g``g_c``sigma``(rho_l -rho_g )} over rho_l^2 right]^1/4]

For churn-turbulent flow, [eqn: V_gj ~=~ U_t]

▶ For bulk boiling (Lahey, 1977),    k3 = 2.9

▶ For churn-turbulent flow (Zuber, 1967),    k3 = 1.41

▶ For slug flow (Zuber, 1969),    [eqn: V_gj ~=~ 0.35 left[ {g TRIANGLE rho D} over rho_f right]^1/2]

▶ For two-phase pseudo jet flow (Zuber, 1969)   [eqn: V_gj ~=~ constant`` left[ { sigma ``g TRIANGLE rho } over rho_g^2 right]^1/4]

▶ For the dispersed two-phase system in vertical motion

[eqn: V_gj ~=~ v_inf (1-alpha)^m]

Fig.4.5 Void concentration profile and flow regime[Lahey}

Fig.4.6 Vgj vs void fraction

Fig.4.4 Co vs void fraction

(4) Numerical Stability in Two-Phase Modeling

For two-fluid model, generally

[eqn: A`` partialU over partialt ~+~ B`` partialU over partialz ~=~ C]                                       (4.2.44)

where U = (α, p, ul, uv, hl, hg)T

Thus,

[eqn: left [ A right ]][eqn: partial over partialt][eqn: bmatrix { alpha`` # p`` # u_l`` # u_v`` # h_l`` # h_g`` }] + [eqn: left[B right] partialU over partialz ~=~ left[ S right]]     =>  [eqn: left| A lambda ~-~ B``` right| ~=~ 0~, ~~lambda_i ~=~ partial over partialt]

a) If all six roots are real and distinct  => Hyperbolic

b) If we have some same roots         => Parabolic

c) If we have complex conjugate pair   => Elliptic

1) [eqn: partialphi over partialt ~+~ nu partialphi over partialz ~=~ 0]                                                  (4.2.45)

Assume a perturbation in φ of the form,

[eqn: delta phi ~=~ delta phi_o`` e^{[i(kz-omegat)]}] ,   where k is real                           (4.2.46)

In order for this system to be well-posed, ω should be real.

Inserting Eq. (4.4.46) into Eq.(4.4.45),

[eqn: -i omega`` delta phi_o`` e^{[i(kz-omegat)]} ~+~ i`` k``nu`` delta phi_o`` e^{[i(kz-omegat)]} ~=~0]

∴ -iω + ikν = 0   or ω = kν

i.e The eigenvalue of this system, v, should be real.

For real [eqn: nu] , stable  when  Cv ([eqn: == {nu TRIANGLE t} over {TRIANGLE z}]) < 1

For complex eigenvalue on a finite difference scheme [Handout 참조]

[eqn: nu ~=~ nu_o (1+i ``epsilon)]

a) For ε < 0 , stable when Cv < 1

b) For ε > 0 , may be unstable for some n even when Cv < 1

2) Charecteristic Analysis and Time Discretization of the HEM Model

[eqn: A`` partialU over partialt ~+~ B`` partialU over partialz ~=~ 0]

where  U = [ρ, u, v, p]  , u ; internal energy

A = [eqn: bmatrix { 1 & 0 & 0 & 0 # 0 & rho & 0 & 0 # 0 & 0 & rho & 0 # 1 & 1/beta^2 & 0& -1/c^2 }] ,  B = [eqn: bmatrix { v & 0 & rho & 0 # 0 & rho v & p & 0 # 0 & 0 & rho v & 1 # 0 & 0 & 0 & 0 }] ,  [eqn: 1/beta^2 ~=~ - partialrho over partialu~ , ~~~1 over c^2 ~=~ partialrho over partialp]

Now, for explicit time discretization and forward differencing, we can perturb the

continuity and momentum equation

[eqn: rho_j^t ~=~ rho_o ~+~ delta rho_j^t]

[eqn: v_j^t ~=~ v_o ~+~ delta v_j^t]

Assume  [eqn: delta`` phi_j^t ~=~ delta``phi^t ``e^ikz_j]

Then, the linearized continuity and momentum equation will be

[eqn: U^t+Δt = GU^t]

where U = [δρ, δv]T    and  G = [eqn: bmatrix {{(1- bar v )} & {- rho {TRIANGLE t} over {TRIANGLE z} ik' } # {- c^2 over rho {TRIANGLE t} over {TRIANGLE z} ik' } & {(1- bar v )} }]

[eqn: bar v ~=~ {v_o TRIANGLE t } over { TRIANGLE z} [1-e^{ik TRIANGLE z}]] ,      [eqn: k' ~=~ 2``sin``(k TRIANGLE z /2)]

For

[eqn: left| G-I ``lambda right| ~=~ 0 ~,~~~~~ lambda_1,2 ~=~ bmatrix { {1- bar v + i C_n} # {1- bar v - i C_n} }]

where

Fig. 4.7 Example of spatial convergence of          void vs length

[eqn: C_n ~ under TRIANGLE ~ {C TRIANGLE t} over { TRIANGLE z} k']

For  v << Cn , λ1,2 > 1 ,  unstable

Thus, the implicit solution can be

recommanded.