4. Interfacial Area Density Transport Equation

☞ A method needs to predict the interfacial area density accurately even in case

of the propagation of the inlet condition and for the flow which varies timely

and spatially.

☞ The dependancy on flow regime can be reduced or negligible.

1) Method to solve the number density by transport equation and to predict the

interfacial area density in consideration of the geometry

- good for diffusion flow, or the flow including boiling in particular

- Many assumptions and simplification are accompanied  in applying to all

the flow regime

2) Method to solve the interfacial area density transport equation directly

① Multi-group model ; Consideration of bubble set

② One-group model

- The source term requires the empirical or theoretical model, and especially          many terms need the information on the distribution of the bubble size.

- difficulity in measurement and phenomenological complexity

(1) One-group Number Density Transport Equation

- basic macroscopic parameters of two-phase diffusive flow ; void fraction,

liquid film, number density, interfacial area density, bubble shape

- The number density is very important to determine the interfacial area

density.

- Correlation for the interfacial area density considering the bubble shape

[eqn: rm a_i = (36 pi )^1/3 `(r_s /r_v )^2 `N_b^1/3 `alpha^2/3]

(3-)

where rs and rv are length dimensions corresponding the surface area and

volume of bubble, which are same values if the bubble is sphere.

One-dimensional bubble number density transport equation by Ishii and            Kocamustafaogullari(1983)

[eqn: rm {partial N_b }over{partial t}+{partial}over{partial x} (N_b vec v `_g ) = phi_bn +phi_wn +phi_dis -phi_cond]

(3-)

① [eqn: phi_bn] ; bubble nucleation rate in flow by boiling

- The bubble nucleation is composed of homogeneous and heterogeneous

shapes.

[eqn: rm phi_bn =phi_"homo" +phi_hetero]

(3-3)

- Actuaaly [eqn: phi_hetero ~>~ phi_"homo"] ,  so the bubble nucleation by the homogeneous

boiling can be negligible.

In particular, it is remarkable in sharp depression that the two-phase

mixture is released through the short nozzle or orifice at high pressure.

② [eqn: phi_{w n}]  ;  bubble nucleation rate at the wall by boiling

∝ nucleation site density, Ncn, and bubble departure frequency, f

[eqn: rm phi_wn =(N_cn f xi_h )/A_c]

(3-4)

where ξh , Ac : heated perimeter and flow area

③ [eqn: phi_dis``,~ phi_cond]  ;  bubble nucleation rate by breakup, bubble sink rate by

coalescence and condensation

- bubble sink ; coalescence of bubbles & condensation at subcooled condition

- neglect coalescence until the transition to slug or churn turbulent flow

( α < 0.3 )

Then, Eq. (3-2) will be rearranged as,

[eqn: rm {partial N_b }over{partial t}+{partial}over{partial x} (N_b vec v `_g ) = LEFT ( {T_f -T_0}over{T_sat -T_0} RIGHT ) LEFT ( {N_cn f xi_h }over{A_c} RIGHT )+phi_hetero +phi_dis]

(3-5)

At first, consider the transport equation of bubble number density distribution

function to observe the physical phenomena related the bubble nucleation and

sink systematically.

[eqn: rm {partial f(v)}over{partial t}+del cdot(f(v) {vec v}_p )= Sum_j=1^4 S_j +S_ph]

(3-6)

where Sph is the source term by the phase change and the first term of the

right side is the change rate of bubble size distribution function by the

breakup and coalescence of bubble. S1 and S2 are related to bubble nucleation

and sink with volume v by breakup and S3 and S4 are related to bubble

nucleation and sink with volume v by coalescence. The related equations are

as follows.

[eqn: rm S _{1} ( vec{x} ,v,t)= int _{v} ^{v _{max}} {} beta (v prime ,v)n(v prime )g(v prime )f( vec{x} ,v prime ,t)dv prime][eqn: rm ~~~~~~~~S _{2} ( vec{x} ,v,t)=-g(v)f( vec{x} ,v,t)~~~~~~~~~~~~~~~~~~~~][eqn: rm S_3 (vec x ,v, t)=int_{v_min}^v/2 &lambda (v-v' ,v' )``h (v-v' ,v' )#&rm times f(vec x ,v-v' ,t) f(vec x ,v' ,t)dv' ~~]

(3-7)

(3-8)

(3-9)

[eqn: rm S_4 (vec x ,v, t)=-int_{v_min}^{v_max -v} &lambda (v ,v' )``h (v ,v' )#&rm times f(vec x ,v' ,t) f(vec x ,v ,t)dv' ~]

(3-10)

where  g(v) is the breakup frequency of bubble Volume v

β(v',v) : bubble size distribution of generated bubbles after breakup of

bubble Volume v'

n(v)	: number of generated bubbles after breakup of bubble Volume v

λ(v',v) : probability of coalescence between the collided bubbles Volume v

and v' ( collision efficiency )

h(v',v) : collision frequency between bubbles of Volume v and v'

Integrating over all the bubble size,

[eqn: rm {partial N}over{partial t}+del cdot(N {vec v}_pm )= Sum_j=1^4 int_{v_min}^{v_max}S_j ``dv+int_{v_min}^{v_max}S_ph ``dv]

(3-11)

where N(x,t) ; number density in unit volume of mixture

[eqn: rm N( vec x ,t) == int_{v_min}^{v_max} f(vec x ,v,t)dv]

(3-12)

vpm(x,t) ; average bubble velocity at specific point

[eqn: rm {vec v}_pm (vec x ,t)== {int_{v_min}^{v_max} f(vec x ,v,t) {vec v}_p (vec x ,v,t) dv}over{int_{v_min}^{v_max} f(vec x ,v,t) dv}]

(3-13)

In addition, Millies, Drew, Lahey(1996) and Kalkach(1994) obtained the numerical     and analytic solutions of bubble size distribution function using Boltzmann          transport equation. The balance equation of bubble size distribution function is      as follows.

[eqn: rm { partial f(v)} over {partial t } + del cdot [f(v)u_g ] =& int_v^inf b(u)f(u)du -int_0^v b(v) u over v f(v)du # rm & + 1 over 2 int_0^v c(u,v-u)f(u)f(v-u)du # rm &-int_0^inf c(u,v)f(u)f(v)du]

(3-14)

where f(v) is bubble number density with volume v.

Assuming the isotropic turbulence, Prince와 Blanch suggested

[eqn: rm b(u)=b_0]

(3-15)

[eqn: rm c(u,v)=c_0]

(3-16)

where

[eqn: rm b_0 =b_1 epsilon _L^9/5 LEFT ( sigma over rho_l RIGHT ) ^-11/5]

(3-17)

[eqn: rm c_0 =c_1 epsilon _L^-3/5 LEFT ( sigma over rho_l RIGHT ) ^7/5]

(3-18)

where [eqn: epsilon_L] is  local turbulent energy diffusion rate of liquid phase and b1, c1 are    the empirical coefficients(=0.0163 and 0.0369). Figure 3-1 shows the results from    numerical analysis on one-dimensional bubble flow using finite differential          method.

In addition, the analytic solution of Eq. (3-14) can be derived on the assumption    of the fully developed flow.

[eqn: rm f(v)= b_0 over c_0 exp LEFT [ - root{b_0 /(c_0 alpha_b )v} RIGHT ]]

(3-18)

The local interfacial area density is derived using Eqs. (3-16), (3-17), (3-18) and    the assumption of sphere bubble as given in Eq. (3-19).

[eqn: rm a_i &=(36pi)^1/3 int_0^inf u^2/3 f(u)du## rm &=5.41 epsilon _L^0.4 LEFT ( sigma over rho_l RIGHT ) ^-0.6 alpha_b^0.83]

(3-19)

Fig 3-. Bubble size distribution

(Millies et al., 1996)

(2) Interfacial Area Transport Equation

- It is convenient to use the number density equation in case of the diffusive        flow such as bubbly flow and mist flow because the source term can be           considered simply.

- However, it is difficult to apply the number density transport equation to all       the flow regime.

- The interfacial area density can be predicted directly by change the number        density transport equation into area density transport equation

① local approach : predict bubble  size distribution by modeling to break-up                           and coalescence processes of bubbles.

1) breakup process : based on the probability of liguird eddies  colliding with                            and destroying the individual bubble.

2) coalescence process : based on the probability that a bubble pair stays                              together long enough for sufficient drainage of the                                liquid film between the bubbles.

② global approach : obtain IAC on a global scale.

1) Ishii : as function of <α> and flow regime

2) Nagel et al. : applied Kolmogoroff's turbulence theory and as function of                           global void fraction and dissipation rate.

- One Group Model in VIPRE-02 code

; express the interfacial area density, momentum and energy through the             interface by using the dynamic flow regime map which presents the               structure of flow

The new variable G, called to be geometrical intensity is introduced as follows.

[eqn: rm G= a_i over alpha]

(3-20)

The interfacial density equation has the following one-dimensional form.

[eqn: rm {partial a_i}over{partial t}+{partial v_g a_i}over{partial x}= LEFT [ da_i over dt RIGHT ]_process#]

(3-21)

The right side of above equation includes the change in bubble density, the        bubble nucleation at the wall, bubble growth in flow, coalescence & breakup        reaction.  Eq. (3-21) can be modified as a form of geometrical intensity using      bubble steady-state continuity equation.

[eqn: rm {partial (a_i / alpha_g )}over{partial t}+v_g{partial (a_i / alpha_g )}over{partial x}&= LEFT [ da_i over dt RIGHT ]_process -{(a_i / alpha_g )} over {alpha_g rho_g}Gamma_g +{(a_i / alpha_g )}over rho_g LEFT [ {partial rho_g }over {partial t } +v_g{partial rho_g}over{partial x}RIGHT ]]

(3-22)

Change in bubble density

[eqn: rm - alpha_g over rho_g dA_b over dV_b LEFT [ {partial rho_g }over {partial t } + v_g{partial rho_g}over{partial x}RIGHT ]]

(3-23)

New bubble nucleation at the wall

[eqn: rm G_0 over rho_g Gamma_w]

(3-24)

Bubble growth in flow

[eqn: rm 1 over rho_g LEFT [ dA_b over dV_b RIGHT ] Gamma_b]

(3-25)

Bubble coalescence

[eqn: rm -f_c LEFT [ a_i over alpha_g - dA_b over dV_b RIGHT ]]

(3-26)

Bubble breakup

[eqn: rm f_b (a_i -A_eq )]

(3-27)

where Γw and Γb are the bubble nucleation rate at the wall and in the flow,        respectively.

Substituting above terms, Eq. (3-22) can be expressed by the equation of          the geometrical intensity, G as shown in Eq. (3-28).

[eqn: rm {partial G}over{partial t}+v_g{partial G}over{partial x}&= 1 over rho_g (G- dA_b over dV_b )({partial rho_g }over {partial t } +v_g{partial rho_g}over{partial x})# rm&+ 1 over {alpha_g rho_g} (G_0 -G)Gamma_w + 1 over {alpha_g rho_g} ( dA_b over dV_b -G) Gamma_b# rm &-f_c (G- dA_b over dV_b )+f_b (G- G_eq )]

(3-28)

The constitutive relations about the geometrical factors such as fc, fb, Geq,          dAb/dVb are needed to solve this equation. fc is a factor related to collision         for bubble coalescence as shown in Eq. (3-29).

[eqn: rm f_c = 1 over 2 {F_cw U_w}over{triangle X_b}]

(3-29)

where Uw is bubble wake velocity and the rest of factors are as follows.

[eqn: rm U_w ={U_inf e^{-(triangle r_b /R_b )^2}}over{A+B(triangle X_b /R_b )+C(triangle X_b /R_b )^2}]

[eqn: rm triangle r_b = min LEFT ( F_r triangle X_b , 1 over 2 D_e -R_b RIGHT )~~~~~~~~~]

[eqn: rm triangle X_b = V_b over A_b LEFT [ {alpha_c -alpha_g}over{alpha_c alpha_g} RIGHT ] ~~~~~~~~~~~~~~~~~~~]

(3-30)

Table 3-1 summarizes the correlations for these geometrical factors. This          model assumes that  the aspect ratio of bubble, ξ is 3.0 and the maximum         bubble radius is limited by the structure of flow path as shown in Eq. (3-31).

[eqn: rm r_max = 1 over 2 D_h f_a]

(3-31)

The non-dimensional variables about geometrical intensity and bubble radius       are defined as following.

[eqn: rm G^* =G~r_max]

(3-32)

[eqn: rm r^* = R_b over r_max]

(3-33)

Consequently, the actual bubble radius, Rb can be obtained using Table 3-1        and above correlations.

Table 3-. Correlation for dynamic flow regime model (VIPRE-02)

Factor	[eqn: rm G^* ge 3 over 2 (2+ xi )]	[eqn: rm G^* < 3 over 2 (2+ xi )]
[eqn: rm dA_b over dV_b]	[eqn: rm {4+5xi `r^* +2 xi^2 (r^* )^2 }over{2`r^* r_max LEFT ( 1+ 2 over 3 xi `r^* RIGHT ) }]	[eqn: rm 2 over r_max]
[eqn: rm V_b]	[eqn: rm {4pi over 3 (r^* )^3 ``r_max^3 }over{1+xi`r^* }]	[eqn: rm {4pi over 3 r_max^3 }over{1+xi}+L_b^* pi ``r_max^3]
[eqn: rm A_b]	[eqn: rm LEFT [ {2(r^* )^2 (2+xi r^* )}over{1+xi `r^*} RIGHT ] pi]	[eqn: rm LEFT [ {2 (2+xi )}over{3+xi }+L_b^* RIGHT ] pi `r_max^2]
[eqn: rm A_F over a_i]	[eqn: rm {1+xi`` r^* }over{2+xi ``r^*}]	[eqn: rm 1 over {{2+xi}over{1+xi}+L_b^*}]
[eqn: rm L_b]	0	[eqn: rm LEFT [ {2 over 3 (2+xi) -G^*}over {3 (G^* -2 )} RIGHT ]~ r_max]
[eqn: rm r^*]	[eqn: rm 3 over {G^* - 2 over 3 xi}]	1

- Kocamustafaogullari and Ishii(1995) multiplied average interfacial area by           the bubble number density transport equation.

[eqn: rm {partial` f(v) A_i (v)}over{partial t}+del cdot(f(v) A_i (v) {vec v}_p )= Sum_j=1^4 `` A_i (v)S_j +`` A_i (v)S_ph]

(3-34)

Integrating this equation over all bubble sizes,

[eqn: rm {partial` a_i }over{partial t}+del cdot( a_i {vec v}_i )= Sum_j=1^4 phi_j +phi_ph]

(3-35)

where

[eqn: rm a_i ( vec x ,t) == int_{v_min}^{v_max} f(vec x ,v,t) A_i (v)dv ~~~~~~~~~~~~]

[eqn: rm {vec v}_i (vec x ,t)== {int_{v_min}^{v_max} f(vec x ,v,t) A_i (v){vec v}_p (vec x ,v,t) dv}over{int_{v_min}^{v_max} f(vec x ,v,t) A_i (v)dv}]

and øj is the change rate of the interfacial area by bubble breakup and             coalescence. øph is the change rate of the interfacial area by the phase             change.

The equations can be rewritten using Si, Sph which is defined in Eqs. (3-7)～       (3-10).

[eqn: rm phi_1 = int_{v_min}^{v_max} &A_i (v) int_v^v_max beta(v' ,v) n(v' ) #& times rm g(v' )`f(vec x ,v' ,t)dv' dv ~~~~~~~~~~]

[eqn: rm phi_2 = -int_{v_min}^{v_max} A_i (v) g(v )f(vec x ,v ,t) dv ~~~~~~~~~~~]

[eqn: rm phi_3 =int_{v_min}^{v_max} &A_i (v) int_{v_min}^v/2 lambda (v-v' ,v' )``h (v-v' ,v' )#&rm times f(vec x ,v-v' ,t) f(vec x ,v' ,t)dv' dv]

[eqn: rm phi_4 =-int_{v_min}^{v_max} &A_i (v) int_{v_min}^{v_max} lambda (v ,v' )``h (v ,v' )#&rm times f(vec x ,v' ,t) f(vec x ,v ,t)dv' dv ~~~~~]

[eqn: rm ~~~~~~~ phi _{ph} = int _{v _{min}} ^{v _{max}} {} A _{i} (v)S _{ph} dv~~~~~~~~~~~~~~~~~~~~~~~~]

Also, the following transport equation of voild fraction can be obtained by          integrating the number density distribution function multiplied by volume.

[eqn: rm {partial` alpha }over{partial t}+del cdot( alpha {vec v} )=int_{v_min}^{v_max} LEFT ( Sum_j=1^4 v S_j +S_ph v RIGHT ) dv]

(3-36)

In case of incompressible gas without phase change,

[eqn: rm {partial` alpha }over{partial t}+del cdot( alpha {vec v}_g )=0]

(3-37)

- Wu and Ishii(1997) derived one-group one-dimensional average interfacial           area transport equation at the system without phase change.

[eqn: rm { partial <a_i >} over {partial t } + partial over {partial z } LEFT ( <a_i >{bar v}_iz RIGHT ) =S_a,rc +S_a,ti +S_a,we - LEFT ( 2a_i over {3alpha rho_g} RIGHT ) {D_g rho_g}over Dt ~~~~]

(3-38)

The last term of the right side means bubble expansion and contraction by the      change in pressure and this term is 0 on the assumption that the bubble is        incompressible.

The interfacial velocity can be expressed as shown in Eq. (3-39) on the           assumption that the bubble has a uniform size,

[eqn: rm {bar v}_iz =={ LEFT < a_i v_iz RIGHT >} over {LEFT < a_i RIGHT > } = { LEFT < alpha v_iz RIGHT > } over { LEFT < alpha RIGHT > } = LEFT < LEFT < v_iz RIGHT > RIGHT >]

(3-39)

The interfacial velocity is same with the vapor velocity if the flow in bubble       is ignored. Thus,

[eqn: rm {bar v}_iz APPROX LEFT < LEFT < v_g RIGHT > RIGHT > = { LEFT < j_g RIGHT > } over { LEFT < alpha RIGHT > }]

(3-40)

In addition, the transport equation of Sauter average diameter for the              incompressible flow without phase change is as follows.

[eqn: rm { partial D } over {partial t } + LEFT < LEFT < v_g RIGHT > RIGHT > { partial D} over {partial z } =S_D,rc +S_D,ti +S_D,we]

(3-41)

Then, from channel-averaged vapor velocity and void fraction,

[eqn: rm a_i (x,t) = 6alpha(x,t) over{ D_sm (x,t)}]

(3-42)

Table . Experimental Case (Wu et al., 1997)

	[eqn: rm LEFT < j_g RIGHT >] (m/s)	[eqn: rm LEFT < j_f RIGHT >] (m/s)	[eqn: rm LEFT < alpha RIGHT >]	D(mm) (L/D=2)
Case 1	0.023	0.77	0.019	3.40
Case 4	0.117	1.58	0.061	4.69

Fig . Bubble Size (case 1)

(Wu et al., 1997)

Fig . Bubble Size (case 4)

(Wu et al. 1997)

Fig . Interfacial Area Density (case 4)

(Wu et al., 1997)

Fig . Interfacial Area Density (case 4)

(Wu et al., 1997)